UP | HOME

Atomic Orbitals

Table of Contents

1 Context

The following arrays are stored in the context:

type   Gaussian ('G') or Slater ('S')
shell_num   Number of shells
prim_num   Total number of primitives
shell_center [shell_num] Id of the nucleus on which each shell is centered
shell_ang_mom [shell_num] Angular momentum of each shell
shell_prim_num [shell_num] Number of primitives in each shell
shell_prim_index [shell_num] Address of the first primitive of each shell in the EXPONENT array
shell_factor [shell_num] Normalization factor for each shell
exponent [prim_num] Array of exponents
coefficient [prim_num] Array of coefficients

For H2 with the following basis set, 

HYDROGEN
S   5
1         3.387000E+01           6.068000E-03
2         5.095000E+00           4.530800E-02
3         1.159000E+00           2.028220E-01
4         3.258000E-01           5.039030E-01
5         1.027000E-01           3.834210E-01
S   1
1         3.258000E-01           1.000000E+00
S   1
1         1.027000E-01           1.000000E+00
P   1
1         1.407000E+00           1.000000E+00
P   1
1         3.880000E-01           1.000000E+00
D   1
1         1.057000E+00           1.0000000

we have: 

type = 'G'
shell_num = 12
prim_num = 20
shell_center = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
shell_ang_mom = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
shell_factor  = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]
shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
shell_prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
exponent = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
             1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
             0.3258, 0.1027, 1.407, 0.388, 1.057]
coefficient = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
                1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
                0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]

1.1 Data structure

typedef struct qmckl_ao_basis_struct {
  int32_t   provided;
  int32_t   uninitialized;
  int64_t   shell_num;
  int64_t   prim_num;
  int64_t * shell_center;
  char    * shell_ang_mom;
  int64_t * shell_prim_num;
  int64_t * shell_prim_index;
  double  * shell_factor;
  double  * exponent    ;
  double  * coefficient ;
  char      type;
} qmckl_ao_basis_struct;

The uninitialized integer contains one bit set to one for each initialization function which has not bee called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and provided == 1.

1.2 Access functions

Access to scalars copies the values at the passed address, and for array values a pointer to the array is returned. 

char      qmckl_get_ao_basis_type             (const qmckl_context context);
int64_t   qmckl_get_ao_basis_shell_num        (const qmckl_context context);
int64_t   qmckl_get_ao_basis_prim_num         (const qmckl_context context);
int64_t*  qmckl_get_ao_basis_shell_center     (const qmckl_context context);
char*     qmckl_get_ao_basis_shell_ang_mom    (const qmckl_context context);
int64_t*  qmckl_get_ao_basis_shell_prim_num   (const qmckl_context context);
int64_t*  qmckl_get_ao_basis_shell_prim_index (const qmckl_context context);
double*   qmckl_get_ao_basis_shell_factor     (const qmckl_context context);
double*   qmckl_get_ao_basis_exponent         (const qmckl_context context);
double*   qmckl_get_ao_basis_coefficient      (const qmckl_context context);

int32_t   qmckl_ao_basis_provided           (const qmckl_context context);

if ( (ctx->ao_basis.uninitialized & mask) != 0) {
  return NULL;
}

char qmckl_get_ao_basis_type (const qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (char) 0;
  }

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

  int32_t mask = 1;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return (char) 0;
  }

  assert (ctx->ao_basis.type != (char) 0);
  return ctx->ao_basis.type;
}


int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (int64_t) 0;
  }

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

  int32_t mask = 1 << 1;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return (int64_t) 0;
  }

  assert (ctx->ao_basis.shell_num > (int64_t) 0);
  return ctx->ao_basis.shell_num;
}


int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return (int64_t) 0;
  }

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

  int32_t mask = 1 << 2;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return (int64_t) 0;
  }

  assert (ctx->ao_basis.prim_num > (int64_t) 0);
  return ctx->ao_basis.prim_num;
}


int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 3;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.shell_center != NULL);
  return ctx->ao_basis.shell_center;
}


char* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 4;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.shell_ang_mom != NULL);
  return ctx->ao_basis.shell_ang_mom;
}


int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 5;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.shell_prim_num != NULL);
  return ctx->ao_basis.shell_prim_num;
}


int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 6;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.shell_prim_index != NULL);
  return ctx->ao_basis.shell_prim_index;
}


double* qmckl_get_ao_basis_shell_factor (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 7;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.shell_factor != NULL);
  return ctx->ao_basis.shell_factor;
}


double* qmckl_get_ao_basis_exponent (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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


  int32_t mask = 1 << 8;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.exponent != NULL);
  return ctx->ao_basis.exponent;
}


double* qmckl_get_ao_basis_coefficient (const qmckl_context context) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return NULL;
  }

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

  int32_t mask = 1 << 9;

  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    return NULL;
  }

  assert (ctx->ao_basis.coefficient != NULL);
  return ctx->ao_basis.coefficient;
}


int32_t qmckl_ao_basis_provided(const qmckl_context context) {

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

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

  return ctx->ao_basis.provided;
}

1.3 Initialization functions

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

qmckl_exit_code  qmckl_set_ao_basis_type             (qmckl_context context, const char      t);
qmckl_exit_code  qmckl_set_ao_basis_shell_num        (qmckl_context context, const int64_t   shell_num);
qmckl_exit_code  qmckl_set_ao_basis_prim_num         (qmckl_context context, const int64_t   prim_num);
qmckl_exit_code  qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index);
qmckl_exit_code  qmckl_set_ao_basis_shell_center     (qmckl_context context, const int64_t * shell_center);
qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom    (qmckl_context context, const char    * shell_ang_mom);
qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num   (qmckl_context context, const int64_t * shell_prim_num);
qmckl_exit_code  qmckl_set_ao_basis_shell_factor     (qmckl_context context, const double  * shell_factor);
qmckl_exit_code  qmckl_set_ao_basis_exponent         (qmckl_context context, const double  * exponent);
qmckl_exit_code  qmckl_set_ao_basis_coefficient      (qmckl_context context, const double  * coefficient);

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

qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

ctx->ao_basis.uninitialized &= ~mask; 

if (ctx->ao_basis.uninitialized == 0) {
  ctx->ao_basis.provided = 1;
}

return QMCKL_SUCCESS;

qmckl_exit_code qmckl_set_ao_basis_type(qmckl_context context, const char t) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  if (t != 'G' && t != 'S') {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_type",
                           NULL);
  }

  int32_t mask = 1;
  ctx->ao_basis.type = t;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code qmckl_set_ao_basis_shell_num(qmckl_context context, const int64_t shell_num) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  if (shell_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num <= 0");
  }

  int64_t prim_num = qmckl_get_ao_basis_prim_num(context);

  if (0L < prim_num && prim_num < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "shell_num > prim_num");
  }

  int32_t mask = 1 << 1;
  ctx->ao_basis.shell_num = shell_num;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_prim_num(qmckl_context context, const int64_t prim_num) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  if (prim_num <= 0) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "prim_num must be positive");
  }

  int64_t shell_num = qmckl_get_ao_basis_shell_num(context);

  if (prim_num < shell_num) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_ao_basis_shell_num",
                           "prim_num < shell_num");
  }

  int32_t mask = 1 << 2;
  ctx->ao_basis.prim_num = prim_num;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_shell_center(qmckl_context context, const int64_t* shell_center) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 3;

  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_center",
                           "shell_num is not set");
  }

  if (ctx->ao_basis.shell_center != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_center);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_shell_center",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_center",
                           NULL);
  }

  memcpy(new_array, shell_center, mem_info.size);

  ctx->ao_basis.shell_center = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const char* shell_ang_mom) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 4;

  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_ang_mom",
                           "shell_num is not set");
  }

  if (ctx->ao_basis.shell_ang_mom != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_ang_mom);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_shell_ang_mom",
                             NULL);
    }
  }


  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(char);
  char* new_array = (char*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_ang_mom",
                           NULL);
  }

  memcpy(new_array, shell_ang_mom, mem_info.size);

  ctx->ao_basis.shell_ang_mom = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num(qmckl_context context, const int64_t* shell_prim_num) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 5;

  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_prim_num",
                           "shell_num is not set");
  }

  if (ctx->ao_basis.shell_prim_num != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_num);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_shell_prim_num",
                             NULL);
    }
  }


  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_prim_num",
                           NULL);
  }

  memcpy(new_array, shell_prim_num, mem_info.size);

  ctx->ao_basis.shell_prim_num = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_shell_prim_index(qmckl_context context, const int64_t* shell_prim_index) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 6;

  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_prim_index",
                           "shell_num is not set");
  }

  if (ctx->ao_basis.shell_prim_index != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_index);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_shell_prim_index",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(int64_t);
  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_prim_index",
                           NULL);
  }

  memcpy(new_array, shell_prim_index, mem_info.size);

  ctx->ao_basis.shell_prim_index = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}


