TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity

Added eN cusp fitting in MOs

Browse Source
This commit is contained in:
Anthony Scemama 2023-03-17 14:54:58 +01:00
parent e10c7584ff
commit 2ecfc55dbc
1 changed files with 253 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

360
org/qmckl_mo.org
View File

@ -71,6 +71,8 @@ cusp is exact.
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_MO_HPF
#define QMCKL_MO_HPF
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval h_private_type)
@ -78,6 +80,7 @@ cusp is exact.
#define QMCKL_MO_HPT
#include <stdbool.h>
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval c_test) :noweb yes
@ -140,11 +143,11 @@ int main() {
Computed data:
|-----------------+--------------------------+----------------------------------------------------------------------------------|
| ~r_cusp_param~ | ~[nucl_num][mo_num][4]~ | 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 |
|-----------------+--------------------------+----------------------------------------------------------------------------------|
|--------------+--------------------------+-----------------------------------------------------------|
| ~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
@ -157,7 +160,7 @@ typedef struct qmckl_mo_basis_struct {
double * restrict mo_vgl;
double * restrict mo_value;
qmckl_tensor r_cusp_param;
qmckl_tensor cusp_param;
uint64_t mo_vgl_date;
uint64_t mo_value_date;
@ -431,36 +434,6 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
double* coord = (double*) malloc(ctx->nucleus.num * 3 * sizeof(double));
// Evaluate MO vgl at r_cusp
qmckl_tensor mo_vgl_at_r_cusp;
{
int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
mo_vgl_at_r_cusp = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
int64_t i=0;
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
coord[i] += r_cusp[inucl];
i+=3;
}
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
QMCKL_FAILURE,
"qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at r_cusp");
}
rc = qmckl_get_mo_basis_mo_vgl(context,
&(qmckl_ten3(mo_vgl_at_r_cusp,0,0,0)),
ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
}
// Set r_cusp
{
assert (ctx->mo_basis.r_cusp == NULL);
@ -480,14 +453,14 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
// Allocate cusp parameters and set them to zero
{
if (ctx->mo_basis.r_cusp_param.size[0] != 0) {
rc = qmckl_tensor_free(context, &(ctx->mo_basis.r_cusp_param));
if (ctx->mo_basis.cusp_param.size[0] != 0) {
rc = qmckl_tensor_free(context, &(ctx->mo_basis.cusp_param));
if (rc != QMCKL_SUCCESS) return rc;
}
int64_t sze[3] = { 4, ctx->mo_basis.mo_num, ctx->nucleus.num };
ctx->mo_basis.r_cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
ctx->mo_basis.r_cusp_param = qmckl_tensor_set(ctx->mo_basis.r_cusp_param, 0.);
int64_t sze[3] = { ctx->mo_basis.mo_num, 4, ctx->nucleus.num };
ctx->mo_basis.cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
ctx->mo_basis.cusp_param = qmckl_tensor_set(ctx->mo_basis.cusp_param, 0.);
}
@ -515,20 +488,24 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
// Evaluate MO vgl at r_cusp without s components
qmckl_tensor mo_vgl_at_r_cusp_no_s;
qmckl_tensor mo_vgl_at_r_cusp_s;
{
int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num };
mo_vgl_at_r_cusp_no_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
int64_t sze[3] = { ctx->mo_basis.mo_num, 3, ctx->nucleus.num };
mo_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
}
{
qmckl_tensor ao_vgl_at_r_cusp_s;
int64_t sze[3] = { ctx->ao_basis.ao_num, 5, ctx->nucleus.num };
ao_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc;
int64_t i=0;
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
coord[i] += r_cusp[inucl];
i+=3;
for (int64_t i=0 ; i<ctx->nucleus.num ; ++i) {
coord[2*ctx->nucleus.num + i] += r_cusp[i];
}
rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3));
if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context,
@ -537,11 +514,28 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
"Unable to set coordinates at r_cusp");
}
rc = qmckl_get_mo_basis_mo_vgl(context,
&(qmckl_ten3(mo_vgl_at_r_cusp_no_s,0,0,0)),
ctx->mo_basis.mo_num * 5 * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc;
rc = qmckl_get_ao_basis_ao_vgl(context,
&(qmckl_ten3(ao_vgl_at_r_cusp_s,0,0,0)),
ctx->ao_basis.ao_num * 5 * ctx->point.num);
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = 0.;
qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) = 0.;
qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) = 0.;
for (int64_t k=0 ; k<ctx->ao_basis.ao_num ; ++k) {
if ( ctx->ao_basis.ao_nucl[k] == inucl && ctx->ao_basis.ao_ang_mom[k] == 0) {
const double ck = ctx->mo_basis.coefficient[k + i*ctx->ao_basis.ao_num];
qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,0,inucl);
qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,3,inucl);
qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,4,inucl);
}
}
}
}
rc = qmckl_tensor_free(context,&ao_vgl_at_r_cusp_s);
if (rc != QMCKL_SUCCESS) return rc;
}
// Compute parameters
@ -551,36 +545,30 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
if (Z < 0.1) continue; // Avoid dummy atoms
const double R = r_cusp[inucl];
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
const double phi = qmckl_ten3(mo_vgl_at_r_cusp,i,0,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,0,inucl);
const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,1,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,1,inucl);
const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,4,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,4,inucl);
const double phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl);
const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl);
const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl);
const double eta = qmckl_mat(mo_value_at_nucl_no_s,i,inucl);
qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl) =
qmckl_ten3(ctx->mo_basis.cusp_param,i,0,inucl) =
-(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl) =
qmckl_ten3(ctx->mo_basis.cusp_param,i,1,inucl) =
(lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.);
qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl) =
qmckl_ten3(ctx->mo_basis.cusp_param,i,2,inucl) =
-(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R);
qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl) =
qmckl_ten3(ctx->mo_basis.cusp_param,i,3,inucl) =
(R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R);
printf("%ld %ld %f %f %f %f\n",i, inucl,
qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl),
qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl));
}
}
}
free(coord);
qmckl_matrix_free(context, &mo_value_at_nucl_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp_s);
// Restore old points
if (old_point_num > 0) {
@ -1115,6 +1103,7 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance,
ctx->mo_basis.r_cusp,
ctx->mo_basis.cusp_param,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_value,
ctx->mo_basis.mo_value);
@ -1195,7 +1184,7 @@ end function qmckl_compute_mo_basis_mo_value_doc_f
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value (
const qmckl_context context,
const int64_t ao_num,
@ -1209,7 +1198,7 @@ end function qmckl_compute_mo_basis_mo_value_doc_f
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
const qmckl_context context,
const int64_t ao_num,
@ -1267,7 +1256,7 @@ qmckl_compute_mo_basis_mo_value (const qmckl_context context,
*** HPC version
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
@ -1582,7 +1571,10 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
ctx->ao_basis.ao_nucl,
ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance,
ctx->nucleus.coord,
ctx->point.coord,
ctx->mo_basis.r_cusp,
ctx->mo_basis.cusp_param,
ctx->mo_basis.coefficient_t,
ctx->ao_basis.ao_vgl,
ctx->mo_basis.mo_vgl);
@ -1678,7 +1670,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl (
const qmckl_context context,
const int64_t ao_num,
@ -1692,7 +1684,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc (
const qmckl_context context,
const int64_t ao_num,
@ -1750,7 +1742,7 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
*** HPC version
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
@ -1910,6 +1902,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
| ~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 |
@ -1920,7 +1913,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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, coefficient_t, ao_value, mo_value) &
r_cusp, cusp_param, coefficient_t, ao_value, mo_value) &
result(info)
use qmckl
implicit none
@ -1930,31 +1923,44 @@ integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, &
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) :: 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)
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
integer*8 :: i, j, k, inucl
double precision :: r
info = QMCKL_SUCCESS
do j=1,point_num
mo_value(:,j) = 0.d0
do i=1,point_num
mo_value(:,i) = 0.d0
do k=1,ao_num
if (ao_value(k,j) == 0.d0) cycle
if (ao_value(k,i) == 0.d0) cycle
inucl = ao_nucl(k)+1
if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
end do
end do
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
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp"))
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp (
const qmckl_context context,
const int64_t nucl_num,
@ -1965,6 +1971,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
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 );
@ -1984,6 +1991,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
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 );
@ -2003,6 +2011,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
ao_ang_mom, &
en_distance, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_value, &
mo_value) &
@ -2020,6 +2029,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num)
real (c_double ) , intent(in) :: en_distance(nucl_num,point_num)
real (c_double ) , intent(in) :: r_cusp(nucl_num)
real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num)
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_value(ao_num,point_num)
real (c_double ) , intent(out) :: mo_value(mo_num,point_num)
@ -2035,6 +2045,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
ao_ang_mom, &
en_distance, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_value, &
mo_value)
@ -2053,26 +2064,34 @@ qmckl_compute_mo_basis_mo_value_cusp (const qmckl_context context,
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
return 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,
coefficient_t, ao_value, mo_value );
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
return 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,
coefficient_t, ao_value, mo_value );
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;
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
@ -2084,6 +2103,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
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 );
@ -2102,6 +2122,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
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)
@ -2158,14 +2179,30 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
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
#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;
}
@ -2192,7 +2229,10 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
| ~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 |
@ -2202,7 +2242,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
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, &
r_cusp, coefficient_t, ao_vgl, mo_vgl) &
nucl_coord, point_coord, r_cusp, cusp_param, coefficient_t, ao_vgl, mo_vgl) &
result(info)
use qmckl
implicit none
@ -2211,14 +2251,20 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
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
@ -2237,16 +2283,45 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
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
#+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp"))
# #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp (
const qmckl_context context,
const int64_t nucl_num,
@ -2256,7 +2331,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
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 );
@ -2275,7 +2353,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
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 );
@ -2294,7 +2375,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
ao_nucl, &
ao_ang_mom, &
en_distance, &
nucl_coord, &
point_coord, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_vgl, &
mo_vgl) &
@ -2311,7 +2395,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
integer (c_int64_t) , intent(in) :: ao_nucl(ao_num)
integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num)
real (c_double ) , intent(in) :: en_distance(nucl_num,point_num)
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
real (c_double ) , intent(in) :: point_coord(point_num,3)
real (c_double ) , intent(in) :: r_cusp(nucl_num)
real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num)
real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num)
real (c_double ) , intent(out) :: mo_vgl(mo_num,5,point_num)
@ -2326,7 +2413,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
ao_nucl, &
ao_ang_mom, &
en_distance, &
nucl_coord, &
point_coord, &
r_cusp, &
cusp_param, &
coefficient_t, &
ao_vgl, &
mo_vgl)
@ -2344,27 +2434,43 @@ qmckl_compute_mo_basis_mo_vgl_cusp (const qmckl_context context,
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
return qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp,
coefficient_t, ao_vgl, mo_vgl );
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
return qmckl_compute_mo_basis_mo_vgl_cusp_doc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp,
coefficient_t, ao_vgl, mo_vgl );
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;
}
#+end_src
*** HPC version
#+begin_src c :tangle (eval h_func) :comments org
#+begin_src c :tangle (eval h_private_func) :comments org
#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
@ -2375,7 +2481,10 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
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 );
@ -2393,7 +2502,10 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
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 )
@ -2512,6 +2624,40 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
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;
}