mirror of
https://github.com/TREX-CoE/trexio.git
synced 2024-08-24 22:21:43 +02:00
Merge branch 'master' into cppcheck
This commit is contained in:
commit
bb90e6b159
@ -200,7 +200,7 @@ program print_energy
|
||||
#+begin_src f90
|
||||
call getarg(1, filename)
|
||||
|
||||
f = trexio_open (filename, 'r', TREXIO_HDF5, rc)
|
||||
f = trexio_open (filename, 'r', TREXIO_AUTO, rc)
|
||||
if (rc /= TREXIO_SUCCESS) then
|
||||
call trexio_string_of_error(rc, err_msg)
|
||||
print *, 'Error opening TREXIO file: '//trim(err_msg)
|
||||
@ -410,6 +410,139 @@ program print_energy
|
||||
end program
|
||||
#+end_src
|
||||
|
||||
** Python
|
||||
:PROPERTIES:
|
||||
:header-args: :tangle print_energy.py
|
||||
:END:
|
||||
|
||||
#+begin_src python
|
||||
import sys
|
||||
import trexio
|
||||
import numpy as np
|
||||
|
||||
BUFSIZE = 100000
|
||||
#+end_src
|
||||
|
||||
This program computes the energy as:
|
||||
|
||||
\[
|
||||
E = E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle j | h | i \rangle\,
|
||||
+\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j
|
||||
\rangle\; \textrm{ with } \; 0 < i,j,k,l \le n
|
||||
\]
|
||||
One needs to read from the TREXIO file:
|
||||
|
||||
- $n$ :: The number of molecular orbitals
|
||||
- $E_{\text{NN}}$ :: The nuclear repulsion energy
|
||||
- $\gamma_{ij}$ :: The one-body reduced density matrix
|
||||
- $\langle j |h| i \rangle$ :: The one-electron Hamiltonian integrals
|
||||
- $\Gamma_{ijkl}$ :: The two-body reduced density matrix
|
||||
- $\langle k l | i j \rangle$ :: The electron repulsion integrals
|
||||
|
||||
*** Obtain the name of the TREXIO file from the command line, and open it for reading
|
||||
|
||||
#+begin_src python
|
||||
filename = sys.argv[1]
|
||||
f = trexio.File(filename, 'r', trexio.TREXIO_AUTO)
|
||||
#+end_src
|
||||
|
||||
*** Read the nuclear repulsion energy
|
||||
|
||||
#+begin_src python
|
||||
E_nn = trexio.read_nucleus_repulsion(f)
|
||||
#+end_src
|
||||
|
||||
*** Read the number of molecular orbitals
|
||||
|
||||
#+begin_src python
|
||||
n = trexio.read_mo_num(f)
|
||||
#+end_src
|
||||
|
||||
*** Read one-electron quantities
|
||||
|
||||
#+begin_src python
|
||||
if not trexio.has_mo_1e_int_core_hamiltonian(f):
|
||||
print("No core hamiltonian in file")
|
||||
sys.exit(-1)
|
||||
|
||||
h0 = trexio.read_mo_1e_int_core_hamiltonian(f)
|
||||
|
||||
if not trexio.has_rdm_1e(f):
|
||||
print("No 1e RDM in file")
|
||||
sys.exit(-1)
|
||||
|
||||
D = trexio.read_rdm_1e(f)
|
||||
#+end_src
|
||||
|
||||
*** Read two-electron quantities
|
||||
|
||||
**** Electron repulsion integrals
|
||||
|
||||
#+begin_src python
|
||||
if not trexio.has_mo_2e_int_eri(f):
|
||||
print("No electron repulsion integrals in file")
|
||||
sys.exit(-1)
|
||||
|
||||
size_max = trexio.read_mo_2e_int_eri_size(f)
|
||||
|
||||
offset = 0
|
||||
icount = BUFSIZE
|
||||
feof = False
|
||||
W = np.zeros( (n,n,n,n) )
|
||||
while not feof:
|
||||
buffer_index, buffer_values, icount, feof = trexio.read_mo_2e_int_eri(f, offset, icount)
|
||||
for m in range(icount):
|
||||
i, j, k, l = buffer_index[m]
|
||||
W[i,j,k,l] = buffer_values[m]
|
||||
W[k,j,i,l] = buffer_values[m]
|
||||
W[i,l,k,j] = buffer_values[m]
|
||||
W[k,l,i,j] = buffer_values[m]
|
||||
W[j,i,l,k] = buffer_values[m]
|
||||
W[j,k,l,i] = buffer_values[m]
|
||||
W[l,i,j,k] = buffer_values[m]
|
||||
W[l,k,j,i] = buffer_values[m]
|
||||
#+end_src
|
||||
|
||||
**** Reduced density matrix
|
||||
|
||||
#+begin_src python
|
||||
if not trexio.has_rdm_2e(f):
|
||||
print("No two-body density matrix in file")
|
||||
|
||||
offset = 0
|
||||
icount = BUFSIZE
|
||||
feof = False
|
||||
G = np.zeros( (n,n,n,n) )
|
||||
while not feof:
|
||||
buffer_index, buffer_values, icount, feof = trexio.read_rdm_2e(f, offset, icount)
|
||||
for m in range(icount):
|
||||
i, j, k, l = buffer_index[m]
|
||||
G[i,j,k,l] = buffer_values[m]
|
||||
|
||||
#+end_src
|
||||
|
||||
*** Compute the energy
|
||||
|
||||
When the orbitals are real, we can use
|
||||
\begin{eqnarray*}
|
||||
E &=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle j | h | i \rangle\,
|
||||
+\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j
|
||||
\rangle \\
|
||||
&=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle i | h | j \rangle\,
|
||||
+\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle i j | k l
|
||||
\rangle \\
|
||||
\end{eqnarray*}
|
||||
|
||||
#+begin_src python
|
||||
G = np.reshape(G, (n*n, n*n) )
|
||||
W = np.reshape(W, (n*n, n*n) )
|
||||
E = E_nn
|
||||
E += 0.5*sum( [ np.dot(G[:,l], W[:,l]) for l in range(n*n) ] )
|
||||
E += sum( [ np.dot(D[:,l], h0[:,l]) for l in range(n) ] )
|
||||
|
||||
print (f"Energy: {E}")
|
||||
#+end_src
|
||||
|
||||
* Reading determinants
|
||||
|
||||
** Fortran
|
||||
|
Loading…
Reference in New Issue
Block a user