qmckl_exit_code  qmckl_set_ao_basis_shell_factor(qmckl_context context, const double* shell_factor) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 7;

  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
  if (shell_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_shell_factor",
                           "shell_num is not set");
  }


  if (ctx->ao_basis.shell_factor != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_factor);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_shell_factor",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = shell_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_shell_factor",
                           NULL);
  }

  memcpy(new_array, shell_factor, mem_info.size);

  ctx->ao_basis.shell_factor = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}

qmckl_exit_code  qmckl_set_ao_basis_exponent(qmckl_context context, const double* exponent) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 8;

  const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
  if (prim_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_exponent",
                           "prim_num is not set");
  }

  if (ctx->ao_basis.exponent != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.exponent);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_exponent",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prim_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_exponent",
                           NULL);
  }

  memcpy(new_array, exponent, mem_info.size);

  ctx->ao_basis.exponent = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}

qmckl_exit_code  qmckl_set_ao_basis_coefficient(qmckl_context context, const double* coefficient) {
  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
   }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;

  int32_t mask = 1 << 9;

  const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
  if (prim_num == 0L) {
    return qmckl_failwith( context,
                           QMCKL_FAILURE,
                           "qmckl_set_ao_basis_coefficient",
                           "prim_num is not set");
  }

  if (ctx->ao_basis.coefficient != NULL) {
    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.coefficient);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context, rc, 
                             "qmckl_set_ao_basis_coefficient",
                             NULL);
    }
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prim_num * sizeof(double);
  double* new_array = (double*) qmckl_malloc(context, mem_info);

  if (new_array == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_set_ao_basis_coefficient",
                           NULL);
  }

  memcpy(new_array, coefficient, mem_info.size);

  ctx->ao_basis.coefficient = new_array;

  ctx->ao_basis.uninitialized &= ~mask; 

  if (ctx->ao_basis.uninitialized == 0) {
    ctx->ao_basis.provided = 1;
  }

  return QMCKL_SUCCESS;
}

1.4 Fortran interfaces

qmcklcontext context in Global state
int64t n in Number of values
double X[n] in Array containing the input values
int32t LMAX[n] in Array containing the maximum power for each value
double P[n][ldp] out Array containing all the powers of X
int64t ldp in Leading dimension of array P

1.5 Test

/* Reference input data */

char    typ = 'G';
#define shell_num   ((int64_t) 12)
#define prim_num    ((int64_t) 20)

int64_t shell_center [shell_num] =
  { 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2 };

char shell_ang_mom [shell_num] =
  { 'S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D' };

double shell_factor [shell_num] =
  { 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1. };

int64_t shell_prim_num [shell_num] =
  {5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1};

int64_t shell_prim_index [shell_num] =
  {1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20};

double exponent [prim_num] =
  { 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
    1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
    0.3258, 0.1027, 1.407, 0.388, 1.057 };

