mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-17 02:41:43 +02:00
Added Python example.
This commit is contained in:
parent
ce1aeb324d
commit
7e854175cc
@ -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
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user