diff --git a/configure.ac b/configure.ac index a434d7a..6f8774e 100644 --- a/configure.ac +++ b/configure.ac @@ -35,7 +35,7 @@ AC_PREREQ([2.69]) -AC_INIT([qmckl],[0.3.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html]) +AC_INIT([qmckl],[0.4.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html]) AC_CONFIG_AUX_DIR([tools]) AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11]) diff --git a/org/qmckl_examples.org b/org/qmckl_examples.org index b01cbac..0408d16 100644 --- a/org/qmckl_examples.org +++ b/org/qmckl_examples.org @@ -82,6 +82,223 @@ 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) ]) ) +linspace = [ None for i in range(3) ] +step = [ None for i in range(3) ] +for a in range(3): + linspace[a], step[a] = np.linspace(rmin[a]-shift[a], + rmax[a]+shift[a], + num=nx[a], + retstep=True) + +dr = step[0] * step[1] * step[2] + #+end_src + + #+RESULTS: + + 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 :exports code +point = [] +for x in linspace[0]: + for y in linspace[1]: + for z in linspace[2]: + point += [ [x, y, z] ] + +point = np.array(point) +point_num = len(point) +qmckl.set_point(context, 'N', point_num, np.reshape(point, (point_num*3))) + #+end_src + + #+RESULTS: + : None + + 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 :exports code +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 :exports code +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 + +* C +** Check numerically that MOs are orthonormal, with cusp fitting + + 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} + \] + \[ + 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 + \] + + We apply the cusp fitting procedure, so the MOs might deviate + slightly from orthonormality. + + #+begin_src c :exports code +#include +#include + +int main(int argc, char** argv) +{ + const char* trexio_filename = "..//share/qmckl/test_data/h2o_5z.h5"; + qmckl_exit_code rc = QMCKL_SUCCESS; + #+end_src + + 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 + determinant for the water molecule in the cc-pV5Z basis set. + + #+begin_src c :exports code + qmckl_context context = qmckl_context_create(); + if (context == NULL) { + fprintf(stderr, "Error creating context\n"); + exit(1); + } + + rc = qmckl_trexio_read(context, trexio_filename, strlen(trexio_filename)); + + if (rc != QMCKL_SUCCESS) { + fprintf(stderr, "Error reading TREXIO file:\n%s\n", qmckl_string_of_error(rc)); + exit(1); + } + + #+end_src + + We impose the electron-nucleus cusp fitting to occur when the + electrons are close to the nuclei. The critical distance + is 0.5 atomic units for hydrogens and 0.1 for the oxygen. + To identify which atom is an oxygen and which are hydrogens, we + fetch the nuclear charges from the context. + + #+begin_src python :exports code + int64_t nucl_num; + + rc = qmckl_get_nucleus_num(context, &nucl_num); + + if (rc != QMCKL_SUCCESS) { + fprintf(stderr, "Error getting nucl_num:\n%s\n", qmckl_string_of_error(rc)); + exit(1); + } + + + double nucl_charge[nucl_num]; + + rc = qmckl_get_nucleus_charge(context, &(nucl_charge[0])); + + if (rc != QMCKL_SUCCESS) { + fprintf(stderr, "Error getting nucl_charge:\n%s\n", qmckl_string_of_error(rc)); + exit(1); + } + + + double r_c[nucl_num]; + + for (size_t i=0 ; i