double coefficient [prim_num] =
  { 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
    1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
    0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0 };
/* --- */

qmckl_exit_code rc;

munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_type (context, typ);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_num (context, shell_num);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_prim_num (context, prim_num);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_center (context, shell_center);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_center  (context, shell_prim_num);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_exponent      (context, exponent);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);

rc = qmckl_set_ao_basis_coefficient   (context, coefficient);
munit_assert_int64(rc, ==, QMCKL_SUCCESS);
munit_assert_int(qmckl_ao_basis_provided(context), ==, 1);

2 Polynomial part

2.1 Powers of \(x-X_i\)

The qmckl_ao_power function computes all the powers of the n input data up to the given maximum value given in input for each of the \(n\) points:

\[ P_{ik} = X_i^k \]

qmcklcontext context in Global state
int64t n in Number of values
double X[n] in Array containing the input values
int32t LMAX[n] in Array containing the maximum power for each value
double P[n][ldp] out Array containing all the powers of X
int64t ldp in Leading dimension of array P

2.1.1 Requirements:

  • context is not QMCKL_NULL_CONTEXT
  • n > 0
  • X is allocated with at least \(n \times 8\) bytes
  • LMAX is allocated with at least \(n \times 4\) bytes
  • P is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytes
  • LDP >= \(\max_i\) LMAX[i]

2.1.2 C Header

qmckl_exit_code qmckl_ao_power (
              const qmckl_context context,
      const int64_t n,
      const double* X,
      const int32_t* LMAX,
      double* const P,
      const int64_t ldp );

2.1.3 Source

integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
  use qmckl
  implicit none
  integer*8 , intent(in)  :: context
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: X(n)
  integer   , intent(in)  :: LMAX(n)
  real*8    , intent(out) :: P(ldp,n)
  integer*8 , intent(in)  :: ldp

  integer*8  :: i,k

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= ldp) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  k = MAXVAL(LMAX)
  if (LDP < k) then
     info = QMCKL_INVALID_ARG_6
     return
  endif

  if (k <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  do i=1,n
     P(1,i) = X(i)
     do k=2,LMAX(i)
        P(k,i) = P(k-1,i) * X(i) 
     end do
  end do

end function qmckl_ao_power_f

2.1.4 C interface

2.1.5 Fortran interface

2.1.6 Test

integer(c_int32_t) function test_qmckl_ao_power(context) bind(C)
  use qmckl
  implicit none

  integer(qmckl_context), intent(in), value :: context

  integer*8                     :: n, LDP 
  integer, allocatable          :: LMAX(:) 
  double precision, allocatable :: X(:), P(:,:)
  integer*8                     :: i,j
  double precision              :: epsilon

  epsilon = qmckl_get_numprec_epsilon(context)

  n = 100;
  LDP = 10;

  allocate(X(n), P(LDP,n), LMAX(n))

  do j=1,n
     X(j) = -5.d0 + 0.1d0 * dble(j)
     LMAX(j) = 1 + int(mod(j, 5),4)
  end do

  test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) 
  if (test_qmckl_ao_power /= QMCKL_SUCCESS) return

  test_qmckl_ao_power = QMCKL_FAILURE

  do j=1,n
     do i=1,LMAX(j)
        if ( X(j)**i == 0.d0 ) then
           if ( P(i,j) /= 0.d0) return
        else
           if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return
        end if
     end do
  end do

  test_qmckl_ao_power = QMCKL_SUCCESS
  deallocate(X,P,LMAX)
end function test_qmckl_ao_power

2.2 Value, Gradient and Laplacian of a polynomial

A polynomial is centered on a nucleus \(\mathbf{R}_i\)

\[ P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c \]

The gradients with respect to electron coordinates are

\begin{eqnarray*} \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ \end{eqnarray*}

and the Laplacian is

\begin{eqnarray*} \left( \frac{\partial }{\partial x^2} + \frac{\partial }{\partial y^2} + \frac{\partial }{\partial z^2} \right) P_l \left(\mathbf{r},\mathbf{R}_i \right) & = & a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}. \end{eqnarray*}

