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

Fixes in Fortran interface

Browse Source
This commit is contained in:
Anthony Scemama 2023-09-11 17:05:41 +02:00
parent c176cd86c3
commit cbc8b9bd58
9 changed files with 231 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
org/qmckl_ao.org
View File

@ -2613,7 +2613,7 @@ qmckl_exit_code qmckl_finalize_basis(qmckl_context context) {
for (int32_t l=0 ; l<32 ; ++l) { for (int32_t l=0 ; l<32 ; ++l) {
lstart[l] = l*(l+1)*(l+2)/6; lstart[l] = l*(l+1)*(l+2)/6;
} }
int64_t ao_idx = 0; int64_t ao_idx = 0;
for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) { for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl]; const int64_t ishell_start = ctx->ao_basis.nucleus_index[inucl];
@ -2976,7 +2976,7 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: primitive_vgl(*) real(c_double), intent(out) :: primitive_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function end function
end interface end interface
@ -3037,7 +3037,7 @@ qmckl_get_ao_basis_shell_vgl (qmckl_context context,
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: shell_vgl(*) real(c_double), intent(out) :: shell_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function end function
end interface end interface
@ -3098,7 +3098,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: ao_vgl(*) real(c_double), intent(out) :: ao_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_ao_basis_ao_vgl end function qmckl_get_ao_basis_ao_vgl
end interface end interface
@ -3164,7 +3164,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: ao_vgl(*) real(c_double), intent(out) :: ao_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_ao_basis_ao_vgl_inplace end function qmckl_get_ao_basis_ao_vgl_inplace
end interface end interface
@ -3225,7 +3225,7 @@ qmckl_get_ao_basis_ao_value (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: ao_value(*) real(c_double) , intent(out) :: ao_value(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_ao_basis_ao_value end function qmckl_get_ao_basis_ao_value
end interface end interface
@ -3291,7 +3291,7 @@ qmckl_get_ao_basis_ao_value_inplace (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: ao_value(*) real(c_double) , intent(out) :: ao_value(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_ao_basis_ao_value_inplace end function qmckl_get_ao_basis_ao_value_inplace
end interface end interface
@ -5791,13 +5791,13 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
{ coord[ipoint], { coord[ipoint],
coord[ipoint + point_num], coord[ipoint + point_num],
coord[ipoint + 2*point_num] }; coord[ipoint + 2*point_num] };
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) { for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
const double n_coord[3] __attribute__((aligned(64))) = const double n_coord[3] __attribute__((aligned(64))) =
{ nucl_coord[inucl], { nucl_coord[inucl],
nucl_coord[inucl + nucl_num], nucl_coord[inucl + nucl_num],
nucl_coord[inucl + 2*nucl_num] }; nucl_coord[inucl + 2*nucl_num] };
/* Test if the point is in the range of the nucleus */ /* Test if the point is in the range of the nucleus */
const double x = e_coord[0] - n_coord[0]; const double x = e_coord[0] - n_coord[0];
const double y = e_coord[1] - n_coord[1]; const double y = e_coord[1] - n_coord[1];
@ -6597,7 +6597,7 @@ qmckl_compute_ao_vgl_hpc_gaussian (
{ coord[ipoint], { coord[ipoint],
coord[ipoint + point_num], coord[ipoint + point_num],
coord[ipoint + 2*point_num] }; coord[ipoint + 2*point_num] };
for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) { for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) {
const double n_coord[3] __attribute__((aligned(64))) = const double n_coord[3] __attribute__((aligned(64))) =
{ nucl_coord[inucl], { nucl_coord[inucl],

16
org/qmckl_blas.org
View File

@ -1098,8 +1098,8 @@ end function qmckl_dgemm_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: TransA character(c_char) , intent(in) , value :: TransA
character , intent(in) , value :: TransB character(c_char) , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k integer (c_int64_t) , intent(in) , value :: k
@ -1133,8 +1133,8 @@ end function qmckl_dgemm_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: TransA character(c_char) , intent(in) , value :: TransA
character , intent(in) , value :: TransB character(c_char) , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k integer (c_int64_t) , intent(in) , value :: k
@ -1386,8 +1386,8 @@ end function qmckl_dgemm_safe_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: TransA character(c_char) , intent(in) , value :: TransA
character , intent(in) , value :: TransB character(c_char) , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k integer (c_int64_t) , intent(in) , value :: k
@ -1424,8 +1424,8 @@ end function qmckl_dgemm_safe_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: TransA character(c_char) , intent(in) , value :: TransA
character , intent(in) , value :: TransB character(c_char) , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k integer (c_int64_t) , intent(in) , value :: k

32
org/qmckl_distance.org
View File

@ -236,8 +236,8 @@ end function qmckl_distance_sq_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -267,8 +267,8 @@ end function qmckl_distance_sq_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -649,8 +649,8 @@ end function qmckl_distance_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -680,8 +680,8 @@ end function qmckl_distance_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -1076,8 +1076,8 @@ end function qmckl_distance_rescaled_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -1108,8 +1108,8 @@ end function qmckl_distance_rescaled_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -1599,8 +1599,8 @@ end function qmckl_distance_rescaled_gl_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)
@ -1631,8 +1631,8 @@ end function qmckl_distance_rescaled_gl_f
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transa character(c_char) , intent(in) , value :: transa
character , intent(in) , value :: transb character(c_char) , intent(in) , value :: transb
integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
real (c_double ) , intent(in) :: A(lda,*) real (c_double ) , intent(in) :: A(lda,*)

10
org/qmckl_electron.org
View File

@ -112,7 +112,7 @@ typedef struct qmckl_walker_struct {
int64_t num; int64_t num;
qmckl_point_struct point; qmckl_point_struct point;
} qmckl_walker; } qmckl_walker;
typedef struct qmckl_electron_struct { typedef struct qmckl_electron_struct {
int64_t num; int64_t num;
int64_t up_num; int64_t up_num;
@ -287,7 +287,7 @@ qmckl_set_electron_coord(qmckl_context context,
{ {
int32_t mask = 0; int32_t mask = 0;
<<pre2>> <<pre2>>
if (transp != 'N' && transp != 'T') { if (transp != 'N' && transp != 'T') {
@ -347,9 +347,9 @@ interface
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
character , intent(in) , value :: transp character(c_char) , intent(in) , value :: transp
integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: walk_num
double precision , intent(in) :: coord(*) real(c_double) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function end function
end interface end interface
@ -723,7 +723,7 @@ qmckl_exit_code qmckl_provide_ee_distance(qmckl_context context)
free(ctx->electron.ee_distance); free(ctx->electron.ee_distance);
ctx->electron.ee_distance = NULL; ctx->electron.ee_distance = NULL;
} }
/* Allocate array */ /* Allocate array */
if (ctx->electron.ee_distance == NULL) { if (ctx->electron.ee_distance == NULL) {

18
org/qmckl_error.org
View File

@ -362,7 +362,7 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
import import
implicit none implicit none
integer (qmckl_exit_code), intent(in), value :: error integer (qmckl_exit_code), intent(in), value :: error
character, intent(out) :: string(<<MAX_STRING_LENGTH()>>) character(c_char), intent(out) :: string(<<MAX_STRING_LENGTH()>>)
end subroutine qmckl_string_of_error end subroutine qmckl_string_of_error
end interface end interface
#+end_src #+end_src
@ -440,7 +440,7 @@ qmckl_set_error(qmckl_context context,
Upon error, the error type and message can be obtained from the Upon error, the error type and message can be obtained from the
context using ~qmckl_get_error~. The message and function name context using ~qmckl_get_error~. The message and function name
is returned in the variables provided. Therefore, passing a is returned in the variables provided. Therefore, passing a
function name and message is mandatory. function name and message is mandatory.
# Header # Header
@ -494,7 +494,7 @@ qmckl_get_error(qmckl_context context,
#+end_src #+end_src
* Failing * Failing
To make a function fail, the ~qmckl_failwith~ function should be To make a function fail, the ~qmckl_failwith~ function should be
called, such that information about the failure is stored in called, such that information about the failure is stored in
the context. The desired exit code is given as an argument, as the context. The desired exit code is given as an argument, as
@ -511,7 +511,7 @@ qmckl_failwith(qmckl_context context,
const char* function, const char* function,
const char* message) ; const char* message) ;
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) :exports none #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_failwith(qmckl_context context, qmckl_failwith(qmckl_context context,
@ -593,7 +593,7 @@ qmckl_last_error(qmckl_context context, char* buffer) {
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
#+end_src #+end_src
** Fortran inteface ** Fortran inteface
#+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes
@ -603,7 +603,7 @@ qmckl_last_error(qmckl_context context, char* buffer) {
import import
implicit none implicit none
integer (c_int64_t) , intent(in), value :: context integer (c_int64_t) , intent(in), value :: context
character, intent(out) :: string(*) character(c_char), intent(out) :: string(*)
end subroutine qmckl_last_error end subroutine qmckl_last_error
end interface end interface
#+end_src #+end_src
@ -625,7 +625,7 @@ qmckl_check(qmckl_context context, qmckl_exit_code rc);
qmckl_exit_code qmckl_exit_code
qmckl_check(qmckl_context context, qmckl_exit_code rc) qmckl_check(qmckl_context context, qmckl_exit_code rc)
{ {
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
char fname[QMCKL_MAX_FUN_LEN]; char fname[QMCKL_MAX_FUN_LEN];
@ -639,7 +639,7 @@ qmckl_check(qmckl_context context, qmckl_exit_code rc)
return rc; return rc;
} }
#+end_src #+end_src
It should be used as: It should be used as:
#+begin_src c #+begin_src c
@ -648,7 +648,7 @@ rc = qmckl_check(context,
); );
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
#+end_src #+end_src
** Fortran inteface ** Fortran inteface
#+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes

162
org/qmckl_examples.org
View File

@ -1,31 +1,72 @@
#+TITLE: Code examples #+TITLE: Code examples
#+SETUPFILE: ../tools/theme.setup #+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org #+INCLUDE: ../tools/lib.org
In this section, we present examples of usage of QMCkl.
For simplicity, we assume that the wave function parameters are stored
in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
* Python In this section, we provide hands-on examples to demonstrate the capabilities
** Check numerically that MOs are orthonormal of the QMCkl library. We furnish code samples in C, Fortran, and Python,
serving as exhaustive tutorials for effectively leveraging QMCkl.
For simplicity, we assume that the wave function parameters are stored in a
[[https://github.com/TREX-CoE/trexio][TREXIO]] file.
* Overlap matrix in the MO basis
The focal point of this example is the numerical evaluation of the overlap
matrix in the MO basis. Utilizing QMCkl, we approximate the MOs at
discrete grid points to compute the overlap matrix \( S_{ij} \) as follows:
\[
S_{ij} = \int \phi_i(\mathbf{r})\, \phi_j(\mathbf{r}) \text{d}\mathbf{r} \approx
\sum_k \phi_i(\mathbf{r}_k)\, \phi_j(\mathbf{r}_k) \delta\mathbf{r}
\]
The code starts by reading a wave function from a TREXIO file. This is
accomplished using the ~qmckl_trexio_read~ function, which populates a
~qmckl_context~ with the necessary wave function parameters. The context
serves as the primary interface for interacting with the QMCkl library,
encapsulating the state and configurations of the system.
Subsequently, the code retrieves various attributes such as the number of
nuclei ~nucl_num~ and coordinates ~nucl_coord~.
These attributes are essential for setting up the integration grid.
The core of the example lies in the numerical computation of the overlap
matrix. To achieve this, the code employs a regular grid in three-dimensional
space, and the grid points are then populated into the ~qmckl_context~ using
the ~qmckl_set_point~ function.
The MO values at these grid points are computed using the
~qmckl_get_mo_basis_mo_value~ function. These values are then used to
calculate the overlap matrix through a matrix multiplication operation
facilitated by the ~qmckl_dgemm~ function.
The code is also instrumented to measure the execution time for the MO
value computation, providing an empirical assessment of the computational
efficiency. Error handling is robustly implemented at each stage to ensure the
reliability of the simulation.
In summary, this example serves as a holistic guide for leveraging the QMCkl
library, demonstrating its ease of use. It provides a concrete starting point
for researchers and developers interested in integrating QMCkl into their QMC
code.
** Python
:PROPERTIES: :PROPERTIES:
:header-args: :tangle mo_ortho.py :header-args: :tangle mo_ortho.py
:END: :END:
In this example, we will compute numerically the overlap In this example, we will compute numerically the overlap
between the molecular orbitals: between the molecular orbitals:
\[ \[
S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r}) S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
\text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k) \text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
\phi_j(\mathbf{r}_k) \delta \mathbf{r} \phi_j(\mathbf{r}_k) \delta \mathbf{r}
\] \]
\[ \[
S_{ij} = \langle \phi_i | \phi_j \rangle S_{ij} = \langle \phi_i | \phi_j \rangle
\sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle \sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
\langle \mathbf{r}_k | \phi_j \rangle \langle \mathbf{r}_k | \phi_j \rangle
\] \]
#+begin_src python :exports code #+begin_src python :exports code
import numpy as np import numpy as np
@ -33,11 +74,11 @@ import qmckl
#+end_src #+end_src
#+RESULTS: #+RESULTS:
First, we create a context for the QMCkl calculation, and load the First, we create a context for the QMCkl calculation, and load the
wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock
determinant for the water molecule in the cc-pV5Z basis set. determinant for the water molecule in the cc-pV5Z basis set.
#+begin_src python :exports code #+begin_src python :exports code
trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5" trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5"
@ -52,7 +93,7 @@ qmckl.trexio_read(context, trexio_filename)
molecule. molecule.
We fetch the nuclear coordinates from the context, We fetch the nuclear coordinates from the context,
#+begin_src python :exports code #+begin_src python :exports code
nucl_num = qmckl.get_nucleus_num(context) nucl_num = qmckl.get_nucleus_num(context)
@ -75,7 +116,7 @@ for i in range(nucl_num):
#+end_example #+end_example
and compute the coordinates of the grid points: and compute the coordinates of the grid points:
#+begin_src python :exports code #+begin_src python :exports code
nx = ( 120, 120, 120 ) nx = ( 120, 120, 120 )
shift = np.array([5.,5.,5.]) shift = np.array([5.,5.,5.])
@ -97,7 +138,7 @@ dr = step[0] * step[1] * step[2]
#+end_src #+end_src
#+RESULTS: #+RESULTS:
Now the grid is ready, we can create the list of grid points Now the grid is ready, we can create the list of grid points
$\mathbf{r}_k$ on which the MOs $\phi_i$ will be evaluated, and $\mathbf{r}_k$ on which the MOs $\phi_i$ will be evaluated, and
transfer them to the QMCkl context: transfer them to the QMCkl context:
@ -116,7 +157,7 @@ qmckl.set_point(context, 'N', point_num, np.reshape(point, (point_num*3)))
#+RESULTS: #+RESULTS:
: None : None
Then, we evaluate all the MOs at the grid points (and time the execution), Then, we evaluate all the MOs at the grid points (and time the execution),
and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i \rangle = and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i \rangle =
\phi_i(\mathbf{r}_k)$. \phi_i(\mathbf{r}_k)$.
@ -168,26 +209,27 @@ print (overlap)
1.18264754e-09 8.97215950e-01]] 1.18264754e-09 8.97215950e-01]]
#+end_example #+end_example
* C ** C
** Check numerically that MOs are orthonormal, with cusp fitting In this example, electron-nucleus cusp fitting is added.
:PROPERTIES: :PROPERTIES:
:header-args: :tangle mo_ortho.c :header-args: :tangle mo_ortho.c
:END: :END:
In this example, we will compute numerically the overlap In this example, we will compute numerically the overlap
between the molecular orbitals: between the molecular orbitals:
\[ \[
S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r}) S_{ij} = \int \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
\text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k) \text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
\phi_j(\mathbf{r}_k) \delta \mathbf{r} \phi_j(\mathbf{r}_k) \delta \mathbf{r}
\] \]
\[ \[
S_{ij} = \langle \phi_i | \phi_j \rangle S_{ij} = \langle \phi_i | \phi_j \rangle
\sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle \sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
\langle \mathbf{r}_k | \phi_j \rangle \langle \mathbf{r}_k | \phi_j \rangle
\] \]
We apply the cusp fitting procedure, so the MOs might deviate We apply the cusp fitting procedure, so the MOs might deviate
slightly from orthonormality. slightly from orthonormality.
@ -206,7 +248,7 @@ int main(int argc, char** argv)
First, we create a context for the QMCkl calculation, and load the First, we create a context for the QMCkl calculation, and load the
wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock wave function stored in =h2o_5z.h5= inside it. It is a Hartree-Fock
determinant for the water molecule in the cc-pV5Z basis set. determinant for the water molecule in the cc-pV5Z basis set.
#+begin_src c :exports code #+begin_src c :exports code
qmckl_context context = qmckl_context_create(); qmckl_context context = qmckl_context_create();
@ -251,18 +293,18 @@ int main(int argc, char** argv)
for (size_t i=0 ; i<nucl_num ; ++i) { for (size_t i=0 ; i<nucl_num ; ++i) {
switch ((int) nucl_charge[i]) { switch ((int) nucl_charge[i]) {
case 1: case 1:
r_cusp[i] = 0.5; r_cusp[i] = 0.5;
break; break;
case 8: case 8:
r_cusp[i] = 0.1; r_cusp[i] = 0.1;
break; break;
} }
} }
rc = qmckl_set_mo_basis_r_cusp(context, &(r_cusp[0]), nucl_num); rc = qmckl_set_mo_basis_r_cusp(context, &(r_cusp[0]), nucl_num);
@ -271,7 +313,7 @@ int main(int argc, char** argv)
exit(1); exit(1);
} }
#+end_src #+end_src
@ -305,7 +347,7 @@ int main(int argc, char** argv)
#+end_example #+end_example
and compute the coordinates of the grid points: and compute the coordinates of the grid points:
#+begin_src c :exports code #+begin_src c :exports code
size_t nx[3] = { 120, 120, 120 }; size_t nx[3] = { 120, 120, 120 };
double shift[3] = {5.,5.,5.}; double shift[3] = {5.,5.,5.};
@ -316,8 +358,8 @@ int main(int argc, char** argv)
for (size_t i=0 ; i<nucl_num ; ++i) { for (size_t i=0 ; i<nucl_num ; ++i) {
for (int j=0 ; j<3 ; ++j) { for (int j=0 ; j<3 ; ++j) {
rmin[j] = nucl_coord[i][j] < rmin[j] ? nucl_coord[i][j] : rmin[j]; rmin[j] = nucl_coord[i][j] < rmin[j] ? nucl_coord[i][j] : rmin[j];
rmax[j] = nucl_coord[i][j] > rmax[j] ? nucl_coord[i][j] : rmax[j]; rmax[j] = nucl_coord[i][j] > rmax[j] ? nucl_coord[i][j] : rmax[j];
} }
} }
@ -335,7 +377,7 @@ int main(int argc, char** argv)
fprintf(stderr, "Allocation failed (linspace)\n"); fprintf(stderr, "Allocation failed (linspace)\n");
exit(1); exit(1);
} }
step[i] = (rmax[i] - rmin[i]) / ((double) (nx[i]-1)); step[i] = (rmax[i] - rmin[i]) / ((double) (nx[i]-1));
for (size_t j=0 ; j<nx[i] ; ++j) { for (size_t j=0 ; j<nx[i] ; ++j) {
@ -343,7 +385,7 @@ int main(int argc, char** argv)
} }
} }
double dr = step[0] * step[1] * step[2]; double dr = step[0] * step[1] * step[2];
#+end_src #+end_src
@ -363,7 +405,7 @@ int main(int argc, char** argv)
for (size_t i=0 ; i<nx[0] ; ++i) { for (size_t i=0 ; i<nx[0] ; ++i) {
for (size_t j=0 ; j<nx[1] ; ++j) { for (size_t j=0 ; j<nx[1] ; ++j) {
for (size_t k=0 ; k<nx[2] ; ++k) { for (size_t k=0 ; k<nx[2] ; ++k) {
point[m] = linspace[0][i]; point[m] = linspace[0][i];
m++; m++;
@ -372,7 +414,7 @@ int main(int argc, char** argv)
point[m] = linspace[2][k]; point[m] = linspace[2][k];
m++; m++;
} }
} }
} }
@ -387,7 +429,7 @@ int main(int argc, char** argv)
#+RESULTS: #+RESULTS:
: None : None
Then, we evaluate all the MOs at the grid points (and time the execution), Then, we evaluate all the MOs at the grid points (and time the execution),
and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i and thus obtain the matrix $M_{ki} = \langle \mathbf{r}_k | \phi_i
\rangle = \phi_i(\mathbf{r}_k)$. \rangle = \phi_i(\mathbf{r}_k)$.
@ -405,7 +447,7 @@ int main(int argc, char** argv)
fprintf(stderr, "Allocation failed (mo_value)\n"); fprintf(stderr, "Allocation failed (mo_value)\n");
exit(1); exit(1);
} }
gettimeofday(&timecheck, NULL); gettimeofday(&timecheck, NULL);
before = (long)timecheck.tv_sec * 1000 + (long)timecheck.tv_usec / 1000; before = (long)timecheck.tv_sec * 1000 + (long)timecheck.tv_usec / 1000;
@ -439,7 +481,7 @@ Execution time : 5.608000 seconds
rc = qmckl_dgemm(context, 'N', 'T', mo_num, mo_num, point_num, dr, rc = qmckl_dgemm(context, 'N', 'T', mo_num, mo_num, point_num, dr,
mo_value, mo_num, mo_value, mo_num, 0.0, mo_value, mo_num, mo_value, mo_num, 0.0,
overlap, mo_num); overlap, mo_num);
for (size_t i=0 ; i<mo_num ; ++i) { for (size_t i=0 ; i<mo_num ; ++i) {
printf("%4ld", i); printf("%4ld", i);
for (size_t j=0 ; j<mo_num ; ++j) { for (size_t j=0 ; j<mo_num ; ++j) {
@ -452,7 +494,7 @@ Execution time : 5.608000 seconds
#+end_src #+end_src
#+begin_example #+begin_example
0 0.988765 0.002336 0.000000 -0.000734 0.000000 0.000530 0.000000 0.000446 0.000000 -0.000000 -0.000511 -0.000000 -0.000267 0.000000 0.000000 0.001007 0.000000 0.000168 -0.000000 -0.000000 -0.000670 -0.000000 0.000000 -0.000251 -0.000261 -0.000000 -0.000000 -0.000000 -0.000397 -0.000000 -0.000810 0.000000 0.000231 -0.000000 -0.000000 0.000000 -0.000000 0 0.988765 0.002336 0.000000 -0.000734 0.000000 0.000530 0.000000 0.000446 0.000000 -0.000000 -0.000511 -0.000000 -0.000267 0.000000 0.000000 0.001007 0.000000 0.000168 -0.000000 -0.000000 -0.000670 -0.000000 0.000000 -0.000251 -0.000261 -0.000000 -0.000000 -0.000000 -0.000397 -0.000000 -0.000810 0.000000 0.000231 -0.000000 -0.000000 0.000000 -0.000000
... ...
200 0.039017 -0.013066 -0.000000 -0.001935 -0.000000 -0.003829 -0.000000 0.000996 -0.000000 0.000000 -0.003733 0.000000 0.000065 -0.000000 -0.000000 -0.002220 -0.000000 -0.001961 0.000000 0.000000 -0.004182 0.000000 -0.000000 -0.000165 -0.002445 0.000000 -0.000000 0.000000 0.001985 0.000000 0.001685 -0.000000 -0.002899 0.000000 0.000000 0.000000 -0.000000 0.002591 0.000000 -0.000000 0.000000 0.002385 0.000000 -0.011086 0.000000 -0.003885 0.000000 -0.000000 0.003602 -0.000000 0.000000 -0.003241 0.000000 0.000000 0.002613 -0.007344 -0.000000 -0.000000 0.000000 0.000017 0.000000 0.000000 0.000000 -0.008719 0.000000 -0.001358 -0.003233 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.003778 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.001190 0.000000 0.000000 -0.000000 0.005578 -0.000000 -0.001502 0.003899 -0.000000 -0.000000 0.000000 -0.003291 -0.001775 -0.000000 -0.002374 0.000000 -0.000000 -0.000000 -0.000000 -0.001290 -0.000000 0.002178 0.000000 0.000000 0.000000 -0.001252 0.000000 -0.000000 -0.000926 0.000000 -0.000000 -0.013130 -0.000000 0.012124 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.025194 0.000343 -0.000000 0.000000 -0.000000 -0.004421 0.000000 0.000000 -0.000599 -0.000000 0.005289 0.000000 -0.000000 0.012826 -0.000000 0.000000 0.006190 0.000000 0.000000 -0.000000 0.000000 -0.000321 0.000000 -0.000000 -0.000000 0.000000 -0.000000 0.001499 -0.006629 0.000000 0.000000 0.000000 -0.000000 0.008737 -0.000000 0.006880 0.000000 -0.001851 -0.000000 -0.000000 0.000000 -0.007464 0.000000 0.010298 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.000540 0.000000 -0.006616 -0.000000 0.000000 -0.002927 -0.000000 0.000000 0.010352 0.000000 -0.003103 -0.000000 -0.007640 -0.000000 -0.000000 0.005302 0.000000 0.000000 -0.000000 -0.000000 -0.010181 0.000000 -0.001108 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.000998 -0.009754 0.013562 0.000000 -0.000000 0.887510 200 0.039017 -0.013066 -0.000000 -0.001935 -0.000000 -0.003829 -0.000000 0.000996 -0.000000 0.000000 -0.003733 0.000000 0.000065 -0.000000 -0.000000 -0.002220 -0.000000 -0.001961 0.000000 0.000000 -0.004182 0.000000 -0.000000 -0.000165 -0.002445 0.000000 -0.000000 0.000000 0.001985 0.000000 0.001685 -0.000000 -0.002899 0.000000 0.000000 0.000000 -0.000000 0.002591 0.000000 -0.000000 0.000000 0.002385 0.000000 -0.011086 0.000000 -0.003885 0.000000 -0.000000 0.003602 -0.000000 0.000000 -0.003241 0.000000 0.000000 0.002613 -0.007344 -0.000000 -0.000000 0.000000 0.000017 0.000000 0.000000 0.000000 -0.008719 0.000000 -0.001358 -0.003233 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.003778 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.001190 0.000000 0.000000 -0.000000 0.005578 -0.000000 -0.001502 0.003899 -0.000000 -0.000000 0.000000 -0.003291 -0.001775 -0.000000 -0.002374 0.000000 -0.000000 -0.000000 -0.000000 -0.001290 -0.000000 0.002178 0.000000 0.000000 0.000000 -0.001252 0.000000 -0.000000 -0.000926 0.000000 -0.000000 -0.013130 -0.000000 0.012124 0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.025194 0.000343 -0.000000 0.000000 -0.000000 -0.004421 0.000000 0.000000 -0.000599 -0.000000 0.005289 0.000000 -0.000000 0.012826 -0.000000 0.000000 0.006190 0.000000 0.000000 -0.000000 0.000000 -0.000321 0.000000 -0.000000 -0.000000 0.000000 -0.000000 0.001499 -0.006629 0.000000 0.000000 0.000000 -0.000000 0.008737 -0.000000 0.006880 0.000000 -0.001851 -0.000000 -0.000000 0.000000 -0.007464 0.000000 0.010298 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 0.000000 0.000540 0.000000 -0.006616 -0.000000 0.000000 -0.002927 -0.000000 0.000000 0.010352 0.000000 -0.003103 -0.000000 -0.007640 -0.000000 -0.000000 0.005302 0.000000 0.000000 -0.000000 -0.000000 -0.010181 0.000000 -0.001108 0.000000 0.000000 -0.000000 0.000000 0.000000 -0.000998 -0.009754 0.013562 0.000000 -0.000000 0.887510
#+end_example #+end_example
@ -465,7 +507,7 @@ Execution time : 5.608000 seconds
error in text format and exits the program. error in text format and exits the program.
#+NAME: qmckl_check_error #+NAME: qmckl_check_error
#+begin_src f90 #+begin_src f90
subroutine qmckl_check_error(rc, message) subroutine qmckl_check_error(rc, message)
use qmckl use qmckl
implicit none implicit none
@ -480,7 +522,7 @@ subroutine qmckl_check_error(rc, message)
end if end if
end subroutine qmckl_check_error end subroutine qmckl_check_error
#+end_src #+end_src
** Computing an atomic orbital on a grid ** Computing an atomic orbital on a grid
:PROPERTIES: :PROPERTIES:
:header-args: :tangle ao_grid.f90 :header-args: :tangle ao_grid.f90
@ -492,14 +534,14 @@ end subroutine qmckl_check_error
atomic units in the borders. atomic units in the borders.
This program uses the ~qmckl_check_error~ function defined above. This program uses the ~qmckl_check_error~ function defined above.
To use this program, run To use this program, run
#+begin_src bash :tangle no :exports code #+begin_src bash :tangle no :exports code
$ ao_grid <trexio_file> <AO_id> <point_num> $ ao_grid <trexio_file> <AO_id> <point_num>
#+end_src #+end_src
#+begin_src f90 :noweb yes #+begin_src f90 :noweb yes
<<qmckl_check_error>> <<qmckl_check_error>>
@ -529,7 +571,7 @@ program ao_grid
Start by fetching the command-line arguments: Start by fetching the command-line arguments:
#+begin_src f90 #+begin_src f90
if (iargc() /= 3) then if (iargc() /= 3) then
print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>' print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>'
call exit(-1) call exit(-1)
@ -549,7 +591,7 @@ program ao_grid
Create the QMCkl context and initialize it with the wave function Create the QMCkl context and initialize it with the wave function
present in the TREXIO file: present in the TREXIO file:
#+begin_src f90 #+begin_src f90
qmckl_ctx = qmckl_context_create() qmckl_ctx = qmckl_context_create()
rc = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename))) rc = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename)))
call qmckl_check_error(rc, 'Read TREXIO') call qmckl_check_error(rc, 'Read TREXIO')
@ -557,8 +599,8 @@ program ao_grid
We need to check that ~ao_id~ is in the range, so we get the total We need to check that ~ao_id~ is in the range, so we get the total
number of AOs from QMCkl: number of AOs from QMCkl:
#+begin_src f90 #+begin_src f90
rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num) rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num)
call qmckl_check_error(rc, 'Getting ao_num') call qmckl_check_error(rc, 'Getting ao_num')
@ -571,7 +613,7 @@ program ao_grid
Now we will compute the limits of the box in which the molecule fits. Now we will compute the limits of the box in which the molecule fits.
For that, we first need to ask QMCkl the coordinates of nuclei. For that, we first need to ask QMCkl the coordinates of nuclei.
#+begin_src f90 #+begin_src f90
rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num) rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num)
call qmckl_check_error(rc, 'Get nucleus num') call qmckl_check_error(rc, 'Get nucleus num')
@ -583,11 +625,11 @@ program ao_grid
We now compute the coordinates of opposite points of the box, and We now compute the coordinates of opposite points of the box, and
the distance between points along the 3 directions: the distance between points along the 3 directions:
#+begin_src f90 #+begin_src f90
rmin(1) = minval( nucl_coord(1,:) ) - 5.d0 rmin(1) = minval( nucl_coord(1,:) ) - 5.d0
rmin(2) = minval( nucl_coord(2,:) ) - 5.d0 rmin(2) = minval( nucl_coord(2,:) ) - 5.d0
rmin(3) = minval( nucl_coord(3,:) ) - 5.d0 rmin(3) = minval( nucl_coord(3,:) ) - 5.d0
rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0 rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0
rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0 rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0
rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0 rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0
@ -597,8 +639,8 @@ program ao_grid
We now produce the list of point coordinates where the AO will be We now produce the list of point coordinates where the AO will be
evaluated: evaluated:
#+begin_src f90 #+begin_src f90
point_num = point_num_x**3 point_num = point_num_x**3
allocate( points(point_num, 3) ) allocate( points(point_num, 3) )
ipoint=0 ipoint=0
@ -619,19 +661,19 @@ program ao_grid
z = z + dr(3) z = z + dr(3)
end do end do
#+end_src #+end_src
We give the points to QMCkl: We give the points to QMCkl:
#+begin_src f90 #+begin_src f90
rc = qmckl_set_point(qmckl_ctx, 'T', point_num, points, size(points)*1_8 ) rc = qmckl_set_point(qmckl_ctx, 'T', point_num, points, size(points)*1_8 )
call qmckl_check_error(rc, 'Setting points') call qmckl_check_error(rc, 'Setting points')
#+end_src #+end_src
We allocate the space required to retrieve the values, gradients and We allocate the space required to retrieve the values, gradients and
Laplacian of all AOs, and ask to retrieve the values of the Laplacian of all AOs, and ask to retrieve the values of the
AOs computed at the point positions. AOs computed at the point positions.
#+begin_src f90 #+begin_src f90
allocate( ao_vgl(ao_num, 5, point_num) ) allocate( ao_vgl(ao_num, 5, point_num) )
rc = qmckl_get_ao_basis_ao_vgl(qmckl_ctx, ao_vgl, ao_num*5_8*point_num) rc = qmckl_get_ao_basis_ao_vgl(qmckl_ctx, ao_vgl, ao_num*5_8*point_num)
call qmckl_check_error(rc, 'Setting points') call qmckl_check_error(rc, 'Setting points')
@ -639,13 +681,13 @@ program ao_grid
We finally print the value and Laplacian of the AO: We finally print the value and Laplacian of the AO:
#+begin_src f90 #+begin_src f90
do ipoint=1, point_num do ipoint=1, point_num
print '(3(F10.6,X),2(E20.10,X))', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint), ao_vgl(ao_id,5,ipoint) print '(3(F10.6,X),2(E20.10,X))', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint), ao_vgl(ao_id,5,ipoint)
end do end do
#+end_src #+end_src
#+begin_src f90 #+begin_src f90
deallocate( nucl_coord, points, ao_vgl ) deallocate( nucl_coord, points, ao_vgl )
end program ao_grid end program ao_grid
#+end_src #+end_src

