Added Python example.

org/examples.org

@@ -7,10 +7,9 @@ For simplicity, we assume that the wave function parameters are stored
 in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
 
 * Python
-
 ** Check numerically that MOs are orthonormal
 
-   In this example, we will compute the numerically the overlap
+   In this example, we will compute numerically the overlap
    between the molecular orbitals:
    
    \[
@@ -18,6 +17,11 @@ in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
    \text{d}\mathbf{r} \sim \sum_{k=1}^{N} \phi_i(\mathbf{r}_k)
    \phi_j(\mathbf{r}_k) \delta \mathbf{r} 
    \]
+   \[
+      S_{ij} = \langle \phi_i | \phi_j \rangle 
+   \sim \sum_{k=1}^{N} \langle \phi_i | \mathbf{r}_k \rangle
+   \langle \mathbf{r}_k | \phi_j \rangle
+   \]
    
 
    #+begin_src python :session 
@@ -27,9 +31,9 @@ import qmckl
 
    #+RESULTS:
    
-
    First, we create a context for the QMCkl calculation, and load the
-   wave function stored in =h2o_5z.h5= inside it:
+   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.
    
    #+begin_src python :session 
 trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5"
@@ -41,12 +45,12 @@ qmckl.trexio_read(context, trexio_filename)
    #+RESULTS:
    : None
 
-   We now define the grid points as a regular grid around the
+   We now define the grid points $\mathbf{r}_k$ as a regular grid around the
    molecule.
 
    We fetch the nuclear coordinates from the context,
    
-   #+begin_src python :session :results output
+   #+begin_src python :session :results output :export both
 nucl_num = qmckl.get_nucleus_num(context)
 
 nucl_charge = qmckl.get_nucleus_charge(context, nucl_num)
@@ -61,21 +65,22 @@ for i in range(nucl_num):
                              nucl_coord[i,2]) )
    #+end_src
 
-   #+RESULTS:
-   : 8  +0.000000 +0.000000 +0.000000
-   : 1  -1.430429 +0.000000 -1.107157
-   : 1  +1.430429 +0.000000 -1.107157
+   #+begin_example
+8  +0.000000 +0.000000 +0.000000
+1  -1.430429 +0.000000 -1.107157
+1  +1.430429 +0.000000 -1.107157
+   #+end_example
 
    and compute the coordinates of the grid points:
    
    #+begin_src python :session
-nx = ( 40, 40, 40 )
+nx = ( 120, 120, 120 )
+shift = np.array([5.,5.,5.])
 point_num = nx[0] * nx[1] * nx[2]
 
 rmin  = np.array( list([ np.min(nucl_coord[:,a]) for a in range(3) ]) )
 rmax  = np.array( list([ np.max(nucl_coord[:,a]) for a in range(3) ]) )
 
-shift = np.array([5.,5.,5.])
 
 linspace = [ None for i in range(3) ]
 step     = [ None for i in range(3) ]
@@ -86,15 +91,13 @@ for a in range(3):
                                        retstep=True)
 
 dr = step[0] * step[1] * step[2]
-dr
    #+end_src
 
    #+RESULTS:
-   : 0.024081249137090373
    
-   Now the grid is ready, we can create the list of grid points on
-   which the MOs will be evaluated, and transfer them to the QMCkl
-   context:
+   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
+   transfer them to the QMCkl context:
 
    #+begin_src python :session
 point = []
@@ -104,13 +107,63 @@ for x in linspace[0]:
             point += [ [x, y, z] ]
 
 point = np.array(point)
-qmckl.set_point(context, 'N', len(point), point)
+point_num = len(point)
+qmckl.set_point(context, 'N', point_num, np.reshape(point, (point_num*3)))
    #+end_src
 
    #+RESULTS:
+   : None
    
-   Then, will first evaluate all the MOs at the grid points, and then we will
-   compute the overlap between all the MOs.
+   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 =
+   \phi_i(\mathbf{r}_k)$.
+
+   #+begin_src python :session :results output :export both
+import time
+
+mo_num = qmckl.get_mo_basis_mo_num(context)
+
+before   = time.time()
+mo_value = qmckl.get_mo_basis_mo_value(context, point_num*mo_num)
+after    = time.time()
+
+mo_value = np.reshape( mo_value, (point_num, mo_num) )
+
+print("Number of MOs: ", mo_num)
+print("Number of grid points: ", point_num)
+print("Execution time : ", (after - before), "seconds")
+
+   #+end_src
+
+   #+begin_example
+Number of MOs:  201
+Number of grid points:  1728000
+Execution time :  3.511528968811035 seconds
+   #+end_example
+
+   and finally we compute the overlap between all the MOs as
+   $M^\dagger M$.
+
+   #+begin_src python :session :results output
+overlap = mo_value.T @ mo_value * dr
+print (overlap)
+   #+end_src
+
+   #+begin_example
+   [[ 9.88693941e-01  2.34719693e-03 -1.50518232e-08 ...  3.12084178e-09
+     -5.81064929e-10  3.70130091e-02]
+    [ 2.34719693e-03  9.99509628e-01  3.18930040e-09 ... -2.46888958e-10
+     -1.06064273e-09 -7.65567973e-03]
+    [-1.50518232e-08  3.18930040e-09  9.99995073e-01 ... -5.84882580e-06
+     -1.21598117e-06  4.59036468e-08]
+    ...
+    [ 3.12084178e-09 -2.46888958e-10 -5.84882580e-06 ...  1.00019107e+00
+     -2.03342837e-04 -1.36954855e-08]
+    [-5.81064929e-10 -1.06064273e-09 -1.21598117e-06 ... -2.03342837e-04
+      9.99262427e-01  1.18264754e-09]
+    [ 3.70130091e-02 -7.65567973e-03  4.59036468e-08 ... -1.36954855e-08
+      1.18264754e-09  8.97215950e-01]]
+   #+end_example
 
 * Fortran
 ** Checking errors
@@ -278,7 +331,7 @@ program ao_grid
    We give the points to QMCkl:
 
    #+begin_src f90 
-  rc = qmckl_set_point(qmckl_ctx, 'T', points, point_num)
+  rc = qmckl_set_point(qmckl_ctx, 'T', point_num, points, size(points)*1_8 )
   call qmckl_check_error(rc, 'Setting points')
    #+end_src
 
@@ -292,11 +345,11 @@ program ao_grid
   call qmckl_check_error(rc, 'Setting points')
    #+end_src
 
-   We finally print the value of the AO:
+   We finally print the value and Laplacian 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)
+     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_src