qmckl_ao_polynomial_vgl computes the values, gradients and Laplacians at a given point in space, of all polynomials with an angular momentum up to lmax.

qmcklcontext context in Global state
double X[3] in Array containing the coordinates of the points
double R[3] in Array containing the x,y,z coordinates of the center
int32t lmax in Maximum angular momentum
int64t n inout Number of computed polynomials
int32t L[n][ldl] out Contains a,b,c for all n results
int64t ldl in Leading dimension of L
double VGL[n][ldv] out Value, gradients and Laplacian of the polynomials
int64t ldv in Leading dimension of array VGL

2.2.1 Requirements:

  • context is not QMCKL_NULL_CONTEXT
  • n > 0
  • lmax >= 0
  • ldl >= 3
  • ldv >= 5
  • X is allocated with at least \(3 \times 8\) bytes
  • R is allocated with at least \(3 \times 8\) bytes
  • n >= (lmax+1)(lmax+2)(lmax+3)/6
  • L is allocated with at least \(3 \times n \times 4\) bytes
  • VGL is allocated with at least \(5 \times n \times 8\) bytes
  • On output, n should be equal to (lmax+1)(lmax+2)(lmax+3)/6
  • On output, the powers are given in the following order (l=a+b+c):
    • Increasing values of l
    • Within a given value of l, alphabetical order of the string made by a*"x" + b*"y" + c*"z" (in Python notation). For example, with a=0, b=2 and c=1 the string is "yyz"

2.2.2 C Header

qmckl_exit_code qmckl_ao_polynomial_vgl (
              const qmckl_context context,
      const double* X,
      const double* R,
      const int32_t lmax,
      int64_t* n,
      int32_t* const L,
      const int64_t ldl,
      double* const VGL,
      const int64_t ldv );

2.2.3 Source

integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
  use qmckl
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer   , intent(in)  :: lmax
  integer*8 , intent(out) :: n
  integer   , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
  integer*8 , intent(in)  :: ldl
  real*8    , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
  integer*8 , intent(in)  :: ldv

  integer*8         :: i,j
  integer           :: a,b,c,d
  real*8            :: Y(3)
  integer           :: lmax_array(3)
  real*8            :: pows(-2:lmax,3)
  integer, external :: qmckl_ao_power_f
  double precision  :: xy, yz, xz
  double precision  :: da, db, dc, dd

  info = 0

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (lmax < 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldl < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (ldv < 5) then
     info = QMCKL_INVALID_ARG_9
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do

  lmax_array(1:3) = lmax
  if (lmax == 0) then
     VGL(1,1) = 1.d0
     vgL(2:5,1) = 0.d0
     l(1:3,1) = 0
     n=1
  else if (lmax > 0) then
     pows(-2:0,1:3) = 1.d0
     do i=1,lmax
        pows(i,1) = pows(i-1,1) * Y(1) 
        pows(i,2) = pows(i-1,2) * Y(2) 
        pows(i,3) = pows(i-1,3) * Y(3) 
     end do

     VGL(1:5,1:4) = 0.d0
     l  (1:3,1:4) = 0

     VGL(1  ,1  ) = 1.d0
     vgl(1:5,2:4) = 0.d0

     l  (1,2) = 1
     vgl(1,2) = pows(1,1)
     vgL(2,2) = 1.d0

     l  (2,3) = 1
     vgl(1,3) = pows(1,2)
     vgL(3,3) = 1.d0

     l  (3,4) = 1
     vgl(1,4) = pows(1,3)
     vgL(4,4) = 1.d0

     n=4
  endif

  ! l>=2
  dd = 2.d0
  do d=2,lmax
     da = dd
     do a=d,0,-1
        db = dd-da
        do b=d-a,0,-1
           c  = d  - a  - b
           dc = dd - da - db
           n = n+1

           l(1,n) = a
           l(2,n) = b
           l(3,n) = c

           xy = pows(a,1) * pows(b,2)
           yz = pows(b,2) * pows(c,3)
           xz = pows(a,1) * pows(c,3)

           vgl(1,n) = xy * pows(c,3)

           xy = dc * xy
           xz = db * xz
           yz = da * yz

           vgl(2,n) = pows(a-1,1) * yz
           vgl(3,n) = pows(b-1,2) * xz
           vgl(4,n) = pows(c-1,3) * xy

           vgl(5,n) = &
                (da-1.d0) * pows(a-2,1) * yz + &
                (db-1.d0) * pows(b-2,2) * xz + &
                (dc-1.d0) * pows(c-2,3) * xy

           db = db - 1.d0
        end do
        da = da - 1.d0
     end do
     dd = dd + 1.d0
  end do

  info = QMCKL_SUCCESS

end function qmckl_ao_polynomial_vgl_f

2.2.4 C interface

2.2.5 Fortran interface

2.2.6 Test

integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
  use qmckl
  implicit none

  integer(c_int64_t), intent(in), value :: context

  integer                       :: lmax, d, i
  integer, allocatable          :: L(:,:)
  integer*8                     :: n, ldl, ldv, j
  double precision              :: X(3), R(3), Y(3)
  double precision, allocatable :: VGL(:,:)
  double precision              :: w
  double precision              :: epsilon

  epsilon = qmckl_get_numprec_epsilon(context)

  X = (/ 1.1 , 2.2 ,  3.3 /)
  R = (/ 0.1 , 1.2 , -2.3 /)
  Y(:) = X(:) - R(:)

  lmax = 4;
  ldl = 3;
  ldv = 100;

  d = (lmax+1)*(lmax+2)*(lmax+3)/6

  allocate (L(ldl,d), VGL(ldv,d))

  test_qmckl_ao_polynomial_vgl = &
       qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)

  if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return
  if (n /= d) return

  do j=1,n
     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     do i=1,3
        if (L(i,j) < 0) return
     end do
     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     if (dabs(1.d0 - VGL(1,j) / (&
          Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j)  &
          )) > epsilon ) return

     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     if (L(1,j) < 1) then
        if (VGL(2,j) /= 0.d0) return
     else
        if (dabs(1.d0 - VGL(2,j) / (&
             L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) &
             )) > epsilon ) return
     end if

     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     if (L(2,j) < 1) then
        if (VGL(3,j) /= 0.d0) return
     else
        if (dabs(1.d0 - VGL(3,j) / (&
             L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) &
             )) > epsilon ) return
     end if

     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     if (L(3,j) < 1) then
        if (VGL(4,j) /= 0.d0) return
     else
        if (dabs(1.d0 - VGL(4,j) / (&
             L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
             )) > epsilon ) return
     end if

     test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE
     w = 0.d0
     if (L(1,j) > 1) then
        w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) 
     end if
     if (L(2,j) > 1) then
        w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j) 
     end if
     if (L(3,j) > 1) then
        w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) 
     end if
     if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return
  end do

  test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS

  deallocate(L,VGL)
end function test_qmckl_ao_polynomial_vgl

int  test_qmckl_ao_polynomial_vgl(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));

3 Radial part

3.1 Gaussian basis functions

qmckl_ao_gaussian_vgl computes the values, gradients and Laplacians at a given point of n Gaussian functions centered at the same point:

\[ v_i = \exp(-a_i |X-R|^2) \] \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \]

context input Global state
X(3) input Array containing the coordinates of the points
R(3) input Array containing the x,y,z coordinates of the center
n input Number of computed Gaussians
A(n) input Exponents of the Gaussians
VGL(ldv,5) output Value, gradients and Laplacian of the Gaussians
ldv input Leading dimension of array VGL

Requirements :

  • context is not 0
  • n > 0
  • ldv >= 5
  • A(i) > 0 for all i
  • X is allocated with at least \(3 \times 8\) bytes
  • R is allocated with at least \(3 \times 8\) bytes
  • A is allocated with at least \(n \times 8\) bytes
  • VGL is allocated with at least \(n \times 5 \times 8\) bytes