40
org/qmckl_jastrow_champ.org
View File

@ -957,7 +957,7 @@ interface
import import
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
double precision, intent(in), value :: kappa_ee real(c_double), intent(in), value :: kappa_ee
end function qmckl_set_jastrow_champ_rescale_factor_ee end function qmckl_set_jastrow_champ_rescale_factor_ee
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_rescale_factor_en (context, & integer(qmckl_exit_code) function qmckl_set_jastrow_champ_rescale_factor_en (context, &
@ -967,7 +967,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(in) :: kappa_en(size_max) real(c_double), intent(in) :: kappa_en(size_max)
end function qmckl_set_jastrow_champ_rescale_factor_en end function qmckl_set_jastrow_champ_rescale_factor_en
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_aord_num (context, & integer(qmckl_exit_code) function qmckl_set_jastrow_champ_aord_num (context, &
@ -1023,7 +1023,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(in) :: a_vector(size_max) real(c_double), intent(in) :: a_vector(size_max)
end function qmckl_set_jastrow_champ_a_vector end function qmckl_set_jastrow_champ_a_vector
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_b_vector(context, & integer(qmckl_exit_code) function qmckl_set_jastrow_champ_b_vector(context, &
@ -1033,7 +1033,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(in) :: b_vector(size_max) real(c_double), intent(in) :: b_vector(size_max)
end function qmckl_set_jastrow_champ_b_vector end function qmckl_set_jastrow_champ_b_vector
integer(qmckl_exit_code) function qmckl_set_jastrow_champ_c_vector(context, & integer(qmckl_exit_code) function qmckl_set_jastrow_champ_c_vector(context, &
@ -1043,7 +1043,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(in) :: c_vector(size_max) real(c_double), intent(in) :: c_vector(size_max)
end function qmckl_set_jastrow_champ_c_vector end function qmckl_set_jastrow_champ_c_vector
end interface end interface
@ -1449,7 +1449,7 @@ interface
import import
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
double precision, intent(out) :: kappa_ee real(c_double), intent(out) :: kappa_ee
end function qmckl_get_jastrow_champ_rescale_factor_ee end function qmckl_get_jastrow_champ_rescale_factor_ee
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_rescale_factor_en (context, & integer(qmckl_exit_code) function qmckl_get_jastrow_champ_rescale_factor_en (context, &
@ -1459,7 +1459,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: kappa_en(size_max) real(c_double), intent(out) :: kappa_en(size_max)
end function qmckl_get_jastrow_champ_rescale_factor_en end function qmckl_get_jastrow_champ_rescale_factor_en
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_aord_num (context, & integer(qmckl_exit_code) function qmckl_get_jastrow_champ_aord_num (context, &
@ -1515,7 +1515,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: a_vector(size_max) real(c_double), intent(out) :: a_vector(size_max)
end function qmckl_get_jastrow_champ_a_vector end function qmckl_get_jastrow_champ_a_vector
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_b_vector(context, & integer(qmckl_exit_code) function qmckl_get_jastrow_champ_b_vector(context, &
@ -1525,7 +1525,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: b_vector(size_max) real(c_double), intent(out) :: b_vector(size_max)
end function qmckl_get_jastrow_champ_b_vector end function qmckl_get_jastrow_champ_b_vector
integer(qmckl_exit_code) function qmckl_get_jastrow_champ_c_vector(context, & integer(qmckl_exit_code) function qmckl_get_jastrow_champ_c_vector(context, &
@ -1535,7 +1535,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in) , value :: context integer (qmckl_context) , intent(in) , value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: c_vector(size_max) real(c_double), intent(out) :: c_vector(size_max)
end function qmckl_get_jastrow_champ_c_vector end function qmckl_get_jastrow_champ_c_vector
end interface end interface
@ -1726,7 +1726,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: asymp_jasb(size_max) real(c_double), intent(out) :: asymp_jasb(size_max)
end function qmckl_get_jastrow_champ_asymp_jasb end function qmckl_get_jastrow_champ_asymp_jasb
end interface end interface
#+end_src #+end_src
@ -2123,7 +2123,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_ee(size_max) real(c_double), intent(out) :: factor_ee(size_max)
end function qmckl_get_jastrow_champ_factor_ee end function qmckl_get_jastrow_champ_factor_ee
end interface end interface
#+end_src #+end_src
@ -2589,7 +2589,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_ee_gl(size_max) real(c_double), intent(out) :: factor_ee_gl(size_max)
end function qmckl_get_jastrow_champ_factor_ee_gl end function qmckl_get_jastrow_champ_factor_ee_gl
end interface end interface
#+end_src #+end_src
@ -3713,7 +3713,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: asymp_jasa(size_max) real(c_double), intent(out) :: asymp_jasa(size_max)
end function qmckl_get_jastrow_champ_asymp_jasa end function qmckl_get_jastrow_champ_asymp_jasa
end interface end interface
#+end_src #+end_src
@ -3973,7 +3973,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_en(size_max) real(c_double), intent(out) :: factor_en(size_max)
end function qmckl_get_jastrow_champ_factor_en end function qmckl_get_jastrow_champ_factor_en
end interface end interface
#+end_src #+end_src
@ -4370,7 +4370,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_en_gl(size_max) real(c_double), intent(out) :: factor_en_gl(size_max)
end function qmckl_get_jastrow_champ_factor_en_gl end function qmckl_get_jastrow_champ_factor_en_gl
end interface end interface
#+end_src #+end_src
@ -8448,7 +8448,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_een(size_max) real(c_double), intent(out) :: factor_een(size_max)
end function qmckl_get_jastrow_champ_factor_een end function qmckl_get_jastrow_champ_factor_een
end interface end interface
#+end_src #+end_src
@ -8982,7 +8982,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: factor_een_gl(size_max) real(c_double), intent(out) :: factor_een_gl(size_max)
end function qmckl_get_jastrow_champ_factor_een_gl end function qmckl_get_jastrow_champ_factor_een_gl
end interface end interface
#+end_src #+end_src
@ -9792,7 +9792,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: value(size_max) real(c_double), intent(out) :: value(size_max)
end function qmckl_get_jastrow_champ_value end function qmckl_get_jastrow_champ_value
end interface end interface
#+end_src #+end_src
@ -10150,7 +10150,7 @@ interface
implicit none implicit none
integer (qmckl_context) , intent(in), value :: context integer (qmckl_context) , intent(in), value :: context
integer(c_int64_t), intent(in), value :: size_max integer(c_int64_t), intent(in), value :: size_max
double precision, intent(out) :: gl(size_max) real(c_double), intent(out) :: gl(size_max)
end function qmckl_get_jastrow_champ_gl end function qmckl_get_jastrow_champ_gl
end interface end interface
#+end_src #+end_src

120
org/qmckl_mo.org
View File

@ -48,7 +48,7 @@ distance is smaller than a certain radius $r_{\text{cusp}}$. At
distances smaller than $r_{\text{cusp}}$, the MOs are replaced by distances smaller than $r_{\text{cusp}}$, the MOs are replaced by
functions that have the correct electron-nucleus cusp condition and functions that have the correct electron-nucleus cusp condition and
that ensure that the values and the gradients of the MOs are that ensure that the values and the gradients of the MOs are
continuous at $r_{\text{cusp}}$. continuous at $r_{\text{cusp}}$.
A radius $r_{\text{cusp}\, A}$ is given for each nucleus $A$, default is A radius $r_{\text{cusp}\, A}$ is given for each nucleus $A$, default is
zero. If an electron is closer to the nucleus $A$ than zero. If an electron is closer to the nucleus $A$ than
@ -282,7 +282,7 @@ qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context,
"qmckl_set_mo_basis_coefficient", "qmckl_set_mo_basis_coefficient",
"Array too small"); "Array too small");
} }
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; 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); mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double);
double* new_array = (double*) qmckl_malloc(context, mem_info); double* new_array = (double*) qmckl_malloc(context, mem_info);
@ -350,7 +350,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
ctx->mo_basis.coefficient_t = new_array; ctx->mo_basis.coefficient_t = new_array;
qmckl_exit_code rc; qmckl_exit_code rc;
if (ctx->mo_basis.mo_vgl != NULL) { if (ctx->mo_basis.mo_vgl != NULL) {
rc = qmckl_free(context, ctx->mo_basis.mo_vgl); rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
@ -373,11 +373,11 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
To activate the cusp adjustment, the user must enter the radius of To activate the cusp adjustment, the user must enter the radius of
the fitting function for each atom. the fitting function for each atom.
This function requires the computation of the value and gradients This function requires the computation of the value and gradients
of the $s$ AOs at the distance equal to the radius, and the values of the $s$ AOs at the distance equal to the radius, and the values
of the non-$s$ AOs at the center. of the non-$s$ AOs at the center.
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
qmckl_exit_code qmckl_exit_code
qmckl_set_mo_basis_r_cusp(qmckl_context context, qmckl_set_mo_basis_r_cusp(qmckl_context context,
@ -391,19 +391,19 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
qmckl_exit_code rc; qmckl_exit_code rc;
qmckl_context_struct* const ctx = (qmckl_context_struct*) context; qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
if (r_cusp == NULL) { if (r_cusp == NULL) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_2, return qmckl_failwith( context, QMCKL_INVALID_ARG_2,
"qmckl_set_mo_basis_r_cusp", "qmckl_set_mo_basis_r_cusp",
"r_cusp: Null pointer"); "r_cusp: Null pointer");
} }
if (size_max < ctx->nucleus.num) { if (size_max < ctx->nucleus.num) {
return qmckl_failwith( context, QMCKL_INVALID_ARG_3, return qmckl_failwith( context, QMCKL_INVALID_ARG_3,
"qmckl_set_mo_basis_r_cusp", "qmckl_set_mo_basis_r_cusp",
"Array too small"); "Array too small");
} }
// Nullify r_cusp // Nullify r_cusp
if (ctx->mo_basis.r_cusp != NULL) { if (ctx->mo_basis.r_cusp != NULL) {
@ -422,7 +422,7 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = old_point_num * 3 * sizeof(double); mem_info.size = old_point_num * 3 * sizeof(double);
old_coord = (double*) qmckl_malloc(context, mem_info); old_coord = (double*) qmckl_malloc(context, mem_info);
rc = qmckl_get_point(context, 'T', old_coord, (old_point_num * 3)); rc = qmckl_get_point(context, 'T', old_coord, (old_point_num * 3));
if (rc != QMCKL_SUCCESS) { if (rc != QMCKL_SUCCESS) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -450,25 +450,25 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
} }
// Allocate cusp parameters and set them to zero // Allocate cusp parameters and set them to zero
{ {
if (ctx->mo_basis.cusp_param.size[0] != 0) { if (ctx->mo_basis.cusp_param.size[0] != 0) {
rc = qmckl_tensor_free(context, &(ctx->mo_basis.cusp_param)); rc = qmckl_tensor_free(context, &(ctx->mo_basis.cusp_param));
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
} }
int64_t sze[3] = { ctx->mo_basis.mo_num, 4, ctx->nucleus.num }; int64_t sze[3] = { ctx->mo_basis.mo_num, 4, ctx->nucleus.num };
ctx->mo_basis.cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0])); ctx->mo_basis.cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0]));
ctx->mo_basis.cusp_param = qmckl_tensor_set(ctx->mo_basis.cusp_param, 0.); ctx->mo_basis.cusp_param = qmckl_tensor_set(ctx->mo_basis.cusp_param, 0.);
} }
// Evaluate MO value at nucleus without s components // Evaluate MO value at nucleus without s components
qmckl_matrix mo_value_at_nucl_no_s; qmckl_matrix mo_value_at_nucl_no_s;
{ {
mo_value_at_nucl_no_s = qmckl_matrix_alloc(context, ctx->mo_basis.mo_num, ctx->nucleus.num); mo_value_at_nucl_no_s = qmckl_matrix_alloc(context, ctx->mo_basis.mo_num, ctx->nucleus.num);
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3); rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
@ -479,26 +479,26 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
"qmckl_set_mo_basis_r_cusp", "qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at the nuclei"); "Unable to set coordinates at the nuclei");
} }
rc = qmckl_get_mo_basis_mo_value(context, rc = qmckl_get_mo_basis_mo_value(context,
&(qmckl_mat(mo_value_at_nucl_no_s,0,0)), &(qmckl_mat(mo_value_at_nucl_no_s,0,0)),
ctx->mo_basis.mo_num * ctx->nucleus.num); ctx->mo_basis.mo_num * ctx->nucleus.num);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
} }
// Evaluate MO vgl at r_cusp without s components // Evaluate MO vgl at r_cusp without s components
qmckl_tensor mo_vgl_at_r_cusp_s; qmckl_tensor mo_vgl_at_r_cusp_s;
{ {
int64_t sze[3] = { ctx->mo_basis.mo_num, 3, ctx->nucleus.num }; int64_t sze[3] = { ctx->mo_basis.mo_num, 3, ctx->nucleus.num };
mo_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0])); mo_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
} }
{ {
qmckl_tensor ao_vgl_at_r_cusp_s; qmckl_tensor ao_vgl_at_r_cusp_s;
int64_t sze[3] = { ctx->ao_basis.ao_num, 5, ctx->nucleus.num }; int64_t sze[3] = { ctx->ao_basis.ao_num, 5, ctx->nucleus.num };
ao_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0])); ao_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0]));
rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3); rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3);
if (rc != QMCKL_SUCCESS) return rc; if (rc != QMCKL_SUCCESS) return rc;
@ -513,11 +513,11 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
"qmckl_set_mo_basis_r_cusp", "qmckl_set_mo_basis_r_cusp",
"Unable to set coordinates at r_cusp"); "Unable to set coordinates at r_cusp");
} }
rc = qmckl_get_ao_basis_ao_vgl(context, rc = qmckl_get_ao_basis_ao_vgl(context,
&(qmckl_ten3(ao_vgl_at_r_cusp_s,0,0,0)), &(qmckl_ten3(ao_vgl_at_r_cusp_s,0,0,0)),
ctx->ao_basis.ao_num * 5 * ctx->point.num); ctx->ao_basis.ao_num * 5 * ctx->point.num);
for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) { for (int64_t inucl=0 ; inucl<ctx->nucleus.num ; ++inucl) {
for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) { for (int64_t i=0 ; i<ctx->mo_basis.mo_num ; ++i) {
qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = 0.; qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = 0.;
@ -565,11 +565,11 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context,
} }
} }
} }
free(coord); free(coord);
qmckl_matrix_free(context, &mo_value_at_nucl_no_s); qmckl_matrix_free(context, &mo_value_at_nucl_no_s);
qmckl_tensor_free(context, &mo_vgl_at_r_cusp_s); qmckl_tensor_free(context, &mo_vgl_at_r_cusp_s);
// Restore old points // Restore old points
if (old_point_num > 0) { if (old_point_num > 0) {
rc = qmckl_set_point(context, 'T', old_point_num, old_coord, (old_point_num * 3)); rc = qmckl_set_point(context, 'T', old_point_num, old_coord, (old_point_num * 3));
@ -688,7 +688,7 @@ qmckl_get_mo_basis_coefficient (const qmckl_context context,
#+begin_src c :comments org :tangle (eval h_func) #+begin_src c :comments org :tangle (eval h_func)
bool qmckl_mo_basis_provided (const qmckl_context context); bool qmckl_mo_basis_provided (const qmckl_context context);
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) :exports none #+begin_src c :comments org :tangle (eval c) :exports none
bool qmckl_mo_basis_provided(const qmckl_context context) { bool qmckl_mo_basis_provided(const qmckl_context context) {
@ -726,7 +726,7 @@ interface
import import
implicit none implicit none
integer (c_int64_t) , intent(in), value :: context integer (c_int64_t) , intent(in), value :: context
double precision , intent(out) :: coefficient(*) real(c_double) , intent(out) :: coefficient(*)
integer (c_int64_t) , intent(in), value :: size_max integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_get_mo_basis_coefficient end function qmckl_get_mo_basis_coefficient
end interface end interface
@ -738,7 +738,7 @@ interface
import import
implicit none implicit none
integer (c_int64_t) , intent(in), value :: context integer (c_int64_t) , intent(in), value :: context
double precision , intent(in) :: r_cusp(*) real(c_double) , intent(in) :: r_cusp(*)
integer (c_int64_t) , intent(in), value :: size_max integer (c_int64_t) , intent(in), value :: size_max
end function qmckl_set_mo_basis_r_cusp end function qmckl_set_mo_basis_r_cusp
end interface end interface
@ -757,7 +757,7 @@ end interface
size ~mo_num~ can be created. If the integer corresponding to an MO is size ~mo_num~ can be created. If the integer corresponding to an MO is
zero, that MO is dropped and will not be included in the zero, that MO is dropped and will not be included in the
calculation. If the integer is non-zero, the MO will be kept. calculation. If the integer is non-zero, the MO will be kept.
#+begin_src c :comments org :tangle (eval h_func) #+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_exit_code
@ -765,7 +765,7 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
const int32_t* keep, const int32_t* keep,
const int64_t size_max); const int64_t size_max);
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) :exports none #+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code qmckl_exit_code
qmckl_mo_basis_select_mo (const qmckl_context context, qmckl_mo_basis_select_mo (const qmckl_context context,
@ -796,8 +796,8 @@ qmckl_mo_basis_select_mo (const qmckl_context context,
"NULL pointer"); "NULL pointer");
} }
const int64_t mo_num = ctx->mo_basis.mo_num; const int64_t mo_num = ctx->mo_basis.mo_num;
const int64_t ao_num = ctx->ao_basis.ao_num; const int64_t ao_num = ctx->ao_basis.ao_num;
if (size_max < mo_num) { if (size_max < mo_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
@ -906,7 +906,7 @@ qmckl_get_mo_basis_mo_value(qmckl_context context,
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_value(*) real(c_double), intent(out) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value end function qmckl_get_mo_basis_mo_value
end interface end interface
@ -972,7 +972,7 @@ qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_value(*) real(c_double), intent(out) :: mo_value(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_value_inplace end function qmckl_get_mo_basis_mo_value_inplace
end interface end interface
@ -1049,12 +1049,12 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
#+end_src #+end_src
#+begin_src c :comments org :tangle (eval c) :noweb yes :exports none #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
if (ctx->mo_basis.mo_vgl_date == ctx->point.date) { if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
// mo_vgl has been computed at this step: Just copy the data. // mo_vgl has been computed at this step: Just copy the data.
double * v = &(ctx->mo_basis.mo_value[0]); double * v = &(ctx->mo_basis.mo_value[0]);
double * vgl = &(ctx->mo_basis.mo_vgl[0]); double * vgl = &(ctx->mo_basis.mo_vgl[0]);
for (int i=0 ; i<ctx->point.num ; ++i) { for (int i=0 ; i<ctx->point.num ; ++i) {
@ -1099,7 +1099,7 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
ctx->ao_basis.ao_num, ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num, ctx->mo_basis.mo_num,
ctx->point.num, ctx->point.num,
ctx->ao_basis.ao_nucl, ctx->ao_basis.ao_nucl,
ctx->ao_basis.ao_ang_mom, ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance, ctx->electron.en_distance,
ctx->mo_basis.r_cusp, ctx->mo_basis.r_cusp,
@ -1177,7 +1177,7 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j) mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
end do end do
end do end do
end function qmckl_compute_mo_basis_mo_value_doc_f end function qmckl_compute_mo_basis_mo_value_doc_f
#+end_src #+end_src
@ -1394,7 +1394,7 @@ qmckl_get_mo_basis_mo_vgl(qmckl_context context,
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_vgl(*) real(c_double), intent(out) :: mo_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_vgl end function qmckl_get_mo_basis_mo_vgl
end interface end interface
@ -1460,7 +1460,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
double precision, intent(out) :: mo_vgl(*) real(c_double), intent(out) :: mo_vgl(*)
integer (c_int64_t) , intent(in) , value :: size_max integer (c_int64_t) , intent(in) , value :: size_max
end function qmckl_get_mo_basis_mo_vgl_inplace end function qmckl_get_mo_basis_mo_vgl_inplace
end interface end interface
@ -1568,7 +1568,7 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
ctx->ao_basis.ao_num, ctx->ao_basis.ao_num,
ctx->mo_basis.mo_num, ctx->mo_basis.mo_num,
ctx->point.num, ctx->point.num,
ctx->ao_basis.ao_nucl, ctx->ao_basis.ao_nucl,
ctx->ao_basis.ao_ang_mom, ctx->ao_basis.ao_ang_mom,
ctx->electron.en_distance, ctx->electron.en_distance,
ctx->nucleus.coord, ctx->nucleus.coord,
@ -1941,7 +1941,7 @@ integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, &
if ( (en_distance(inucl,i) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle if ( (en_distance(inucl,i) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
mo_value(:,i) = mo_value(:,i) + coefficient_t(:,k) * ao_value(k,i) mo_value(:,i) = mo_value(:,i) + coefficient_t(:,k) * ao_value(k,i)
end do ! k end do ! k
do inucl=1,nucl_num do inucl=1,nucl_num
r = en_distance(inucl,i) r = en_distance(inucl,i)
if (r > r_cusp(inucl)) cycle if (r > r_cusp(inucl)) cycle
@ -1974,7 +1974,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
const qmckl_tensor cusp_param, const qmckl_tensor cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_value, const double* ao_value,
double* const mo_value ); double* const mo_value );
#+end_src #+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc")) #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
@ -1994,7 +1994,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f
const double* cusp_param, const double* cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_value, const double* ao_value,
double* const mo_value ); double* const mo_value );
#+end_src #+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc")) #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp_doc"))
@ -2070,17 +2070,17 @@ qmckl_compute_mo_basis_mo_value_cusp (const qmckl_context context,
double* const mo_value ) double* const mo_value )
{ {
qmckl_exit_code rc; qmckl_exit_code rc;
#ifdef HAVE_HPC #ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, rc = qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp, ao_nucl, ao_ang_mom, en_distance, r_cusp,
cusp_param_tensor, coefficient_t, ao_value, mo_value ); cusp_param_tensor, coefficient_t, ao_value, mo_value );
#else #else
double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor); double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor);
rc = qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, rc = qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, r_cusp, ao_nucl, ao_ang_mom, en_distance, r_cusp,
cusp_param, coefficient_t, ao_value, mo_value ); cusp_param, coefficient_t, ao_value, mo_value );
qmckl_free(context, cusp_param); qmckl_free(context, cusp_param);
#endif #endif
@ -2106,7 +2106,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
const qmckl_tensor cusp_param, const qmckl_tensor cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_value, const double* ao_value,
double* const mo_value ); double* const mo_value );
#endif #endif
#+end_src #+end_src
@ -2179,7 +2179,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
for (int64_t m=n ; m < nidx ; m+=1) { for (int64_t m=n ; m < nidx ; m+=1) {
const double* restrict ck = coefficient_t + idx[m]*mo_num; const double* restrict ck = coefficient_t + idx[m]*mo_num;
const double a1 = av1[m]; const double a1 = av1[m];
#ifdef HAVE_OPENMP #ifdef HAVE_OPENMP
#pragma omp simd #pragma omp simd
#endif #endif
@ -2283,7 +2283,7 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5 mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
end do end do
end do end do
! Cusp adjustment ! Cusp adjustment
do inucl=1,nucl_num do inucl=1,nucl_num
r = en_distance(inucl,j) r = en_distance(inucl,j)
@ -2299,7 +2299,7 @@ integer function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f(context, &
cusp_param(i,3,inucl) + r* cusp_param(i,4,inucl) )) cusp_param(i,3,inucl) + r* cusp_param(i,4,inucl) ))
c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) + & c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) + &
r * 3.d0 * cusp_param(i,4,inucl) r * 3.d0 * cusp_param(i,4,inucl)
mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1 mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1
mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1 mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1
@ -2319,7 +2319,7 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
#+end_src #+end_src
# #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp")) # #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp"))
#+RESULTS: #+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org #+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp ( qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp (
@ -2337,7 +2337,7 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
const qmckl_tensor cusp_param, const qmckl_tensor cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_vgl, const double* ao_vgl,
double* const mo_vgl ); double* const mo_vgl );
#+end_src #+end_src
#+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc")) #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
@ -2359,7 +2359,7 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f
const double* cusp_param, const double* cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_vgl, const double* ao_vgl,
double* const mo_vgl ); double* const mo_vgl );
#+end_src #+end_src
#+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc")) #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp_doc"))
@ -2443,7 +2443,7 @@ qmckl_compute_mo_basis_mo_vgl_cusp (const qmckl_context context,
double* const mo_vgl ) double* const mo_vgl )
{ {
qmckl_exit_code rc; qmckl_exit_code rc;
#ifdef HAVE_HPC #ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, rc = qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
ao_nucl, ao_ang_mom, en_distance, nucl_coord_matrix, ao_nucl, ao_ang_mom, en_distance, nucl_coord_matrix,
@ -2487,7 +2487,7 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
const qmckl_tensor cusp_param, const qmckl_tensor cusp_param,
const double* coefficient_t, const double* coefficient_t,
const double* ao_vgl, const double* ao_vgl,
double* const mo_vgl ); double* const mo_vgl );
#endif #endif
#+end_src #+end_src
@ -2670,7 +2670,7 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
may get out of the range of double precision. A simple fix is to may get out of the range of double precision. A simple fix is to
rescale the MO coefficients to put back the determinants in the rescale the MO coefficients to put back the determinants in the
correct range. correct range.
#+begin_src c :comments org :tangle (eval h_func) :noweb yes #+begin_src c :comments org :tangle (eval h_func) :noweb yes
qmckl_exit_code qmckl_exit_code
qmckl_mo_basis_rescale(qmckl_context context, qmckl_mo_basis_rescale(qmckl_context context,
@ -2817,7 +2817,7 @@ const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
const int64_t point_num = walk_num*elec_num; const int64_t point_num = walk_num*elec_num;
rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num);
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
@ -2945,7 +2945,7 @@ for (int i=0 ; i< point_num; ++i) {
rc = qmckl_context_touch(context); rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num); rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), point_num*chbrclf_mo_num);
assert (rc == QMCKL_SUCCESS); assert (rc == QMCKL_SUCCESS);
@ -3054,7 +3054,7 @@ printf("\n");
assert(mo_coefficient_new[0][i] == mo_coefficient[i + ao_num*2]); assert(mo_coefficient_new[0][i] == mo_coefficient[i + ao_num*2]);
assert(mo_coefficient_new[1][i] == mo_coefficient[i + ao_num*5]); assert(mo_coefficient_new[1][i] == mo_coefficient[i + ao_num*5]);
} }
} }

2
tools/lib.org
View File

@ -46,7 +46,7 @@ f_of_c_d = { '' : ''
, 'uint64_t' : 'integer (c_int64_t)' , 'uint64_t' : 'integer (c_int64_t)'
, 'float' : 'real (c_float )' , 'float' : 'real (c_float )'
, 'double' : 'real (c_double )' , 'double' : 'real (c_double )'
, 'char' : 'character' , 'char' : 'character(c_char )'
} }
#+END_SRC #+END_SRC