mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-06-30 00:44:52 +02:00
Merge branch 'master' of github.com:TREX-CoE/qmckl
This commit is contained in:
commit
fd38bb09b9
18
.github/workflows/test-build.yml
vendored
18
.github/workflows/test-build.yml
vendored
|
@ -15,7 +15,15 @@ jobs:
|
|||
steps:
|
||||
- uses: actions/checkout@v2
|
||||
- name: install dependencies
|
||||
run: sudo apt-get install emacs autoconf
|
||||
run: sudo apt-get install emacs autoconf libhdf5-dev
|
||||
- name: install trexio
|
||||
run: |
|
||||
wget https://github.com/TREX-CoE/trexio/releases/download/v1.0/trexio-1.0.0.tar.gz
|
||||
tar -zxf trexio-1.0.0.tar.gz
|
||||
cd trexio-1.0.0
|
||||
./configure
|
||||
make -j 8
|
||||
sudo make install
|
||||
- name: Run test
|
||||
run: |
|
||||
./autogen.sh
|
||||
|
@ -41,6 +49,14 @@ jobs:
|
|||
# brew install gfortran-10
|
||||
# brew install emacs
|
||||
# brew install autoconf automake libtool
|
||||
# - name: install trexio
|
||||
# run: |
|
||||
# wget https://github.com/TREX-CoE/trexio/releases/download/v1.0/trexio-1.0.0.tar.gz
|
||||
# tar -zxf trexio-1.0.0.tar.gz
|
||||
# cd trexio-1.0.0
|
||||
# ./configure
|
||||
# make -j 8
|
||||
# sudo make install
|
||||
# - name: Run test
|
||||
# run: |
|
||||
# ./autogen.sh
|
||||
|
|
465
org/qmckl_ao.org
465
org/qmckl_ao.org
|
@ -37,6 +37,9 @@ In this section we describe first how the basis set is stored in the
|
|||
context, and then we present the kernels used to compute the values,
|
||||
gradients and Laplacian of the atomic basis functions.
|
||||
|
||||
* TODO [0/1] Missing features :noexport:
|
||||
- [ ] Error messages to tell what is missing when initializing
|
||||
|
||||
* Headers :noexport:
|
||||
#+begin_src elisp :noexport :results none
|
||||
(org-babel-lob-ingest "../tools/lib.org")
|
||||
|
@ -115,19 +118,19 @@ int main() {
|
|||
|
||||
Computed data:
|
||||
|
||||
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------|
|
||||
|--------------------------+----------------------------+-----------------------------------------------------------------------------------------------|
|
||||
| ~coefficient_normalized~ | ~[prim_num]~ | Normalized primitive coefficients |
|
||||
| ~nucleus_prim_index~ | ~[nucl_num]~ | Index of the first primitive for each nucleus |
|
||||
| ~nucleus_max_ang_mom~ | ~[nucl_num]~ | Maximum angular momentum for each nucleus |
|
||||
| ~nucleus_range~ | ~[nucl_num]~ | Distance beyond which all the AOs are zero |
|
||||
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------|
|
||||
| ~primitive_vgl~ | ~[5][walk_num][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
|--------------------------+----------------------------+-----------------------------------------------------------------------------------------------|
|
||||
| ~primitive_vgl~ | ~[5][elec_num][prim_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~primitive_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~shell_vgl~ | ~[5][walk_num][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~shell_vgl~ | ~[5][elec_num][shell_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~shell_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions |
|
||||
| ~ao_vgl~ | ~[5][walk_num][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~ao_vgl~ | ~[5][elec_num][ao_num]~ | Value, gradients, Laplacian of the primitives at electron positions |
|
||||
| ~ao_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the AOs at electron positions |
|
||||
|--------------------------+--------------------------------------+-----------------------------------------------------------------------------------------------|
|
||||
|--------------------------+----------------------------+-----------------------------------------------------------------------------------------------|
|
||||
| ~nucl_shell_index~ | ~[nucl_num]~ | Index of the first shell for each nucleus |
|
||||
| ~exponent_sorted~ | ~[prim_num]~ | Array of exponents for sorted primitives |
|
||||
| ~coeff_norm_sorted~ | ~[prim_num]~ | Array of normalized coefficients for sorted primitives |
|
||||
|
@ -135,6 +138,7 @@ int main() {
|
|||
|
||||
For H_2 with the following basis set,
|
||||
|
||||
#+NAME: basis
|
||||
#+BEGIN_EXAMPLE
|
||||
HYDROGEN
|
||||
S 5
|
||||
|
@ -157,6 +161,7 @@ D 1
|
|||
|
||||
we have:
|
||||
|
||||
#+NAME: params
|
||||
#+BEGIN_EXAMPLE
|
||||
type = 'G'
|
||||
shell_num = 12
|
||||
|
@ -1834,7 +1839,7 @@ qmckl_exit_code qmckl_get_ao_basis_primitive_vgl(qmckl_context context, double*
|
|||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
||||
size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num;
|
||||
memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double));
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
|
@ -1873,7 +1878,7 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
|
|||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = ctx->ao_basis.prim_num * 5 * ctx->electron.num *
|
||||
ctx->electron.walk_num * sizeof(double);
|
||||
sizeof(double);
|
||||
double* primitive_vgl = (double*) qmckl_malloc(context, mem_info);
|
||||
|
||||
if (primitive_vgl == NULL) {
|
||||
|
@ -1891,7 +1896,6 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
|
|||
ctx->ao_basis.prim_num,
|
||||
ctx->electron.num,
|
||||
ctx->nucleus.num,
|
||||
ctx->electron.walk_num,
|
||||
ctx->ao_basis.nucleus_prim_index,
|
||||
ctx->electron.coord_new,
|
||||
ctx->nucleus.coord,
|
||||
|
@ -1926,16 +1930,15 @@ qmckl_exit_code qmckl_provide_ao_basis_primitive_vgl(qmckl_context context)
|
|||
| int64_t | prim_num | in | Number of primitives |
|
||||
| int64_t | elec_num | in | Number of electrons |
|
||||
| int64_t | nucl_num | in | Number of nuclei |
|
||||
| int64_t | walk_num | in | Number of walkers |
|
||||
| int64_t | nucleus_prim_index[nucl_num] | in | Index of the 1st primitive of each nucleus |
|
||||
| double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates |
|
||||
| double | elec_coord[3][elec_num] | in | Electron coordinates |
|
||||
| double | nucl_coord[3][elec_num] | in | Nuclear coordinates |
|
||||
| double | expo[prim_num] | in | Exponents of the primitives |
|
||||
| double | primitive_vgl[5][walk_num][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
|
||||
| double | primitive_vgl[5][elec_num][prim_num] | out | Value, gradients and Laplacian of the primitives |
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
|
||||
prim_num, elec_num, nucl_num, walk_num, &
|
||||
prim_num, elec_num, nucl_num, &
|
||||
nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) &
|
||||
result(info)
|
||||
use qmckl
|
||||
|
@ -1944,14 +1947,13 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
|
|||
integer*8 , intent(in) :: prim_num
|
||||
integer*8 , intent(in) :: nucl_num
|
||||
integer*8 , intent(in) :: elec_num
|
||||
integer*8 , intent(in) :: walk_num
|
||||
integer*8 , intent(in) :: nucleus_prim_index(nucl_num+1)
|
||||
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
double precision , intent(in) :: elec_coord(elec_num,3)
|
||||
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
||||
double precision , intent(in) :: expo(prim_num)
|
||||
double precision , intent(out) :: primitive_vgl(prim_num,elec_num,walk_num,5)
|
||||
double precision , intent(out) :: primitive_vgl(prim_num,elec_num,5)
|
||||
|
||||
integer*8 :: inucl, iprim, iwalk, ielec
|
||||
integer*8 :: inucl, iprim, ielec
|
||||
double precision :: x, y, z, two_a, ar2, r2, v, cutoff
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
@ -1962,11 +1964,10 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
|
|||
do inucl=1,nucl_num
|
||||
! C is zero-based, so shift bounds by one
|
||||
do iprim = nucleus_prim_index(inucl)+1, nucleus_prim_index(inucl+1)
|
||||
do iwalk = 1, walk_num
|
||||
do ielec = 1, elec_num
|
||||
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
|
||||
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
|
||||
z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3)
|
||||
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
|
||||
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
|
||||
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
|
||||
|
||||
r2 = x*x + y*y + z*z
|
||||
ar2 = expo(iprim)*r2
|
||||
|
@ -1975,16 +1976,15 @@ integer function qmckl_compute_ao_basis_primitive_gaussian_vgl_f(context, &
|
|||
v = dexp(-ar2)
|
||||
two_a = -2.d0 * expo(iprim) * v
|
||||
|
||||
primitive_vgl(iprim, ielec, iwalk, 1) = v
|
||||
primitive_vgl(iprim, ielec, iwalk, 2) = two_a * x
|
||||
primitive_vgl(iprim, ielec, iwalk, 3) = two_a * y
|
||||
primitive_vgl(iprim, ielec, iwalk, 4) = two_a * z
|
||||
primitive_vgl(iprim, ielec, iwalk, 5) = two_a * (3.d0 - 2.d0*ar2)
|
||||
primitive_vgl(iprim, ielec, 1) = v
|
||||
primitive_vgl(iprim, ielec, 2) = two_a * x
|
||||
primitive_vgl(iprim, ielec, 3) = two_a * y
|
||||
primitive_vgl(iprim, ielec, 4) = two_a * z
|
||||
primitive_vgl(iprim, ielec, 5) = two_a * (3.d0 - 2.d0*ar2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end function qmckl_compute_ao_basis_primitive_gaussian_vgl_f
|
||||
#+end_src
|
||||
|
@ -1995,7 +1995,6 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
|
|||
const int64_t prim_num,
|
||||
const int64_t elec_num,
|
||||
const int64_t nucl_num,
|
||||
const int64_t walk_num,
|
||||
const int64_t* nucleus_prim_index,
|
||||
const double* elec_coord,
|
||||
const double* nucl_coord,
|
||||
|
@ -2008,7 +2007,15 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
|
|||
#+RESULTS:
|
||||
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
||||
integer(c_int32_t) function qmckl_compute_ao_basis_primitive_gaussian_vgl &
|
||||
(context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl) &
|
||||
(context, &
|
||||
prim_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
nucleus_prim_index, &
|
||||
elec_coord, &
|
||||
nucl_coord, &
|
||||
expo, &
|
||||
primitive_vgl) &
|
||||
bind(C) result(info)
|
||||
|
||||
use, intrinsic :: iso_c_binding
|
||||
|
@ -2018,16 +2025,23 @@ qmckl_exit_code qmckl_compute_ao_basis_primitive_gaussian_vgl(
|
|||
integer (c_int64_t) , intent(in) , value :: prim_num
|
||||
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||
integer (c_int64_t) , intent(in) , value :: nucl_num
|
||||
integer (c_int64_t) , intent(in) , value :: walk_num
|
||||
integer (c_int64_t) , intent(in) :: nucleus_prim_index(nucl_num)
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3)
|
||||
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
||||
real (c_double ) , intent(in) :: expo(prim_num)
|
||||
real (c_double ) , intent(out) :: primitive_vgl(elec_num,walk_num,5,prim_num)
|
||||
real (c_double ) , intent(out) :: primitive_vgl(prim_num,elec_num,5)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_compute_ao_basis_primitive_gaussian_vgl_f
|
||||
info = qmckl_compute_ao_basis_primitive_gaussian_vgl_f &
|
||||
(context, prim_num, elec_num, nucl_num, walk_num, nucleus_prim_index, elec_coord, nucl_coord, expo, primitive_vgl)
|
||||
(context, &
|
||||
prim_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
nucleus_prim_index, &
|
||||
elec_coord, &
|
||||
nucl_coord, &
|
||||
expo, &
|
||||
primitive_vgl)
|
||||
|
||||
end function qmckl_compute_ao_basis_primitive_gaussian_vgl
|
||||
#+end_src
|
||||
|
@ -2056,40 +2070,25 @@ def lf(a,x,y):
|
|||
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
|
||||
|
||||
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
||||
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
||||
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
|
||||
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
|
||||
|
||||
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
||||
#double prim_vgl[prim_num][5][elec_num];
|
||||
a = 0.9059; x = elec_26_w1 ; y = nucl_1
|
||||
print ( "[7][0][0][26] : %e"% f(a,x,y))
|
||||
print ( "[7][1][0][26] : %e"% df(a,x,y,1))
|
||||
print ( "[7][2][0][26] : %e"% df(a,x,y,2))
|
||||
print ( "[7][3][0][26] : %e"% df(a,x,y,3))
|
||||
print ( "[7][4][0][26] : %e"% lf(a,x,y))
|
||||
|
||||
a = 0.32578; x = elec_15_w2 ; y = nucl_2
|
||||
print ( "[39][0][1][15] : %e"% f(a,x,y))
|
||||
print ( "[39][1][1][15] : %e"% df(a,x,y,1))
|
||||
print ( "[39][2][1][15] : %e"% df(a,x,y,2))
|
||||
print ( "[39][3][1][15] : %e"% df(a,x,y,3))
|
||||
print ( "[39][4][1][15] : %e"% lf(a,x,y))
|
||||
print ( "[7][0][26] : %e"% f(a,x,y))
|
||||
print ( "[7][1][26] : %e"% df(a,x,y,1))
|
||||
print ( "[7][2][26] : %e"% df(a,x,y,2))
|
||||
print ( "[7][3][26] : %e"% df(a,x,y,3))
|
||||
print ( "[7][4][26] : %e"% lf(a,x,y))
|
||||
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_example
|
||||
[7][0][0][26] : 1.050157e-03
|
||||
[7][1][0][26] : -7.501497e-04
|
||||
[7][2][0][26] : -3.825069e-03
|
||||
[7][3][0][26] : 3.495056e-03
|
||||
[7][4][0][26] : 2.040013e-02
|
||||
[39][0][1][15] : 1.083038e-03
|
||||
[39][1][1][15] : 2.378275e-03
|
||||
[39][2][1][15] : 2.143086e-03
|
||||
[39][3][1][15] : 4.332750e-04
|
||||
[39][4][1][15] : 7.514605e-03
|
||||
#+end_example
|
||||
: [7][0][26] : 1.050157e-03
|
||||
: [7][1][26] : -7.501497e-04
|
||||
: [7][2][26] : -3.825069e-03
|
||||
: [7][3][26] : 3.495056e-03
|
||||
: [7][4][26] : 2.040013e-02
|
||||
|
||||
*** Test
|
||||
|
||||
|
@ -2116,22 +2115,16 @@ assert(rc == QMCKL_SUCCESS);
|
|||
|
||||
|
||||
|
||||
double prim_vgl[5][walk_num][elec_num][prim_num];
|
||||
double prim_vgl[5][elec_num][prim_num];
|
||||
|
||||
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0][0]));
|
||||
rc = qmckl_get_ao_basis_primitive_vgl(context, &(prim_vgl[0][0][0]));
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
assert( fabs(prim_vgl[0][0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[1][0][26][7] - (-7.5014974095310560E-004)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[2][0][26][7] - (-3.8250692897610380E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[3][0][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[4][0][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
|
||||
|
||||
assert( fabs(prim_vgl[0][1][15][39] - ( 1.0825844173157661E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[1][1][15][39] - ( 2.3774237611651531E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[2][1][15][39] - ( 2.1423191526963063E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[3][1][15][39] - ( 4.3312003523048492E-004)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[0][26][7] - ( 1.0501570432064878E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[1][26][7] - (-7.5014974095310560E-004)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[2][26][7] - (-3.8250692897610380E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[3][26][7] - ( 3.4950559194080275E-003)) < 1.e-14 );
|
||||
assert( fabs(prim_vgl[4][26][7] - ( 2.0392163767356572E-002)) < 1.e-14 );
|
||||
|
||||
}
|
||||
|
||||
|
@ -2141,13 +2134,11 @@ assert( fabs(prim_vgl[4][1][15][39] - ( 7.5174404780004771E-003)) < 1.e-14 );
|
|||
*** Ideas for improvement
|
||||
|
||||
#+begin_src c
|
||||
// m : walkers
|
||||
// j : electrons
|
||||
// l : primitives
|
||||
|
||||
k=0;
|
||||
for (m=0 ; m<walk_num ; ++m) {
|
||||
for (j=0 ; j<elec_num ; ++j) {
|
||||
for (j=0 ; j<elec_num ; ++j) {
|
||||
for (i=0 ; i<nucl_num ; ++i) {
|
||||
|
||||
r2 = nucl_elec_dist[i][j];
|
||||
|
@ -2164,7 +2155,6 @@ for (m=0 ; m<walk_num ; ++m) {
|
|||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// sort(tmp) in increasing ar2;
|
||||
// Identify first ar2 above numerical accuracy threshold
|
||||
|
@ -2194,7 +2184,7 @@ qmckl_exit_code qmckl_get_ao_basis_shell_vgl(qmckl_context context, double* cons
|
|||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
||||
size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num;
|
||||
memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double));
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
|
@ -2253,8 +2243,7 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
|
|||
if (ctx->ao_basis.shell_vgl == NULL) {
|
||||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num *
|
||||
ctx->electron.walk_num * sizeof(double);
|
||||
mem_info.size = ctx->ao_basis.shell_num * 5 * ctx->electron.num * sizeof(double);
|
||||
double* shell_vgl = (double*) qmckl_malloc(context, mem_info);
|
||||
|
||||
if (shell_vgl == NULL) {
|
||||
|
@ -2273,7 +2262,6 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
|
|||
ctx->ao_basis.shell_num,
|
||||
ctx->electron.num,
|
||||
ctx->nucleus.num,
|
||||
ctx->electron.walk_num,
|
||||
ctx->ao_basis.nucleus_shell_num,
|
||||
ctx->ao_basis.nucleus_index,
|
||||
ctx->ao_basis.shell_prim_index,
|
||||
|
@ -2313,20 +2301,19 @@ qmckl_exit_code qmckl_provide_ao_basis_shell_vgl(qmckl_context context)
|
|||
| ~int64_t~ | ~shell_num~ | in | Number of shells |
|
||||
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
|
||||
| ~int64_t~ | ~nucl_num~ | in | Number of nuclei |
|
||||
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
|
||||
| ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells for each nucleus |
|
||||
| ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
||||
| ~int64_t~ | ~shell_prim_index[shell_num]~ | in | Index of the 1st primitive of each shell |
|
||||
| ~int64_t~ | ~shell_prim_num[shell_num]~ | in | Number of primitives per shell |
|
||||
| ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates |
|
||||
| ~double~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates |
|
||||
| ~double~ | ~nucl_coord[3][elec_num]~ | in | Nuclear coordinates |
|
||||
| ~double~ | ~expo[prim_num]~ | in | Exponents of the primitives |
|
||||
| ~double~ | ~coef_normalized[prim_num]~ | in | Coefficients of the primitives |
|
||||
| ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
|
||||
| ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | out | Value, gradients and Laplacian of the shells |
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
|
||||
prim_num, shell_num, elec_num, nucl_num, walk_num, &
|
||||
prim_num, shell_num, elec_num, nucl_num, &
|
||||
nucleus_shell_num, nucleus_index, shell_prim_index, shell_prim_num, &
|
||||
elec_coord, nucl_coord, expo, coef_normalized, shell_vgl) &
|
||||
result(info)
|
||||
|
@ -2337,18 +2324,17 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
|
|||
integer*8 , intent(in) :: shell_num
|
||||
integer*8 , intent(in) :: nucl_num
|
||||
integer*8 , intent(in) :: elec_num
|
||||
integer*8 , intent(in) :: walk_num
|
||||
integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
|
||||
integer*8 , intent(in) :: nucleus_index(nucl_num)
|
||||
integer*8 , intent(in) :: shell_prim_index(shell_num)
|
||||
integer*8 , intent(in) :: shell_prim_num(shell_num)
|
||||
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
double precision , intent(in) :: elec_coord(elec_num,3)
|
||||
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
||||
double precision , intent(in) :: expo(prim_num)
|
||||
double precision , intent(in) :: coef_normalized(prim_num)
|
||||
double precision , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
||||
double precision , intent(out) :: shell_vgl(shell_num,elec_num,5)
|
||||
|
||||
integer*8 :: inucl, iprim, iwalk, ielec, ishell
|
||||
integer*8 :: inucl, iprim, ielec, ishell
|
||||
double precision :: x, y, z, two_a, ar2, r2, v, cutoff
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
@ -2359,23 +2345,22 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
|
|||
|
||||
do inucl=1,nucl_num
|
||||
|
||||
do iwalk = 1, walk_num
|
||||
do ielec = 1, elec_num
|
||||
|
||||
x = elec_coord(ielec,1,iwalk) - nucl_coord(inucl,1)
|
||||
y = elec_coord(ielec,2,iwalk) - nucl_coord(inucl,2)
|
||||
z = elec_coord(ielec,3,iwalk) - nucl_coord(inucl,3)
|
||||
x = elec_coord(ielec,1) - nucl_coord(inucl,1)
|
||||
y = elec_coord(ielec,2) - nucl_coord(inucl,2)
|
||||
z = elec_coord(ielec,3) - nucl_coord(inucl,3)
|
||||
|
||||
r2 = x*x + y*y + z*z
|
||||
|
||||
do ishell=nucleus_index(inucl)+1, nucleus_index(inucl)+nucleus_shell_num(inucl)
|
||||
! C is zero-based, so shift bounds by one
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 1) = 0.d0
|
||||
shell_vgl(ishell, ielec, iwalk, 2) = 0.d0
|
||||
shell_vgl(ishell, ielec, iwalk, 3) = 0.d0
|
||||
shell_vgl(ishell, ielec, iwalk, 4) = 0.d0
|
||||
shell_vgl(ishell, ielec, iwalk, 5) = 0.d0
|
||||
shell_vgl(ishell, ielec, 1) = 0.d0
|
||||
shell_vgl(ishell, ielec, 2) = 0.d0
|
||||
shell_vgl(ishell, ielec, 3) = 0.d0
|
||||
shell_vgl(ishell, ielec, 4) = 0.d0
|
||||
shell_vgl(ishell, ielec, 5) = 0.d0
|
||||
|
||||
do iprim = shell_prim_index(ishell)+1, shell_prim_index(ishell)+shell_prim_num(ishell)
|
||||
|
||||
|
@ -2387,20 +2372,20 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
|
|||
v = coef_normalized(iprim) * dexp(-ar2)
|
||||
two_a = -2.d0 * expo(iprim) * v
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 1) = &
|
||||
shell_vgl(ishell, ielec, iwalk, 1) + v
|
||||
shell_vgl(ishell, ielec, 1) = &
|
||||
shell_vgl(ishell, ielec, 1) + v
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 2) = &
|
||||
shell_vgl(ishell, ielec, iwalk, 2) + two_a * x
|
||||
shell_vgl(ishell, ielec, 2) = &
|
||||
shell_vgl(ishell, ielec, 2) + two_a * x
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 3) = &
|
||||
shell_vgl(ishell, ielec, iwalk, 3) + two_a * y
|
||||
shell_vgl(ishell, ielec, 3) = &
|
||||
shell_vgl(ishell, ielec, 3) + two_a * y
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 4) = &
|
||||
shell_vgl(ishell, ielec, iwalk, 4) + two_a * z
|
||||
shell_vgl(ishell, ielec, 4) = &
|
||||
shell_vgl(ishell, ielec, 4) + two_a * z
|
||||
|
||||
shell_vgl(ishell, ielec, iwalk, 5) = &
|
||||
shell_vgl(ishell, ielec, iwalk, 5) + two_a * (3.d0 - 2.d0*ar2)
|
||||
shell_vgl(ishell, ielec, 5) = &
|
||||
shell_vgl(ishell, ielec, 5) + two_a * (3.d0 - 2.d0*ar2)
|
||||
|
||||
end do
|
||||
|
||||
|
@ -2408,7 +2393,6 @@ integer function qmckl_compute_ao_basis_shell_gaussian_vgl_f(context, &
|
|||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
||||
#+end_src
|
||||
|
@ -2423,7 +2407,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|||
const int64_t shell_num,
|
||||
const int64_t elec_num,
|
||||
const int64_t nucl_num,
|
||||
const int64_t walk_num,
|
||||
const int64_t* nucleus_shell_num,
|
||||
const int64_t* nucleus_index,
|
||||
const int64_t* shell_prim_index,
|
||||
|
@ -2445,7 +2428,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|||
shell_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
walk_num, &
|
||||
nucleus_shell_num, &
|
||||
nucleus_index, &
|
||||
shell_prim_index, &
|
||||
|
@ -2465,16 +2447,15 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|||
integer (c_int64_t) , intent(in) , value :: shell_num
|
||||
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||
integer (c_int64_t) , intent(in) , value :: nucl_num
|
||||
integer (c_int64_t) , intent(in) , value :: walk_num
|
||||
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
|
||||
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
|
||||
integer (c_int64_t) , intent(in) :: shell_prim_index(shell_num)
|
||||
integer (c_int64_t) , intent(in) :: shell_prim_num(shell_num)
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3)
|
||||
real (c_double ) , intent(in) :: nucl_coord(elec_num,3)
|
||||
real (c_double ) , intent(in) :: expo(prim_num)
|
||||
real (c_double ) , intent(in) :: coef_normalized(prim_num)
|
||||
real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
||||
real (c_double ) , intent(out) :: shell_vgl(shell_num,elec_num,5)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
||||
info = qmckl_compute_ao_basis_shell_gaussian_vgl_f &
|
||||
|
@ -2483,7 +2464,6 @@ end function qmckl_compute_ao_basis_shell_gaussian_vgl_f
|
|||
shell_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
walk_num, &
|
||||
nucleus_shell_num, &
|
||||
nucleus_index, &
|
||||
shell_prim_index, &
|
||||
|
@ -2525,7 +2505,7 @@ elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
|||
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
|
||||
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
|
||||
|
||||
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
||||
#double prim_vgl[prim_num][5][elec_num];
|
||||
x = elec_26_w1 ; y = nucl_1
|
||||
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
|
||||
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
|
||||
|
@ -2538,40 +2518,20 @@ a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
|
|||
( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ),
|
||||
( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )]
|
||||
|
||||
print ( "[1][0][0][26] : %25.15e"% f(a,x,y))
|
||||
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1))
|
||||
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2))
|
||||
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3))
|
||||
print ( "[1][4][0][26] : %25.15e"% lf(a,x,y))
|
||||
|
||||
x = elec_15_w2 ; y = nucl_2
|
||||
a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01),
|
||||
(5.095000E+00, 4.530800E-02 *2.4169531573445120e+00),
|
||||
(1.159000E+00, 2.028220E-01 *7.9610924849766440e-01),
|
||||
(3.258000E-01, 5.039030E-01 *3.0734305383061117e-01),
|
||||
(1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)]
|
||||
|
||||
print ( "[0][1][15][14] : %25.15e"% f(a,x,y))
|
||||
print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1))
|
||||
print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
|
||||
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
|
||||
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
|
||||
print ( "[1][0][26] : %25.15e"% f(a,x,y))
|
||||
print ( "[1][1][26] : %25.15e"% df(a,x,y,1))
|
||||
print ( "[1][2][26] : %25.15e"% df(a,x,y,2))
|
||||
print ( "[1][3][26] : %25.15e"% df(a,x,y,3))
|
||||
print ( "[1][4][26] : %25.15e"% lf(a,x,y))
|
||||
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_example
|
||||
[1][0][0][26] : 3.564393437193867e-02
|
||||
[1][1][0][26] : -6.030177988891605e-03
|
||||
[1][2][0][26] : -3.074832579871845e-02
|
||||
[1][3][0][26] : 2.809546963133958e-02
|
||||
[1][4][0][26] : 1.903338597841753e-02
|
||||
[0][1][15][14] : 5.928089771361000e-03
|
||||
[1][1][15][14] : 4.355862298893037e-03
|
||||
[2][1][15][14] : 3.925108924950765e-03
|
||||
[3][1][15][14] : 7.935527764416084e-04
|
||||
[4][1][15][14] : 2.697495005143935e-03
|
||||
#+end_example
|
||||
: [1][0][26] : 3.564393437193867e-02
|
||||
: [1][1][26] : -6.030177988891605e-03
|
||||
: [1][2][26] : -3.074832579871845e-02
|
||||
: [1][3][26] : 2.809546963133958e-02
|
||||
: [1][4][26] : 1.903338597841753e-02
|
||||
|
||||
*** Test
|
||||
|
||||
|
@ -2597,34 +2557,23 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
|
||||
double shell_vgl[5][walk_num][elec_num][shell_num];
|
||||
double shell_vgl[5][elec_num][shell_num];
|
||||
|
||||
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0][0]));
|
||||
rc = qmckl_get_ao_basis_shell_vgl(context, &(shell_vgl[0][0][0]));
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
printf(" shell_vgl[1][0][0][26] %25.15e\n", shell_vgl[0][0][26][1]);
|
||||
printf(" shell_vgl[1][1][0][26] %25.15e\n", shell_vgl[1][0][26][1]);
|
||||
printf(" shell_vgl[1][2][0][26] %25.15e\n", shell_vgl[2][0][26][1]);
|
||||
printf(" shell_vgl[1][3][0][26] %25.15e\n", shell_vgl[3][0][26][1]);
|
||||
printf(" shell_vgl[1][4][0][26] %25.15e\n", shell_vgl[4][0][26][1]);
|
||||
printf(" shell_vgl[1][0][26] %25.15e\n", shell_vgl[0][26][1]);
|
||||
printf(" shell_vgl[1][1][26] %25.15e\n", shell_vgl[1][26][1]);
|
||||
printf(" shell_vgl[1][2][26] %25.15e\n", shell_vgl[2][26][1]);
|
||||
printf(" shell_vgl[1][3][26] %25.15e\n", shell_vgl[3][26][1]);
|
||||
printf(" shell_vgl[1][4][26] %25.15e\n", shell_vgl[4][26][1]);
|
||||
|
||||
printf(" shell_vgl[14][0][1][15] %25.15e\n", shell_vgl[0][1][15][14]);
|
||||
printf(" shell_vgl[14][1][1][15] %25.15e\n", shell_vgl[1][1][15][14]);
|
||||
printf(" shell_vgl[14][2][1][15] %25.15e\n", shell_vgl[2][1][15][14]);
|
||||
printf(" shell_vgl[14][3][1][15] %25.15e\n", shell_vgl[3][1][15][14]);
|
||||
printf(" shell_vgl[14][4][1][15] %25.15e\n", shell_vgl[4][1][15][14]);
|
||||
assert( fabs(shell_vgl[0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[1][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[2][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[3][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[4][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
|
||||
|
||||
assert( fabs(shell_vgl[0][0][26][1] - ( 3.564393437193868e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[1][0][26][1] - (-6.030177987072189e-03)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[2][0][26][1] - (-3.074832579537582e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[3][0][26][1] - ( 2.809546963519935e-02)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[4][0][26][1] - ( 1.896046117183968e-02)) < 1.e-14 );
|
||||
|
||||
assert( fabs(shell_vgl[0][1][15][14] - ( 5.928089771361000e-03)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[1][1][15][14] - ( 4.355862296021654e-03)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[2][1][15][14] - ( 3.925108924923650e-03)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[3][1][15][14] - ( 7.935527784022099e-04)) < 1.e-14 );
|
||||
assert( fabs(shell_vgl[4][1][15][14] - ( 2.708246573703548e-03)) < 1.e-14 );
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
@ -3240,7 +3189,7 @@ qmckl_exit_code qmckl_get_ao_vgl(qmckl_context context, double* const ao_vgl) {
|
|||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num * ctx->electron.walk_num;
|
||||
size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num;
|
||||
memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double));
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
|
@ -3307,8 +3256,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|||
if (ctx->ao_basis.ao_vgl == NULL) {
|
||||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num *
|
||||
ctx->electron.walk_num * sizeof(double);
|
||||
mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->electron.num * sizeof(double);
|
||||
double* ao_vgl = (double*) qmckl_malloc(context, mem_info);
|
||||
|
||||
if (ao_vgl == NULL) {
|
||||
|
@ -3325,7 +3273,6 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|||
ctx->ao_basis.shell_num,
|
||||
ctx->electron.num,
|
||||
ctx->nucleus.num,
|
||||
ctx->electron.walk_num,
|
||||
ctx->electron.coord_new,
|
||||
ctx->nucleus.coord,
|
||||
ctx->ao_basis.nucleus_index,
|
||||
|
@ -3360,8 +3307,7 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|||
| ~int64_t~ | ~shell_num~ | in | Number of shells |
|
||||
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
|
||||
| ~int64_t~ | ~nucl_num~ | in | Number of nuclei |
|
||||
| ~int64_t~ | ~walk_num~ | in | Number of walkers |
|
||||
| ~double~ | ~elec_coord[walk_num][3][elec_num]~ | in | Electron coordinates |
|
||||
| ~double~ | ~elec_coord[3][elec_num]~ | in | Electron coordinates |
|
||||
| ~double~ | ~nucl_coord[3][nucl_num]~ | in | Nuclear coordinates |
|
||||
| ~int64_t~ | ~nucleus_index[nucl_num]~ | in | Index of the 1st shell of each nucleus |
|
||||
| ~int64_t~ | ~nucleus_shell_num[nucl_num]~ | in | Number of shells per nucleus |
|
||||
|
@ -3369,12 +3315,12 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context)
|
|||
| ~int32_t~ | ~nucleus_max_ang_mom[nucl_num]~ | in | Maximum angular momentum per nucleus |
|
||||
| ~int32_t~ | ~shell_ang_mom[shell_num]~ | in | Angular momentum of each shell |
|
||||
| ~double~ | ~ao_factor[ao_num]~ | in | Normalization factor of the AOs |
|
||||
| ~double~ | ~shell_vgl[5][walk_num][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells |
|
||||
| ~double~ | ~ao_vgl[5][walk_num][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
|
||||
| ~double~ | ~shell_vgl[5][elec_num][shell_num]~ | in | Value, gradients and Laplacian of the shells |
|
||||
| ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | out | Value, gradients and Laplacian of the AOs |
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
integer function qmckl_compute_ao_vgl_f(context, &
|
||||
ao_num, shell_num, elec_num, nucl_num, walk_num, &
|
||||
ao_num, shell_num, elec_num, nucl_num, &
|
||||
elec_coord, nucl_coord, nucleus_index, nucleus_shell_num, &
|
||||
nucleus_range, nucleus_max_ang_mom, shell_ang_mom, &
|
||||
ao_factor, shell_vgl, ao_vgl) &
|
||||
|
@ -3386,8 +3332,7 @@ integer function qmckl_compute_ao_vgl_f(context, &
|
|||
integer*8 , intent(in) :: shell_num
|
||||
integer*8 , intent(in) :: elec_num
|
||||
integer*8 , intent(in) :: nucl_num
|
||||
integer*8 , intent(in) :: walk_num
|
||||
double precision , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
double precision , intent(in) :: elec_coord(elec_num,3)
|
||||
double precision , intent(in) :: nucl_coord(nucl_num,3)
|
||||
integer*8 , intent(in) :: nucleus_index(nucl_num)
|
||||
integer*8 , intent(in) :: nucleus_shell_num(nucl_num)
|
||||
|
@ -3395,13 +3340,13 @@ integer function qmckl_compute_ao_vgl_f(context, &
|
|||
integer , intent(in) :: nucleus_max_ang_mom(nucl_num)
|
||||
integer , intent(in) :: shell_ang_mom(shell_num)
|
||||
double precision , intent(in) :: ao_factor(ao_num)
|
||||
double precision , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
||||
double precision , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
||||
double precision , intent(in) :: shell_vgl(shell_num,elec_num,5)
|
||||
double precision , intent(out) :: ao_vgl(ao_num,elec_num,5)
|
||||
|
||||
double precision :: e_coord(3), n_coord(3)
|
||||
integer*8 :: n_poly
|
||||
integer :: l, il, k
|
||||
integer*8 :: ielec, inucl, ishell, iwalk
|
||||
integer*8 :: ielec, inucl, ishell
|
||||
integer :: lstart(0:20)
|
||||
double precision :: x, y, z, r2
|
||||
double precision :: cutoff
|
||||
|
@ -3423,11 +3368,10 @@ integer function qmckl_compute_ao_vgl_f(context, &
|
|||
! TODO : Use numerical precision here
|
||||
cutoff = -dlog(1.d-15)
|
||||
|
||||
do iwalk = 1,walk_num
|
||||
do ielec = 1, elec_num
|
||||
e_coord(1) = elec_coord(ielec,1,iwalk)
|
||||
e_coord(2) = elec_coord(ielec,2,iwalk)
|
||||
e_coord(3) = elec_coord(ielec,3,iwalk)
|
||||
e_coord(1) = elec_coord(ielec,1)
|
||||
e_coord(2) = elec_coord(ielec,2)
|
||||
e_coord(3) = elec_coord(ielec,3)
|
||||
k=1
|
||||
do inucl=1,nucl_num
|
||||
n_coord(1) = nucl_coord(inucl,1)
|
||||
|
@ -3455,35 +3399,35 @@ integer function qmckl_compute_ao_vgl_f(context, &
|
|||
l = shell_ang_mom(ishell)
|
||||
do il = lstart(l), lstart(l+1)-1
|
||||
! Value
|
||||
ao_vgl(k,ielec,iwalk,1) = &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,1) * ao_factor(k)
|
||||
ao_vgl(k,ielec,1) = &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,1) * ao_factor(k)
|
||||
|
||||
! Grad_x
|
||||
ao_vgl(k,ielec,iwalk,2) = ( &
|
||||
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,2) &
|
||||
ao_vgl(k,ielec,2) = ( &
|
||||
poly_vgl(2,il) * shell_vgl(ishell,ielec,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,2) &
|
||||
) * ao_factor(k)
|
||||
|
||||
! Grad_y
|
||||
ao_vgl(k,ielec,iwalk,3) = ( &
|
||||
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,3) &
|
||||
ao_vgl(k,ielec,3) = ( &
|
||||
poly_vgl(3,il) * shell_vgl(ishell,ielec,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,3) &
|
||||
) * ao_factor(k)
|
||||
|
||||
! Grad_z
|
||||
ao_vgl(k,ielec,iwalk,4) = ( &
|
||||
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,4) &
|
||||
ao_vgl(k,ielec,4) = ( &
|
||||
poly_vgl(4,il) * shell_vgl(ishell,ielec,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,4) &
|
||||
) * ao_factor(k)
|
||||
|
||||
! Lapl_z
|
||||
ao_vgl(k,ielec,iwalk,5) = ( &
|
||||
poly_vgl(5,il) * shell_vgl(ishell,ielec,iwalk,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,iwalk,5) + &
|
||||
ao_vgl(k,ielec,5) = ( &
|
||||
poly_vgl(5,il) * shell_vgl(ishell,ielec,1) + &
|
||||
poly_vgl(1,il) * shell_vgl(ishell,ielec,5) + &
|
||||
2.d0 * ( &
|
||||
poly_vgl(2,il) * shell_vgl(ishell,ielec,iwalk,2) + &
|
||||
poly_vgl(3,il) * shell_vgl(ishell,ielec,iwalk,3) + &
|
||||
poly_vgl(4,il) * shell_vgl(ishell,ielec,iwalk,4) ) &
|
||||
poly_vgl(2,il) * shell_vgl(ishell,ielec,2) + &
|
||||
poly_vgl(3,il) * shell_vgl(ishell,ielec,3) + &
|
||||
poly_vgl(4,il) * shell_vgl(ishell,ielec,4) ) &
|
||||
) * ao_factor(k)
|
||||
|
||||
k = k+1
|
||||
|
@ -3491,7 +3435,6 @@ integer function qmckl_compute_ao_vgl_f(context, &
|
|||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(poly_vgl, powers)
|
||||
end function qmckl_compute_ao_vgl_f
|
||||
|
@ -3507,7 +3450,6 @@ end function qmckl_compute_ao_vgl_f
|
|||
const int64_t shell_num,
|
||||
const int64_t elec_num,
|
||||
const int64_t nucl_num,
|
||||
const int64_t walk_num,
|
||||
const double* elec_coord,
|
||||
const double* nucl_coord,
|
||||
const int64_t* nucleus_index,
|
||||
|
@ -3530,7 +3472,6 @@ end function qmckl_compute_ao_vgl_f
|
|||
shell_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
walk_num, &
|
||||
elec_coord, &
|
||||
nucl_coord, &
|
||||
nucleus_index, &
|
||||
|
@ -3551,8 +3492,7 @@ end function qmckl_compute_ao_vgl_f
|
|||
integer (c_int64_t) , intent(in) , value :: shell_num
|
||||
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||
integer (c_int64_t) , intent(in) , value :: nucl_num
|
||||
integer (c_int64_t) , intent(in) , value :: walk_num
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num)
|
||||
real (c_double ) , intent(in) :: elec_coord(elec_num,3)
|
||||
real (c_double ) , intent(in) :: nucl_coord(nucl_num,3)
|
||||
integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num)
|
||||
integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num)
|
||||
|
@ -3560,8 +3500,8 @@ end function qmckl_compute_ao_vgl_f
|
|||
integer (c_int32_t) , intent(in) :: nucleus_max_ang_mom(nucl_num)
|
||||
integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num)
|
||||
real (c_double ) , intent(in) :: ao_factor(ao_num)
|
||||
real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,walk_num,5)
|
||||
real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,walk_num,5)
|
||||
real (c_double ) , intent(in) :: shell_vgl(shell_num,elec_num,5)
|
||||
real (c_double ) , intent(out) :: ao_vgl(ao_num,elec_num,5)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_compute_ao_vgl_f
|
||||
info = qmckl_compute_ao_vgl_f &
|
||||
|
@ -3570,7 +3510,6 @@ end function qmckl_compute_ao_vgl_f
|
|||
shell_num, &
|
||||
elec_num, &
|
||||
nucl_num, &
|
||||
walk_num, &
|
||||
elec_coord, &
|
||||
nucl_coord, &
|
||||
nucleus_index, &
|
||||
|
@ -3587,6 +3526,7 @@ end function qmckl_compute_ao_vgl_f
|
|||
|
||||
#+begin_src python :results output :exports none
|
||||
import numpy as np
|
||||
from math import sqrt
|
||||
|
||||
def f(a,x,y):
|
||||
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
|
||||
|
@ -3612,7 +3552,7 @@ elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
|||
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
||||
nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] )
|
||||
|
||||
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
||||
#double ao_vgl[prim_num][5][elec_num];
|
||||
x = elec_26_w1 ; y = nucl_1
|
||||
a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
|
||||
( 121.170000, 0.012672 * 7.2836806319891484e+03),
|
||||
|
@ -3623,40 +3563,41 @@ a = [( 403.830000, 0.001473 * 5.9876577632594533e+04),
|
|||
( 1.763600, 0.273774 * 4.4423176930919421e+00),
|
||||
( 0.706190, 0.074397 * 8.9541051939952665e-01)]
|
||||
|
||||
print ( "[0][0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
|
||||
print ( "[1][0][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )
|
||||
norm = sqrt(3.)
|
||||
print ( "[0][26][219] : %25.15e"%(f(a,x,y) * (x[0] - y[0])**2) )
|
||||
print ( "[1][26][219] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + 2.*f(a,x,y) * (x[0] - y[0])) )
|
||||
|
||||
print ( "[0][0][26][220] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1])) )
|
||||
print ( "[1][0][26][220] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + f(a,x,y) * (x[1] - y[1])) )
|
||||
print ( "[0][26][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) ))
|
||||
print ( "[1][26][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) )
|
||||
|
||||
print ( "[0][0][26][221] : %25.15e"%(f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
|
||||
print ( "[1][0][26][221] : %25.15e"%(df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + f(a,x,y) * (x[2] - y[2])) )
|
||||
print ( "[0][26][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) )
|
||||
print ( "[1][26][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) )
|
||||
|
||||
print ( "[0][0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
|
||||
print ( "[1][0][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
|
||||
print ( "[0][26][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) )
|
||||
print ( "[1][26][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) )
|
||||
|
||||
print ( "[0][0][26][223] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
|
||||
print ( "[1][0][26][223] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
|
||||
print ( "[0][26][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) )
|
||||
print ( "[1][26][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) )
|
||||
|
||||
print ( "[0][0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
|
||||
print ( "[1][0][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
|
||||
print ( "[0][26][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) )
|
||||
print ( "[1][26][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) )
|
||||
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_example
|
||||
[0][0][26][219] : 1.020302912653649e-08
|
||||
[1][0][26][219] : -4.153046808203204e-08
|
||||
[0][0][26][220] : 8.756380857379661e-09
|
||||
[1][0][26][220] : -4.460176677299534e-08
|
||||
[0][0][26][221] : -2.705688401075445e-09
|
||||
[1][0][26][221] : 1.378177639720419e-08
|
||||
[0][0][26][222] : 7.514847283937212e-09
|
||||
[1][0][26][222] : -4.025905373647693e-08
|
||||
[0][0][26][223] : -2.322059246071533e-09
|
||||
[1][0][26][223] : 1.243989457599443e-08
|
||||
[0][0][26][224] : 7.175074806631758e-10
|
||||
[1][0][26][224] : -3.843880138733679e-09
|
||||
[0][26][219] : 1.020302912653649e-08
|
||||
[1][26][219] : -4.153046808203204e-08
|
||||
[0][26][220] : 1.516649653540510e-08
|
||||
[1][26][220] : -7.725252615816528e-08
|
||||
[0][26][221] : -4.686389780112468e-09
|
||||
[1][26][221] : 2.387073693851122e-08
|
||||
[0][26][222] : 7.514847283937212e-09
|
||||
[1][26][222] : -4.025905373647693e-08
|
||||
[0][26][223] : -4.021924592380977e-09
|
||||
[1][26][223] : 2.154652944642284e-08
|
||||
[0][26][224] : 7.175074806631758e-10
|
||||
[1][26][224] : -3.843880138733679e-09
|
||||
#+end_example
|
||||
|
||||
*** Test
|
||||
|
@ -3684,38 +3625,38 @@ rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
|||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
|
||||
double ao_vgl[5][walk_num][elec_num][ao_num];
|
||||
double ao_vgl[5][elec_num][ao_num];
|
||||
|
||||
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
|
||||
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0]));
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
printf("\n");
|
||||
printf(" ao_vgl ao_vgl[0][0][26][219] %25.15e\n", ao_vgl[0][0][26][219]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][219] %25.15e\n", ao_vgl[1][0][26][219]);
|
||||
printf(" ao_vgl ao_vgl[0][0][26][220] %25.15e\n", ao_vgl[0][0][26][220]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][220] %25.15e\n", ao_vgl[1][0][26][220]);
|
||||
printf(" ao_vgl ao_vgl[0][0][26][221] %25.15e\n", ao_vgl[0][0][26][221]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][221] %25.15e\n", ao_vgl[1][0][26][221]);
|
||||
printf(" ao_vgl ao_vgl[0][0][26][222] %25.15e\n", ao_vgl[0][0][26][222]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][222] %25.15e\n", ao_vgl[1][0][26][222]);
|
||||
printf(" ao_vgl ao_vgl[0][0][26][223] %25.15e\n", ao_vgl[0][0][26][223]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][223] %25.15e\n", ao_vgl[1][0][26][223]);
|
||||
printf(" ao_vgl ao_vgl[0][0][26][224] %25.15e\n", ao_vgl[0][0][26][224]);
|
||||
printf(" ao_vgl ao_vgl[1][0][26][224] %25.15e\n", ao_vgl[1][0][26][224]);
|
||||
printf(" ao_vgl ao_vgl[0][26][219] %25.15e\n", ao_vgl[0][26][219]);
|
||||
printf(" ao_vgl ao_vgl[1][26][219] %25.15e\n", ao_vgl[1][26][219]);
|
||||
printf(" ao_vgl ao_vgl[0][26][220] %25.15e\n", ao_vgl[0][26][220]);
|
||||
printf(" ao_vgl ao_vgl[1][26][220] %25.15e\n", ao_vgl[1][26][220]);
|
||||
printf(" ao_vgl ao_vgl[0][26][221] %25.15e\n", ao_vgl[0][26][221]);
|
||||
printf(" ao_vgl ao_vgl[1][26][221] %25.15e\n", ao_vgl[1][26][221]);
|
||||
printf(" ao_vgl ao_vgl[0][26][222] %25.15e\n", ao_vgl[0][26][222]);
|
||||
printf(" ao_vgl ao_vgl[1][26][222] %25.15e\n", ao_vgl[1][26][222]);
|
||||
printf(" ao_vgl ao_vgl[0][26][223] %25.15e\n", ao_vgl[0][26][223]);
|
||||
printf(" ao_vgl ao_vgl[1][26][223] %25.15e\n", ao_vgl[1][26][223]);
|
||||
printf(" ao_vgl ao_vgl[0][26][224] %25.15e\n", ao_vgl[0][26][224]);
|
||||
printf(" ao_vgl ao_vgl[1][26][224] %25.15e\n", ao_vgl[1][26][224]);
|
||||
printf("\n");
|
||||
|
||||
assert( fabs(ao_vgl[0][0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][219] - ( -4.928035238010602e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][0][26][220] - ( 8.756345547784206e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][220] - ( -4.460158690983819e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][0][26][221] - ( -2.705677490544664e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][221] - ( 1.378172082017231e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][222] - ( -4.025889138635182e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][0][26][223] - ( -2.322049882502961e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][223] - ( 1.243984441042288e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][0][26][224] - ( -3.843864637762753e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][219] - ( 1.020298798341620e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][219] - (-4.928035238010602e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][220] - ( 1.516643537739178e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][220] - (-7.725221462603871e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][221] - (-4.686370882518819e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][221] - ( 2.387064067626827e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][222] - ( 7.514816980753531e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][222] - (-4.025889138635182e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][223] - (-4.021908374204471e-09)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][223] - ( 2.154644255710413e-08)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[0][26][224] - ( 7.175045873560788e-10)) < 1.e-14 );
|
||||
assert( fabs(ao_vgl[1][26][224] - (-3.843864637762753e-09)) < 1.e-14 );
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
|
322
org/qmckl_blas.org
Normal file
322
org/qmckl_blas.org
Normal file
|
@ -0,0 +1,322 @@
|
|||
#+TITLE: BLAS functions
|
||||
#+SETUPFILE: ../tools/theme.setup
|
||||
#+INCLUDE: ../tools/lib.org
|
||||
|
||||
* Headers :noexport:
|
||||
#+begin_src elisp :noexport :results none
|
||||
(org-babel-lob-ingest "../tools/lib.org")
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
|
||||
#include "qmckl.h"
|
||||
#include "assert.h"
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif
|
||||
int main() {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
|
||||
#+end_src
|
||||
|
||||
* Matrix operations
|
||||
|
||||
** ~qmckl_dgemm~
|
||||
|
||||
Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.
|
||||
TODO: Add description about the external library dependence.
|
||||
|
||||
#+NAME: qmckl_dgemm_args
|
||||
| qmckl_context | context | in | Global state |
|
||||
| bool | TransA | in | Number of rows of the input matrix |
|
||||
| bool | TransB | in | Number of rows of the input matrix |
|
||||
| int64_t | m | in | Number of rows of the input matrix |
|
||||
| int64_t | n | in | Number of columns of the input matrix |
|
||||
| int64_t | k | in | Number of columns of the input matrix |
|
||||
| double | alpha | in | Number of columns of the input matrix |
|
||||
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
|
||||
| int64_t | lda | in | Leading dimension of array ~A~ |
|
||||
| double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ |
|
||||
| int64_t | ldb | in | Leading dimension of array ~B~ |
|
||||
| double | beta | in | Array containing the $n \times m$ matrix $B$ |
|
||||
| double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
|
||||
| int64_t | ldc | in | Leading dimension of array ~B~ |
|
||||
|
||||
*** Requirements
|
||||
|
||||
- ~context~ is not ~QMCKL_NULL_CONTEXT~
|
||||
- ~m > 0~
|
||||
- ~n > 0~
|
||||
- ~k > 0~
|
||||
- ~lda >= m~
|
||||
- ~ldb >= n~
|
||||
- ~ldc >= n~
|
||||
- ~A~ is allocated with at least $m \times k \times 8$ bytes
|
||||
- ~B~ is allocated with at least $k \times n \times 8$ bytes
|
||||
- ~C~ is allocated with at least $m \times n \times 8$ bytes
|
||||
|
||||
*** C header
|
||||
|
||||
#+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||
|
||||
#+RESULTS:
|
||||
#+BEGIN_src c :tangle (eval h_func) :comments org
|
||||
qmckl_exit_code qmckl_dgemm (
|
||||
const qmckl_context context,
|
||||
const bool TransA,
|
||||
const bool TransB,
|
||||
const int64_t m,
|
||||
const int64_t n,
|
||||
const int64_t k,
|
||||
const double alpha,
|
||||
const double* A,
|
||||
const int64_t lda,
|
||||
const double* B,
|
||||
const int64_t ldb,
|
||||
const double beta,
|
||||
double* const C,
|
||||
const int64_t ldc );
|
||||
#+END_src
|
||||
|
||||
|
||||
*** Source
|
||||
#+begin_src f90 :tangle (eval f)
|
||||
integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
|
||||
result(info)
|
||||
use qmckl
|
||||
implicit none
|
||||
integer(qmckl_context) , intent(in) :: context
|
||||
logical*8 , intent(in) :: TransA, TransB
|
||||
integer*8 , intent(in) :: m, n, k
|
||||
real*8 , intent(in) :: alpha, beta
|
||||
integer*8 , intent(in) :: lda
|
||||
real*8 , intent(in) :: A(m,k)
|
||||
integer*8 , intent(in) :: ldb
|
||||
real*8 , intent(in) :: B(k,n)
|
||||
integer*8 , intent(in) :: ldc
|
||||
real*8 , intent(out) :: C(m,n)
|
||||
real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:)
|
||||
|
||||
integer*8 :: i,j,l, LDA_2, LDB_2
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
|
||||
if (TransA) then
|
||||
allocate(AT(k,m))
|
||||
do i = 1, m
|
||||
do j = 1, k
|
||||
AT(j,i) = A(i,j)
|
||||
end do
|
||||
end do
|
||||
LDA_2 = M
|
||||
else
|
||||
LDA_2 = LDA
|
||||
endif
|
||||
|
||||
if (TransB) then
|
||||
allocate(BT(n,k))
|
||||
do i = 1, k
|
||||
do j = 1, n
|
||||
BT(j,i) = B(i,j)
|
||||
end do
|
||||
end do
|
||||
LDB_2 = K
|
||||
else
|
||||
LDB_2 = LDB
|
||||
endif
|
||||
|
||||
if (context == QMCKL_NULL_CONTEXT) then
|
||||
info = QMCKL_INVALID_CONTEXT
|
||||
return
|
||||
endif
|
||||
|
||||
if (m <= 0_8) then
|
||||
info = QMCKL_INVALID_ARG_4
|
||||
return
|
||||
endif
|
||||
|
||||
if (n <= 0_8) then
|
||||
info = QMCKL_INVALID_ARG_5
|
||||
return
|
||||
endif
|
||||
|
||||
if (k <= 0_8) then
|
||||
info = QMCKL_INVALID_ARG_6
|
||||
return
|
||||
endif
|
||||
|
||||
if (LDA_2 .ne. m) then
|
||||
info = QMCKL_INVALID_ARG_9
|
||||
return
|
||||
endif
|
||||
|
||||
if (LDB_2 .ne. k) then
|
||||
info = QMCKL_INVALID_ARG_10
|
||||
return
|
||||
endif
|
||||
|
||||
if (LDC .ne. m) then
|
||||
info = QMCKL_INVALID_ARG_13
|
||||
return
|
||||
endif
|
||||
|
||||
if (TransA) then
|
||||
C = matmul(AT,B)
|
||||
else if (TransB) then
|
||||
C = matmul(A,BT)
|
||||
else
|
||||
C = matmul(A,B)
|
||||
endif
|
||||
end function qmckl_dgemm_f
|
||||
#+end_src
|
||||
|
||||
*** C interface :noexport:
|
||||
|
||||
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||
|
||||
#+RESULTS:
|
||||
#+BEGIN_src f90 :tangle (eval f) :comments org :exports none
|
||||
integer(c_int32_t) function qmckl_dgemm &
|
||||
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
||||
bind(C) result(info)
|
||||
|
||||
use, intrinsic :: iso_c_binding
|
||||
implicit none
|
||||
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
logical*8 , intent(in) , value :: TransA
|
||||
logical*8 , intent(in) , value :: TransB
|
||||
integer (c_int64_t) , intent(in) , value :: m
|
||||
integer (c_int64_t) , intent(in) , value :: n
|
||||
integer (c_int64_t) , intent(in) , value :: k
|
||||
real (c_double ) , intent(in) , value :: alpha
|
||||
integer (c_int64_t) , intent(in) , value :: lda
|
||||
real (c_double ) , intent(in) :: A(lda,*)
|
||||
integer (c_int64_t) , intent(in) , value :: ldb
|
||||
real (c_double ) , intent(in) :: B(ldb,*)
|
||||
real (c_double ) , intent(in) , value :: beta
|
||||
integer (c_int64_t) , intent(in) , value :: ldc
|
||||
real (c_double ) , intent(out) :: C(ldc,*)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_dgemm_f
|
||||
info = qmckl_dgemm_f &
|
||||
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
|
||||
|
||||
end function qmckl_dgemm
|
||||
#+END_src
|
||||
|
||||
|
||||
#+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
||||
interface
|
||||
integer(c_int32_t) function qmckl_dgemm &
|
||||
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
|
||||
bind(C)
|
||||
use, intrinsic :: iso_c_binding
|
||||
import
|
||||
implicit none
|
||||
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
logical*8 , intent(in) , value :: TransA
|
||||
logical*8 , intent(in) , value :: TransB
|
||||
integer (c_int64_t) , intent(in) , value :: m
|
||||
integer (c_int64_t) , intent(in) , value :: n
|
||||
integer (c_int64_t) , intent(in) , value :: k
|
||||
real (c_double ) , intent(in) , value :: alpha
|
||||
integer (c_int64_t) , intent(in) , value :: lda
|
||||
real (c_double ) , intent(in) :: A(lda,*)
|
||||
integer (c_int64_t) , intent(in) , value :: ldb
|
||||
real (c_double ) , intent(in) :: B(ldb,*)
|
||||
real (c_double ) , intent(in) , value :: beta
|
||||
integer (c_int64_t) , intent(in) , value :: ldc
|
||||
real (c_double ) , intent(out) :: C(ldc,*)
|
||||
|
||||
end function qmckl_dgemm
|
||||
end interface
|
||||
#+end_src
|
||||
|
||||
*** Test :noexport:
|
||||
#+begin_src f90 :tangle (eval f_test)
|
||||
integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
|
||||
use qmckl
|
||||
implicit none
|
||||
integer(qmckl_context), intent(in), value :: context
|
||||
|
||||
double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:)
|
||||
integer*8 :: m, n, k, LDA, LDB, LDC
|
||||
integer*8 :: i,j,l
|
||||
logical*8 :: TransA, TransB
|
||||
double precision :: x, alpha, beta
|
||||
|
||||
TransA = .False.
|
||||
TransB = .False.
|
||||
m = 1_8
|
||||
k = 4_8
|
||||
n = 6_8
|
||||
LDA = m
|
||||
LDB = k
|
||||
LDC = m
|
||||
|
||||
allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
|
||||
|
||||
A = 0.d0
|
||||
B = 0.d0
|
||||
C = 0.d0
|
||||
D = 0.d0
|
||||
alpha = 1.0d0
|
||||
beta = 0.0d0
|
||||
do j=1,k
|
||||
do i=1,m
|
||||
A(i,j) = -10.d0 + dble(i+j)
|
||||
end do
|
||||
end do
|
||||
|
||||
do j=1,n
|
||||
do i=1,k
|
||||
B(i,j) = -10.d0 + dble(i+j)
|
||||
end do
|
||||
end do
|
||||
|
||||
test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC)
|
||||
|
||||
if (test_qmckl_dgemm /= QMCKL_SUCCESS) return
|
||||
|
||||
test_qmckl_dgemm = QMCKL_FAILURE
|
||||
|
||||
x = 0.d0
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
do l=1,k
|
||||
D(i,j) = D(i,j) + A(i,l)*B(l,j)
|
||||
end do
|
||||
x = x + (D(i,j) - C(i,j))**2
|
||||
end do
|
||||
end do
|
||||
|
||||
if (dabs(x) <= 1.d-15) then
|
||||
test_qmckl_dgemm = QMCKL_SUCCESS
|
||||
endif
|
||||
|
||||
deallocate(A,B,C,D)
|
||||
end function test_qmckl_dgemm
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments link :tangle (eval c_test)
|
||||
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
|
||||
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
|
||||
#+end_src
|
||||
|
||||
* End of files :noexport:
|
||||
|
||||
#+begin_src c :comments link :tangle (eval c_test)
|
||||
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
|
||||
return 0;
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
||||
|
||||
# -*- mode: org -*-
|
||||
# vim: syntax=c
|
|
@ -31,10 +31,12 @@ int main() {
|
|||
#include "qmckl_nucleus_private_type.h"
|
||||
#include "qmckl_electron_private_type.h"
|
||||
#include "qmckl_ao_private_type.h"
|
||||
#include "qmckl_mo_private_type.h"
|
||||
#include "qmckl_jastrow_private_type.h"
|
||||
#include "qmckl_nucleus_private_func.h"
|
||||
#include "qmckl_electron_private_func.h"
|
||||
#include "qmckl_ao_private_func.h"
|
||||
#include "qmckl_mo_private_func.h"
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c)
|
||||
|
@ -119,6 +121,7 @@ typedef struct qmckl_context_struct {
|
|||
qmckl_nucleus_struct nucleus;
|
||||
qmckl_electron_struct electron;
|
||||
qmckl_ao_basis_struct ao_basis;
|
||||
qmckl_mo_basis_struct mo_basis;
|
||||
qmckl_jastrow_struct jastrow;
|
||||
|
||||
/* To be implemented:
|
||||
|
|
|
@ -532,11 +532,11 @@ qmckl_set_electron_num(qmckl_context context,
|
|||
"up_num <= 0");
|
||||
}
|
||||
|
||||
if (down_num <= 0) {
|
||||
if (down_num < 0) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_INVALID_ARG_3,
|
||||
"qmckl_set_electron_num",
|
||||
"down_num <= 0");
|
||||
"down_num < 0");
|
||||
}
|
||||
|
||||
int32_t mask = 1 << 0;
|
||||
|
|
872
org/qmckl_mo.org
Normal file
872
org/qmckl_mo.org
Normal file
|
@ -0,0 +1,872 @@
|
|||
#+TITLE: Molecular Orbitals
|
||||
#+SETUPFILE: ../tools/theme.setup
|
||||
#+INCLUDE: ../tools/lib.org
|
||||
|
||||
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 SCF method)
|
||||
the MOs are defined as follows:
|
||||
|
||||
\[
|
||||
\phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r})
|
||||
\]
|
||||
|
||||
|
||||
In this section we demonstrate how to use the QMCkl specific DGEMM
|
||||
function to calculate the MOs.
|
||||
|
||||
|
||||
* Headers :noexport:
|
||||
#+begin_src elisp :noexport :results none
|
||||
(org-babel-lob-ingest "../tools/lib.org")
|
||||
#+end_src
|
||||
|
||||
|
||||
#+begin_src c :tangle (eval h_private_type)
|
||||
#ifndef QMCKL_MO_HPT
|
||||
#define QMCKL_MO_HPT
|
||||
|
||||
#include <stdbool.h>
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c_test) :noweb yes
|
||||
#include "qmckl.h"
|
||||
#include "assert.h"
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
#include "chbrclf.h"
|
||||
#include "qmckl_ao_private_func.h"
|
||||
#include "qmckl_mo_private_func.h"
|
||||
|
||||
int main() {
|
||||
qmckl_context context;
|
||||
context = qmckl_context_create();
|
||||
|
||||
qmckl_exit_code rc;
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c)
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif
|
||||
|
||||
#ifdef HAVE_STDINT_H
|
||||
#include <stdint.h>
|
||||
#elif HAVE_INTTYPES_H
|
||||
#include <inttypes.h>
|
||||
#endif
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <stdbool.h>
|
||||
#include <assert.h>
|
||||
|
||||
#include "qmckl.h"
|
||||
#include "qmckl_context_private_type.h"
|
||||
#include "qmckl_memory_private_type.h"
|
||||
#include "qmckl_memory_private_func.h"
|
||||
#include "qmckl_ao_private_type.h"
|
||||
#include "qmckl_ao_private_func.h"
|
||||
#include "qmckl_mo_private_type.h"
|
||||
#include "qmckl_mo_private_func.h"
|
||||
#+end_src
|
||||
|
||||
* Context
|
||||
|
||||
The following arrays are stored in the context:
|
||||
|
||||
|---------------------+--------------------+--------------------------------------------------------------|
|
||||
| ~type~ | | Gaussian (~'G'~) or Slater (~'S'~) |
|
||||
| ~mo_num~ | | Number of MOs |
|
||||
| ~coefficient~ | ~[mo_num, ao_num]~ | Orbital coefficients |
|
||||
|
||||
Computed data:
|
||||
|
||||
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||
| ~mo_vgl~ | ~[5][elec_num][mo_num]~ | Value, gradients, Laplacian of the MOs at electron positions |
|
||||
| ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at electron positions |
|
||||
|---------------+-----------------------------------+----------------------------------------------------------------------------------------|
|
||||
|
||||
** Data structure
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_type)
|
||||
typedef struct qmckl_mo_basis_struct {
|
||||
int64_t mo_num;
|
||||
double * coefficient;
|
||||
|
||||
double * mo_vgl;
|
||||
int64_t mo_vgl_date;
|
||||
|
||||
int32_t uninitialized;
|
||||
bool provided;
|
||||
char type;
|
||||
} qmckl_mo_basis_struct;
|
||||
#+end_src
|
||||
|
||||
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.
|
||||
|
||||
** Access functions
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
|
||||
char qmckl_get_mo_basis_type (const qmckl_context context);
|
||||
int64_t qmckl_get_mo_basis_mo_num (const qmckl_context context);
|
||||
double* qmckl_get_mo_basis_coefficient (const qmckl_context context);
|
||||
#+end_src
|
||||
|
||||
When all the data for the AOs have been provided, the following
|
||||
function returns ~true~.
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func)
|
||||
bool qmckl_mo_basis_provided (const qmckl_context context);
|
||||
#+end_src
|
||||
|
||||
#+NAME:post
|
||||
#+begin_src c :exports none
|
||||
if ( (ctx->mo_basis.uninitialized & mask) != 0) {
|
||||
return NULL;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
char qmckl_get_mo_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->mo_basis.uninitialized & mask) != 0) {
|
||||
return (char) 0;
|
||||
}
|
||||
|
||||
assert (ctx->mo_basis.type != (char) 0);
|
||||
return ctx->mo_basis.type;
|
||||
}
|
||||
|
||||
int64_t qmckl_get_mo_basis_mo_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->mo_basis.uninitialized & mask) != 0) {
|
||||
return (int64_t) 0;
|
||||
}
|
||||
|
||||
assert (ctx->mo_basis.mo_num > (int64_t) 0);
|
||||
return ctx->mo_basis.mo_num;
|
||||
}
|
||||
|
||||
bool qmckl_mo_basis_provided(const qmckl_context context) {
|
||||
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return false;
|
||||
}
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
return ctx->mo_basis.provided;
|
||||
}
|
||||
|
||||
|
||||
#+end_src
|
||||
|
||||
** Initialization functions
|
||||
|
||||
To set the basis set, all the following functions need to be
|
||||
called.
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func)
|
||||
qmckl_exit_code qmckl_set_mo_basis_type (qmckl_context context, const char t);
|
||||
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);
|
||||
#+end_src
|
||||
|
||||
#+NAME:pre2
|
||||
#+begin_src c :exports none
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_NULL_CONTEXT;
|
||||
}
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
#+end_src
|
||||
|
||||
#+NAME:post2
|
||||
#+begin_src c :exports none
|
||||
ctx->mo_basis.uninitialized &= ~mask;
|
||||
ctx->mo_basis.provided = (ctx->mo_basis.uninitialized == 0);
|
||||
if (ctx->mo_basis.provided) {
|
||||
qmckl_exit_code rc_ = qmckl_finalize_basis(context);
|
||||
if (rc_ != QMCKL_SUCCESS) return rc_;
|
||||
}
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_set_mo_basis_type(qmckl_context context, const char t) {
|
||||
<<pre2>>
|
||||
|
||||
if (t != 'G' && t != 'S') {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_INVALID_ARG_2,
|
||||
"qmckl_set_mo_basis_type",
|
||||
NULL);
|
||||
}
|
||||
|
||||
int32_t mask = 1;
|
||||
ctx->mo_basis.type = t;
|
||||
|
||||
<<post2>>
|
||||
}
|
||||
|
||||
qmckl_exit_code qmckl_set_mo_basis_mo_num(qmckl_context context, const int64_t mo_num) {
|
||||
<<pre2>>
|
||||
|
||||
if (mo_num <= 0) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_INVALID_ARG_2,
|
||||
"qmckl_set_mo_basis_mo_num",
|
||||
"mo_num <= 0");
|
||||
}
|
||||
|
||||
int32_t mask = 1 << 1;
|
||||
ctx->mo_basis.mo_num = mo_num;
|
||||
|
||||
<<post2>>
|
||||
}
|
||||
|
||||
qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) {
|
||||
<<pre2>>
|
||||
|
||||
int32_t mask = 1 << 2;
|
||||
|
||||
if (ctx->mo_basis.coefficient != NULL) {
|
||||
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
|
||||
if (rc != QMCKL_SUCCESS) {
|
||||
return qmckl_failwith( context, rc,
|
||||
"qmckl_set_mo_basis_coefficient",
|
||||
NULL);
|
||||
}
|
||||
}
|
||||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_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_mo_basis_coefficient",
|
||||
NULL);
|
||||
}
|
||||
|
||||
memcpy(new_array, coefficient, mem_info.size);
|
||||
|
||||
ctx->mo_basis.coefficient = new_array;
|
||||
|
||||
<<post2>>
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
||||
When the basis set is completely entered, other data structures are
|
||||
computed to accelerate the calculations.
|
||||
|
||||
* Computation
|
||||
|
||||
** Computation of MOs
|
||||
|
||||
*** Get
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func) :noweb yes
|
||||
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl);
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_get_mo_basis_vgl(qmckl_context context, double* const mo_vgl) {
|
||||
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_NULL_CONTEXT;
|
||||
}
|
||||
|
||||
qmckl_exit_code rc;
|
||||
|
||||
rc = qmckl_provide_ao_vgl(context);
|
||||
if (rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
rc = qmckl_provide_mo_vgl(context);
|
||||
if (rc != QMCKL_SUCCESS) return rc;
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
size_t sze = 5 * ctx->electron.num * ctx->mo_basis.mo_num;
|
||||
memcpy(mo_vgl, ctx->mo_basis.mo_vgl, sze * sizeof(double));
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
||||
interface
|
||||
integer(c_int32_t) function qmckl_get_mo_basis_vgl (context, mo_vgl) &
|
||||
bind(C)
|
||||
use, intrinsic :: iso_c_binding
|
||||
import
|
||||
implicit none
|
||||
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
double precision, intent(out) :: mo_vgl(*)
|
||||
end function
|
||||
end interface
|
||||
#+end_src
|
||||
|
||||
*** Provide
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context);
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
|
||||
qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
|
||||
{
|
||||
|
||||
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
||||
return QMCKL_NULL_CONTEXT;
|
||||
}
|
||||
|
||||
qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
|
||||
assert (ctx != NULL);
|
||||
|
||||
if (!ctx->ao_basis.provided) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_NOT_PROVIDED,
|
||||
"qmckl_ao_basis",
|
||||
NULL);
|
||||
}
|
||||
|
||||
if(!(ctx->electron.provided)) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_NOT_PROVIDED,
|
||||
"qmckl_electron",
|
||||
NULL);
|
||||
}
|
||||
|
||||
if (!ctx->mo_basis.provided) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_NOT_PROVIDED,
|
||||
"qmckl_mo_basis",
|
||||
NULL);
|
||||
}
|
||||
|
||||
/* Compute if necessary */
|
||||
if (ctx->electron.coord_new_date > ctx->mo_basis.mo_vgl_date) {
|
||||
|
||||
/* Allocate array */
|
||||
if (ctx->mo_basis.mo_vgl == NULL) {
|
||||
|
||||
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
|
||||
mem_info.size = 5 * ctx->electron.num * ctx->mo_basis.mo_num * sizeof(double);
|
||||
double* mo_vgl = (double*) qmckl_malloc(context, mem_info);
|
||||
|
||||
if (mo_vgl == NULL) {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_ALLOCATION_FAILED,
|
||||
"qmckl_mo_basis_vgl",
|
||||
NULL);
|
||||
}
|
||||
ctx->mo_basis.mo_vgl = mo_vgl;
|
||||
}
|
||||
|
||||
qmckl_exit_code rc;
|
||||
if (ctx->mo_basis.type == 'G') {
|
||||
rc = qmckl_compute_mo_basis_gaussian_vgl(context,
|
||||
ctx->ao_basis.ao_num,
|
||||
ctx->mo_basis.mo_num,
|
||||
ctx->electron.num,
|
||||
ctx->mo_basis.coefficient,
|
||||
ctx->ao_basis.ao_vgl,
|
||||
ctx->mo_basis.mo_vgl);
|
||||
} else {
|
||||
return qmckl_failwith( context,
|
||||
QMCKL_FAILURE,
|
||||
"compute_mo_basis_vgl",
|
||||
"Not yet implemented");
|
||||
}
|
||||
if (rc != QMCKL_SUCCESS) {
|
||||
return rc;
|
||||
}
|
||||
|
||||
ctx->mo_basis.mo_vgl_date = ctx->date;
|
||||
}
|
||||
|
||||
return QMCKL_SUCCESS;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Compute
|
||||
:PROPERTIES:
|
||||
:Name: qmckl_compute_mo_basis_gaussian_vgl
|
||||
:CRetType: qmckl_exit_code
|
||||
:FRetType: qmckl_exit_code
|
||||
:END:
|
||||
|
||||
#+NAME: qmckl_mo_basis_gaussian_vgl_args
|
||||
| ~qmckl_context~ | ~context~ | in | Global state |
|
||||
| ~int64_t~ | ~ao_num~ | in | Number of AOs |
|
||||
| ~int64_t~ | ~mo_num~ | in | Number of MOs |
|
||||
| ~int64_t~ | ~elec_num~ | in | Number of electrons |
|
||||
| ~double~ | ~coef_normalized[mo_num][ao_num]~ | in | AO to MO transformation matrix |
|
||||
| ~double~ | ~ao_vgl[5][elec_num][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
|
||||
| ~double~ | ~mo_vgl[5][elec_num][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
|
||||
|
||||
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
|
||||
integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, &
|
||||
ao_num, mo_num, elec_num, &
|
||||
coef_normalized, 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) :: elec_num
|
||||
double precision , intent(in) :: ao_vgl(ao_num,elec_num,5)
|
||||
double precision , intent(in) :: coef_normalized(ao_num,mo_num)
|
||||
double precision , intent(out) :: mo_vgl(mo_num,elec_num,5)
|
||||
logical*8 :: TransA, TransB
|
||||
double precision :: alpha, beta
|
||||
integer :: info_qmckl_dgemm_value
|
||||
integer :: info_qmckl_dgemm_Gx
|
||||
integer :: info_qmckl_dgemm_Gy
|
||||
integer :: info_qmckl_dgemm_Gz
|
||||
integer :: info_qmckl_dgemm_lap
|
||||
integer*8 :: M, N, K, LDA, LDB, LDC, i,j
|
||||
|
||||
integer*8 :: inucl, iprim, iwalk, ielec, ishell
|
||||
double precision :: x, y, z, two_a, ar2, r2, v, cutoff
|
||||
TransA = .False.
|
||||
TransB = .False.
|
||||
alpha = 1.0d0
|
||||
beta = 0.0d0
|
||||
|
||||
info = QMCKL_SUCCESS
|
||||
info_qmckl_dgemm_value = QMCKL_SUCCESS
|
||||
info_qmckl_dgemm_Gx = QMCKL_SUCCESS
|
||||
info_qmckl_dgemm_Gy = QMCKL_SUCCESS
|
||||
info_qmckl_dgemm_Gz = QMCKL_SUCCESS
|
||||
info_qmckl_dgemm_lap = QMCKL_SUCCESS
|
||||
|
||||
! Don't compute exponentials when the result will be almost zero.
|
||||
! TODO : Use numerical precision here
|
||||
cutoff = -dlog(1.d-15)
|
||||
M = 1_8
|
||||
N = mo_num * 1_8
|
||||
K = ao_num * 1_8
|
||||
LDA = M
|
||||
LDB = K
|
||||
LDC = M
|
||||
|
||||
do ielec = 1, elec_num
|
||||
! Value
|
||||
info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||
ao_vgl(:, ielec, 1), LDA, &
|
||||
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||
beta, &
|
||||
mo_vgl(:,ielec,1),LDC)
|
||||
! Grad_x
|
||||
info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||
ao_vgl(:, ielec, 2), LDA, &
|
||||
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||
beta, &
|
||||
mo_vgl(:,ielec,2),LDC)
|
||||
! Grad_y
|
||||
info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||
ao_vgl(:, ielec, 3), LDA, &
|
||||
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||
beta, &
|
||||
mo_vgl(:,ielec,3),LDC)
|
||||
! Grad_z
|
||||
info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, &
|
||||
ao_vgl(:, ielec, 4), LDA, &
|
||||
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||
beta, &
|
||||
mo_vgl(:,ielec,4),LDC)
|
||||
! Lapl_z
|
||||
info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, &
|
||||
ao_vgl(:, ielec, 5), LDA, &
|
||||
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
|
||||
beta, &
|
||||
mo_vgl(:,ielec,5),LDC)
|
||||
end do
|
||||
|
||||
if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. &
|
||||
info_qmckl_dgemm_Gx .eq. QMCKL_SUCCESS .and. &
|
||||
info_qmckl_dgemm_Gy .eq. QMCKL_SUCCESS .and. &
|
||||
info_qmckl_dgemm_Gz .eq. QMCKL_SUCCESS .and. &
|
||||
info_qmckl_dgemm_lap .eq. QMCKL_SUCCESS ) then
|
||||
info = QMCKL_SUCCESS
|
||||
else
|
||||
info = QMCKL_FAILURE
|
||||
end if
|
||||
|
||||
end function qmckl_compute_mo_basis_gaussian_vgl_f
|
||||
#+end_src
|
||||
|
||||
#+CALL: generate_c_header(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl"))
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_src c :tangle (eval h_func) :comments org
|
||||
qmckl_exit_code qmckl_compute_mo_basis_gaussian_vgl (
|
||||
const qmckl_context context,
|
||||
const int64_t ao_num,
|
||||
const int64_t mo_num,
|
||||
const int64_t elec_num,
|
||||
const double* coef_normalized,
|
||||
const double* ao_vgl,
|
||||
double* const mo_vgl );
|
||||
#+end_src
|
||||
|
||||
|
||||
#+CALL: generate_c_interface(table=qmckl_mo_basis_gaussian_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_gaussian_vgl"))
|
||||
|
||||
#+RESULTS:
|
||||
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
||||
integer(c_int32_t) function qmckl_compute_mo_basis_gaussian_vgl &
|
||||
(context, ao_num, mo_num, elec_num, coef_normalized, ao_vgl, mo_vgl) &
|
||||
bind(C) result(info)
|
||||
|
||||
use, intrinsic :: iso_c_binding
|
||||
implicit none
|
||||
|
||||
integer (c_int64_t) , intent(in) , value :: context
|
||||
integer (c_int64_t) , intent(in) , value :: ao_num
|
||||
integer (c_int64_t) , intent(in) , value :: mo_num
|
||||
integer (c_int64_t) , intent(in) , value :: elec_num
|
||||
real (c_double ) , intent(in) :: coef_normalized(ao_num,mo_num)
|
||||
real (c_double ) , intent(in) :: ao_vgl(ao_num,elec_num,5)
|
||||
real (c_double ) , intent(out) :: mo_vgl(mo_num,elec_num,5)
|
||||
|
||||
integer(c_int32_t), external :: qmckl_compute_mo_basis_gaussian_vgl_f
|
||||
info = qmckl_compute_mo_basis_gaussian_vgl_f &
|
||||
(context, ao_num, mo_num, elec_num, coef_normalized, ao_vgl, mo_vgl)
|
||||
|
||||
end function qmckl_compute_mo_basis_gaussian_vgl
|
||||
#+end_src
|
||||
|
||||
|
||||
|
||||
*** Test
|
||||
|
||||
#+begin_src python :results output :exports none
|
||||
import numpy as np
|
||||
|
||||
def f(a,x,y):
|
||||
return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] )
|
||||
|
||||
def df(a,x,y,n):
|
||||
h0 = 1.e-6
|
||||
if n == 1: h = np.array([h0,0.,0.])
|
||||
elif n == 2: h = np.array([0.,h0,0.])
|
||||
elif n == 3: h = np.array([0.,0.,h0])
|
||||
return ( f(a,x+h,y) - f(a,x-h,y) ) / (2.*h0)
|
||||
|
||||
def d2f(a,x,y,n):
|
||||
h0 = 1.e-6
|
||||
if n == 1: h = np.array([h0,0.,0.])
|
||||
elif n == 2: h = np.array([0.,h0,0.])
|
||||
elif n == 3: h = np.array([0.,0.,h0])
|
||||
return ( f(a,x+h,y) - 2.*f(a,x,y) + f(a,x-h,y) ) / h0**2
|
||||
|
||||
def lf(a,x,y):
|
||||
return d2f(a,x,y,1) + d2f(a,x,y,2) + d2f(a,x,y,3)
|
||||
|
||||
elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] )
|
||||
elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] )
|
||||
nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] )
|
||||
nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] )
|
||||
|
||||
#double prim_vgl[prim_num][5][walk_num][elec_num];
|
||||
x = elec_26_w1 ; y = nucl_1
|
||||
a = [( 8.236000E+03, -1.130000E-04 * 6.1616545431994848e+02 ),
|
||||
( 1.235000E+03, -8.780000E-04 * 1.4847738511079908e+02 ),
|
||||
( 2.808000E+02, -4.540000E-03 * 4.8888635917437597e+01 ),
|
||||
( 7.927000E+01, -1.813300E-02 * 1.8933972232608955e+01 ),
|
||||
( 2.559000E+01, -5.576000E-02 * 8.1089160941724145e+00 ),
|
||||
( 8.997000E+00, -1.268950E-01 * 3.7024003863155635e+00 ),
|
||||
( 3.319000E+00, -1.703520E-01 * 1.7525302846177560e+00 ),
|
||||
( 9.059000E-01, 1.403820E-01 * 6.6179013183966806e-01 ),
|
||||
( 3.643000E-01, 5.986840E-01 * 3.3419848027174592e-01 ),
|
||||
( 1.285000E-01, 3.953890E-01 * 1.5296336817449557e-01 )]
|
||||
|
||||
print ( "[1][0][0][26] : %25.15e"% f(a,x,y))
|
||||
print ( "[1][1][0][26] : %25.15e"% df(a,x,y,1))
|
||||
print ( "[1][2][0][26] : %25.15e"% df(a,x,y,2))
|
||||
print ( "[1][3][0][26] : %25.15e"% df(a,x,y,3))
|
||||
print ( "[1][4][0][26] : %25.15e"% lf(a,x,y))
|
||||
|
||||
x = elec_15_w2 ; y = nucl_2
|
||||
a = [(3.387000E+01, 6.068000E-03 *1.0006253235944540e+01),
|
||||
(5.095000E+00, 4.530800E-02 *2.4169531573445120e+00),
|
||||
(1.159000E+00, 2.028220E-01 *7.9610924849766440e-01),
|
||||
(3.258000E-01, 5.039030E-01 *3.0734305383061117e-01),
|
||||
(1.027000E-01, 3.834210E-01 *1.2929684417481876e-01)]
|
||||
|
||||
print ( "[0][1][15][14] : %25.15e"% f(a,x,y))
|
||||
print ( "[1][1][15][14] : %25.15e"% df(a,x,y,1))
|
||||
print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
|
||||
print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
|
||||
print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
|
||||
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c_test) :exports none
|
||||
{
|
||||
#define walk_num chbrclf_walk_num
|
||||
#define elec_num chbrclf_elec_num
|
||||
#define shell_num chbrclf_shell_num
|
||||
#define ao_num chbrclf_ao_num
|
||||
|
||||
int64_t elec_up_num = chbrclf_elec_up_num;
|
||||
int64_t elec_dn_num = chbrclf_elec_dn_num;
|
||||
double* elec_coord = &(chbrclf_elec_coord[0][0][0]);
|
||||
const int64_t nucl_num = chbrclf_nucl_num;
|
||||
const double* nucl_charge = chbrclf_charge;
|
||||
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
|
||||
|
||||
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_electron_walk_num (context, walk_num);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
assert(qmckl_electron_provided(context));
|
||||
|
||||
rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_nucleus_num (context, nucl_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]));
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_nucleus_charge(context, nucl_charge);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
assert(qmckl_nucleus_provided(context));
|
||||
|
||||
const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]);
|
||||
const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]);
|
||||
const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]);
|
||||
const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]);
|
||||
const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]);
|
||||
const double * shell_factor = &(chbrclf_basis_shell_factor[0]);
|
||||
const double * exponent = &(chbrclf_basis_exponent[0]);
|
||||
const double * coefficient = &(chbrclf_basis_coefficient[0]);
|
||||
const double * prim_factor = &(chbrclf_basis_prim_factor[0]);
|
||||
const double * ao_factor = &(chbrclf_basis_ao_factor[0]);
|
||||
|
||||
const char typ = 'G';
|
||||
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_type (context, typ);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_shell_factor (context, shell_factor);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_exponent (context, exponent);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_coefficient (context, coefficient);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
assert(!qmckl_ao_basis_provided(context));
|
||||
|
||||
rc = qmckl_set_ao_basis_prim_factor (context, prim_factor);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
rc = qmckl_set_ao_basis_ao_factor (context, ao_factor);
|
||||
assert(rc == QMCKL_SUCCESS);
|
||||
|
||||
assert(qmckl_ao_basis_provided(context));
|
||||
|
||||
|
||||
double ao_vgl[5][walk_num][elec_num][chbrclf_ao_num];
|
||||
|
||||
rc = qmckl_get_ao_vgl(context, &(ao_vgl[0][0][0][0]));
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
/* Set up MO data */
|
||||
rc = qmckl_set_mo_basis_type(context, typ);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
const int64_t mo_num = chbrclf_mo_num;
|
||||
rc = qmckl_set_mo_basis_mo_num(context, mo_num);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
const double * mo_coefficient = &(chbrclf_mo_coef[0]);
|
||||
|
||||
rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
assert(qmckl_mo_basis_provided(context));
|
||||
|
||||
double mo_vgl[5][elec_num][chbrclf_mo_num];
|
||||
rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0]));
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
// Test overlap of MO
|
||||
//double point_x[100];
|
||||
//double point_y[100];
|
||||
//double point_z[100];
|
||||
//int32_t npoints=100;
|
||||
//// obtain points
|
||||
//double dr = 20./(npoints-1);
|
||||
//double dr3 = dr*dr*dr;
|
||||
//
|
||||
//for (int i=0;i<npoints;++i) {
|
||||
// point_x[i] = -10. + dr*i;
|
||||
// point_y[i] = -10. + dr*i;
|
||||
// point_z[i] = -10. + dr*i;
|
||||
//}
|
||||
//
|
||||
//double ovlmo1 = 0.0;
|
||||
//// Calculate overlap
|
||||
//for (int i=0;i<npoints;++i) {
|
||||
// printf(".");
|
||||
// fflush(stdout);
|
||||
// for (int j=0;j<npoints;++j) {
|
||||
// for (int k=0;k<npoints;++k) {
|
||||
// // Set point
|
||||
// elec_coord[0] = point_x[i];
|
||||
// elec_coord[1] = point_y[j];
|
||||
// elec_coord[2] = point_z[k];
|
||||
// rc = qmckl_set_electron_coord (context, 'N', elec_coord);
|
||||
// assert(rc == QMCKL_SUCCESS);
|
||||
//
|
||||
// // Calculate value of MO (1st electron)
|
||||
// double mo_vgl[5][walk_num][elec_num][chbrclf_mo_num];
|
||||
// rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0]));
|
||||
// assert (rc == QMCKL_SUCCESS);
|
||||
// ovlmo1 += mo_vgl[0][0][0][0]*mo_vgl[0][0][0][0]*dr3;
|
||||
// }
|
||||
// }
|
||||
//}
|
||||
//printf("OVL MO1 = %10.15f\n",ovlmo1);
|
||||
|
||||
|
||||
printf("\n");
|
||||
printf(" mo_vgl mo_vgl[0][26][219] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][219] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf(" mo_vgl mo_vgl[0][26][220] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][220] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf(" mo_vgl mo_vgl[0][26][221] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][221] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf(" mo_vgl mo_vgl[0][26][222] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][222] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf(" mo_vgl mo_vgl[0][26][223] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[0][2][3]);
|
||||
printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[1][2][3]);
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
#+end_src
|
||||
|
||||
|
||||
* End of files :noexport:
|
||||
|
||||
#+begin_src c :tangle (eval h_private_type)
|
||||
#endif
|
||||
#+end_src
|
||||
|
||||
*** Test
|
||||
#+begin_src c :tangle (eval c_test)
|
||||
rc = qmckl_context_destroy(context);
|
||||
assert (rc == QMCKL_SUCCESS);
|
||||
|
||||
return 0;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Compute file names
|
||||
#+begin_src emacs-lisp
|
||||
; The following is required to compute the file names
|
||||
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq f (concat pwd name "_f.f90"))
|
||||
(setq fh (concat pwd name "_fh.f90"))
|
||||
(setq c (concat pwd name ".c"))
|
||||
(setq h (concat name ".h"))
|
||||
(setq h_private (concat name "_private.h"))
|
||||
(setq c_test (concat pwd "test_" name ".c"))
|
||||
(setq f_test (concat pwd "test_" name "_f.f90"))
|
||||
|
||||
; Minted
|
||||
(require 'ox-latex)
|
||||
(setq org-latex-listings 'minted)
|
||||
(add-to-list 'org-latex-packages-alist '("" "listings"))
|
||||
(add-to-list 'org-latex-packages-alist '("" "color"))
|
||||
|
||||
#+end_src
|
||||
|
||||
|
||||
# -*- mode: org -*-
|
||||
# vim: syntax=c
|
||||
|
||||
|
|
@ -262,7 +262,7 @@ qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int r
|
|||
~qmckl_get_numprec_range~ returns the value of the numerical range in the context.
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
||||
int32_t qmckl_context_get_range(const qmckl_context context);
|
||||
int32_t qmckl_get_numprec_get_range(const qmckl_context context);
|
||||
#+end_src
|
||||
|
||||
# Source
|
||||
|
|
File diff suppressed because one or more lines are too long
58924
org/qmckl_tests.org
58924
org/qmckl_tests.org
File diff suppressed because it is too large
Load Diff
|
@ -6,8 +6,10 @@ qmckl_numprec.org
|
|||
qmckl_nucleus.org
|
||||
qmckl_electron.org
|
||||
qmckl_ao.org
|
||||
qmckl_mo.org
|
||||
qmckl_jastrow.org
|
||||
qmckl_sherman_morrison_woodbury.org
|
||||
qmckl_distance.org
|
||||
qmckl_utils.org
|
||||
qmckl_blas.org
|
||||
qmckl_tests.org
|
||||
|
|
|
@ -9,6 +9,7 @@
|
|||
(setq org-confirm-babel-evaluate nil)
|
||||
(global-font-lock-mode t)
|
||||
(setq org-src-fontify-natively t)
|
||||
(setq python-indent-guess-indent-offset-verbose nil) ;; Remove warning : Can’t guess python-indent-offset
|
||||
|
||||
;(require 'ox-latex)
|
||||
;(setq org-latex-listings t)
|
||||
|
|
|
@ -39,6 +39,7 @@
|
|||
(setq org-hide-emphasis-markers nil)
|
||||
(setq org-pretty-entities nil)
|
||||
(setq org-confirm-babel-evaluate nil) ;; Do not ask for confirmation all the time!!
|
||||
(setq python-indent-guess-indent-offset-verbose nil) ;; Remove warning : Can’t guess python-indent-offset
|
||||
|
||||
(org-babel-do-load-languages
|
||||
'org-babel-load-languages
|
||||
|
|
|
@ -39,6 +39,7 @@
|
|||
f_of_c_d = { '' : ''
|
||||
, 'qmckl_context' : 'integer (c_int64_t)'
|
||||
, 'qmckl_exit_code' : 'integer (c_int32_t)'
|
||||
, 'bool' : 'logical*8'
|
||||
, 'int32_t' : 'integer (c_int32_t)'
|
||||
, 'int64_t' : 'integer (c_int64_t)'
|
||||
, 'uint32_t' : 'integer (c_int32_t)'
|
||||
|
@ -61,6 +62,7 @@ ctypeid_d = { '' : ''
|
|||
, 'real' : 'real(c_float)'
|
||||
, 'real*8' : 'real(c_double)'
|
||||
, 'character' : 'character(c_char)'
|
||||
, 'character' : 'character(c_char)'
|
||||
}
|
||||
#+END_SRC
|
||||
|
||||
|
|
Loading…
Reference in New Issue
Block a user