qmckl_exit_code
qmckl_ao_gaussian_vgl(const qmckl_context context,
                      const double *X,
                      const double *R,
                      const int64_t *n,
                      const int64_t *A,
                      const double *VGL,
                      const int64_t ldv);

integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info)
  use qmckl
  implicit none
  integer*8 , intent(in)  :: context
  real*8    , intent(in)  :: X(3), R(3)
  integer*8 , intent(in)  :: n
  real*8    , intent(in)  :: A(n)
  real*8    , intent(out) :: VGL(ldv,5)
  integer*8 , intent(in)  :: ldv

  integer*8         :: i,j
  real*8            :: Y(3), r2, t, u, v

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (n <= 0) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (ldv < n) then
     info = QMCKL_INVALID_ARG_7
     return
  endif


  do i=1,3
     Y(i) = X(i) - R(i)
  end do
  r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3)

  do i=1,n
     VGL(i,1) = dexp(-A(i) * r2)
  end do

  do i=1,n
     VGL(i,5) = A(i) * VGL(i,1)
  end do

  t = -2.d0 * ( X(1) - R(1) )
  u = -2.d0 * ( X(2) - R(2) )
  v = -2.d0 * ( X(3) - R(3) )

  do i=1,n
     VGL(i,2) = t * VGL(i,5)
     VGL(i,3) = u * VGL(i,5)
     VGL(i,4) = v * VGL(i,5)
  end do

  t = 4.d0 * r2
  do i=1,n
     VGL(i,5) = (t * A(i) - 6.d0) *  VGL(i,5)
  end do

end function qmckl_ao_gaussian_vgl_f

integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C)
  use qmckl
  implicit none

  integer(c_int64_t), intent(in), value :: context

  integer*8                     :: n, ldv, j, i
  double precision              :: X(3), R(3), Y(3), r2
  double precision, allocatable :: VGL(:,:), A(:)
  double precision              :: epsilon

  epsilon = qmckl_get_numprec_epsilon(context)

  X = (/ 1.1 , 2.2 ,  3.3 /)
  R = (/ 0.1 , 1.2 , -2.3 /)
  Y(:) = X(:) - R(:)
  r2 = Y(1)**2 + Y(2)**2 + Y(3)**2

  n = 10;
  ldv = 100;

  allocate (A(n), VGL(ldv,5))
  do i=1,n
     A(i) = 0.0013 * dble(ishft(1,i))
  end do


  test_qmckl_ao_gaussian_vgl = &
       qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv)
  if (test_qmckl_ao_gaussian_vgl /= 0) return

  test_qmckl_ao_gaussian_vgl = -1

  do i=1,n
     test_qmckl_ao_gaussian_vgl = -11
     if (dabs(1.d0 - VGL(i,1) / (&
          dexp(-A(i) * r2) &
          )) > epsilon ) return

     test_qmckl_ao_gaussian_vgl = -12
     if (dabs(1.d0 - VGL(i,2) / (&
          -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) &
          )) > epsilon ) return

     test_qmckl_ao_gaussian_vgl = -13
     if (dabs(1.d0 - VGL(i,3) / (&
          -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) &
          )) > epsilon ) return

     test_qmckl_ao_gaussian_vgl = -14
     if (dabs(1.d0 - VGL(i,4) / (&
          -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) &
          )) > epsilon ) return

     test_qmckl_ao_gaussian_vgl = -15
     if (dabs(1.d0 - VGL(i,5) / (&
          A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) &
          )) > epsilon ) return
  end do

  test_qmckl_ao_gaussian_vgl = 0

  deallocate(VGL)
end function test_qmckl_ao_gaussian_vgl

3.2 TODO Slater basis functions

3.3 TODO Radial functions on a grid

4 Combining radial and polynomial parts

Author: TREX CoE

Created: 2021-04-20 Tue 22:01

Validate