From 26fe759209c821a90eb3f1e4fdc86a75542445a5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama <scemama@irsamc.ups-tlse.fr> Date: Sun, 27 Feb 2022 12:35:58 +0100 Subject: [PATCH] Added examples.org --- org/examples.org | 190 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 190 insertions(+) create mode 100644 org/examples.org diff --git a/org/examples.org b/org/examples.org new file mode 100644 index 0000000..6f8f4d0 --- /dev/null +++ b/org/examples.org @@ -0,0 +1,190 @@ +#+TITLE: Code examples +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +In this section, we present examples of usage of QMCkl. +For simplicity, we assume that the wave function parameters are stores +in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file. + +* Checking errors + + All QMCkl functions return an error code. A convenient way to handle + errors is to write an error-checking function that displays the + error in text format and exits the program. + + #+NAME: qmckl_check_error + #+begin_src f90 +subroutine qmckl_check_error(rc, message) + use qmckl + implicit none + integer(qmckl_exit_code), intent(in) :: rc + character(len=*) , intent(in) :: message + character(len=128) :: str_buffer + if (rc /= QMCKL_SUCCESS) then + print *, message + call qmckl_string_of_error(rc, str_buffer) + print *, str_buffer + call exit(rc) + end if +end subroutine qmckl_check_error + #+end_src + +* Computing an atomic orbital on a grid + :PROPERTIES: + :header-args: :tangle ao_grid.f90 + :END: + + The following program, in Fortran, computes the values of an atomic + orbital on a regular 3-dimensional grid. The 100^3 grid points are + automatically defined, such that the molecule fits in a box with 5 + atomic units in the borders. + + This program uses the ~qmckl_check_error~ function defined above. + + To use this program, run + + #+begin_src bash :tangle no +$ ao_grid <trexio_file> <AO_id> <point_num> + #+end_src + + + #+begin_src f90 :noweb yes +<<qmckl_check_error>> + +program ao_grid + use qmckl + implicit none + + integer(qmckl_context) :: qmckl_ctx ! QMCkl context + integer(qmckl_exit_code) :: rc ! Exit code of QMCkl functions + + character(len=128) :: trexio_filename + character(len=128) :: str_buffer + integer :: ao_id + integer :: point_num_x + + integer(c_int64_t) :: nucl_num + double precision, allocatable :: nucl_coord(:,:) + + integer(c_int64_t) :: point_num + integer(c_int64_t) :: ao_num + integer(c_int64_t) :: ipoint, i, j, k + double precision :: x, y, z, dr(3) + double precision :: rmin(3), rmax(3) + double precision, allocatable :: points(:,:) + double precision, allocatable :: ao_vgl(:,:,:) + #+end_src + + Start by fetching the command-line arguments: + + #+begin_src f90 + if (iargc() /= 3) then + print *, 'Syntax: ao_grid <trexio_file> <AO_id> <point_num>' + call exit(-1) + end if + call getarg(1, trexio_filename) + call getarg(2, str_buffer) + read(str_buffer, *) ao_id + call getarg(3, str_buffer) + read(str_buffer, *) point_num_x + + if (point_num_x < 0 .or. point_num_x > 300) then + print *, 'Error: 0 < point_num < 300' + call exit(-1) + end if + #+end_src + + Create the QMCkl context and initialize it with the wave function + present in the TREXIO file: + + #+begin_src f90 + qmckl_ctx = qmckl_context_create() + rc = qmckl_trexio_read(qmckl_ctx, trexio_filename, 1_8*len(trim(trexio_filename))) + call qmckl_check_error(rc, 'Read TREXIO') + #+end_src + + 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. + + #+begin_src f90 + rc = qmckl_get_nucleus_num(qmckl_ctx, nucl_num) + call qmckl_check_error(rc, 'Get nucleus num') + + allocate( nucl_coord(3, nucl_num) ) + rc = qmckl_get_nucleus_coord(qmckl_ctx, 'N', nucl_coord, 3_8*nucl_num) + call qmckl_check_error(rc, 'Get nucleus coord') + #+end_src + + We now compute the coordinates of opposite points of the box, and + the distance between points along the 3 directions: + + #+begin_src f90 + rmin(1) = minval( nucl_coord(1,:) ) - 5.d0 + rmin(2) = minval( nucl_coord(2,:) ) - 5.d0 + rmin(3) = minval( nucl_coord(3,:) ) - 5.d0 + + rmax(1) = maxval( nucl_coord(1,:) ) + 5.d0 + rmax(2) = maxval( nucl_coord(2,:) ) + 5.d0 + rmax(3) = maxval( nucl_coord(3,:) ) + 5.d0 + + dr(1:3) = (rmax(1:3) - rmin(1:3)) / dble(point_num_x-1) + #+end_src + + We now produce the list of point coordinates where the AO will be + evaluated: + + #+begin_src f90 + point_num = point_num_x**3 + allocate( points(point_num, 3) ) + ipoint=0 + z = rmin(3) + do k=1,point_num_x + y = rmin(2) + do j=1,point_num_x + x = rmin(1) + do i=1,point_num_x + ipoint = ipoint+1 + points(ipoint,1) = x + points(ipoint,2) = y + points(ipoint,3) = z + x = x + dr(1) + end do + y = y + dr(2) + end do + z = z + dr(3) + end do + #+end_src + + We give the points to QMCkl: + + #+begin_src f90 + rc = qmckl_set_point(qmckl_ctx, 'T', points, point_num) + call qmckl_check_error(rc, 'Setting points') + #+end_src + + We allocate the space required to retrieve the values, gradients and + Laplacian of all AOs, and ask to retrieve the values of the + AOs computed at the point positions. For that, we first need to know + the number of AOs: + + #+begin_src f90 + rc = qmckl_get_ao_basis_ao_num(qmckl_ctx, ao_num) + call qmckl_check_error(rc, 'Getting ao_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) + call qmckl_check_error(rc, 'Setting points') + #+end_src + + We finally print the value of the AO: + + #+begin_src f90 + do ipoint=1, point_num + print '(3(F16.10,X),E20.10)', points(ipoint, 1:3), ao_vgl(ao_id,1,ipoint) + end do + #+end_src + + #+begin_src f90 + deallocate( nucl_coord, points, ao_vgl ) +end program ao_grid + #+end_src