diff --git a/org/examples.org b/org/examples.org index eef87b4..8650c40 100644 --- a/org/examples.org +++ b/org/examples.org @@ -7,16 +7,20 @@ 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: \[ 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) - \phi_j(\mathbf{r}_k) \delta \mathbf{r} + \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 \] @@ -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