Contents

More advanced usage of the VeloxChem API

More advanced usage of the VeloxChem API

Now that we know how to create a molecule, generate a basis set, and run an SCF calculation through the API, we can venture into more complex workflows. For example, computing integrals and visualizing molecular orbitals. We will keep using our water molecule as an example system.

import veloxchem as vlx

h2o_xyz = """3
water                                                                                                                          
O    0.000000000000        0.000000000000        0.000000000000                         
H    0.000000000000        0.740848095288        0.582094932012                         
H    0.000000000000       -0.740848095288        0.582094932012
"""

h2o = vlx.Molecule.from_xyz_string(h2o_xyz)
print(h2o.get_string())

basis = vlx.MolecularBasis.read(h2o, "aug-cc-pvdz")
print(basis.get_string(h2o))

* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.
Molecular Geometry (Angstroms)
================================

  Atom         Coordinate X          Coordinate Y          Coordinate Z  

  O           0.000000000000        0.000000000000        0.000000000000
  H           0.000000000000        0.740848095288        0.582094932012
  H           0.000000000000       -0.740848095288        0.582094932012


Molecular Basis (Atomic Basis)
================================

Basis: AUG-CC-PVDZ                                    

  Atom Contracted GTOs          Primitive GTOs           

  H   (3S,2P)                  (5S,2P)                  
  O   (4S,3P,2D)               (18S,5P,2D)              

Contracted Basis Functions : 41                       
Primitive Basis Functions  : 65

Visualizing molecular orbitals and densities

As we’ve seen in the previous episode, molecular orbitals and AO-basis densities are saved in the scf_tensors dictionary attached to an SCF driver object. This is available whenever the SCF calculation completes.

scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(h2o, basis)

                                                                                                                          
                                            Self Consistent Field Driver Setup                                            
                                           ====================================                                           
                                                                                                                          
                   Wave Function Model             : Spin-Restricted Hartree-Fock                                         
                   Initial Guess Model             : Superposition of Atomic Densities                                    
                   Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace                     
                   Max. Number of Iterations       : 50                                                                   
                   Max. Number of Error Vectors    : 10                                                                   
                   Convergence Threshold           : 1.0e-06                                                              
                   ERI Screening Scheme            : Cauchy Schwarz + Density                                             
                   ERI Screening Mode              : Dynamic                                                              
                   ERI Screening Threshold         : 1.0e-12                                                              
                   Linear Dependence Threshold     : 1.0e-06                                                              
                                                                                                                          
* Info * Nuclear repulsion energy: 9.3436381580 au                                                                        
                                                                                                                          
* Info * Overlap matrix computed in 0.00 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.05 sec.                                                                   
                                                                                                                          
* Info * SAD initial guess computed in 0.00 sec.                                                                          
                                                                                                                          
* Info * Starting Reduced Basis SCF calculation...                                                                        
* Info * ...done. SCF energy in reduced basis set: -75.990361787820 au. Time: 1.26 sec.                                   
                                                                                                                          
* Info * Overlap matrix computed in 0.08 sec.                                                                             
                                                                                                                          
* Info * Kinetic energy matrix computed in 0.00 sec.                                                                      
                                                                                                                          
* Info * Nuclear potential matrix computed in 0.00 sec.                                                                   
                                                                                                                          
* Info * Orthogonalization matrix computed in 0.00 sec.                                                                   
                                                                                                                          
                                                                                                                          
               Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change               
               --------------------------------------------------------------------------------------------               
                  1       -76.040225772729    0.0000000000      0.10527360      0.01410687      0.00000000                
                  2       -76.041581009463   -0.0013552367      0.02396700      0.00364108      0.03764551                
                  3       -76.041683167061   -0.0001021576      0.00872358      0.00122472      0.00801861                
                  4       -76.041694355980   -0.0000111889      0.00512350      0.00094645      0.00386334                
                  5       -76.041697479316   -0.0000031233      0.00046249      0.00008300      0.00125597                
                  6       -76.041697547982   -0.0000000687      0.00007707      0.00001089      0.00028877                
                  7       -76.041697549783   -0.0000000018      0.00000995      0.00000158      0.00004453                
                  8       -76.041697549817   -0.0000000000      0.00000212      0.00000032      0.00000530                
                  9       -76.041697549818   -0.0000000000      0.00000039      0.00000005      0.00000148                
                                                                                                                          
               *** SCF converged in 9 iterations. Time: 7.19 sec.                                                         
                                                                                                                          
               Spin-Restricted Hartree-Fock:                                                                              
               -----------------------------                                                                              
               Total Energy                       :      -76.0416975498 au                                                
               Electronic Energy                  :      -85.3853357078 au                                                
               Nuclear Repulsion Energy           :        9.3436381580 au                                                
               ------------------------------------                                                                       
               Gradient Norm                      :        0.0000003901 au                                                
                                                                                                                          
                                                                                                                          
               Ground State Information                                                                                   
               ------------------------                                                                                   
               Charge of Molecule            :  0.0                                                                       
               Multiplicity (2S+1)           :  1.0                                                                       
               Magnetic Quantum Number (M_S) :  0.0                                                                       
                                                                                                                          
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   1:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.0 Energy:  -20.57504 au                                                                      
               (   1 O   1s  :    -1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.0 Energy:   -1.36523 au                                                                      
               (   1 O   2s  :     0.45) (   1 O   3s  :     0.36) (   2 H   1s  :     0.21)                              
               (   3 H   1s  :     0.21)                                                                                  
                                                                                                                          
               Molecular Orbital No.   3:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.0 Energy:   -0.72596 au                                                                      
               (   1 O   1p-1:    -0.50) (   1 O   2p-1:    -0.19) (   2 H   1s  :    -0.35)                              
               (   3 H   1s  :     0.35)                                                                                  
                                                                                                                          
               Molecular Orbital No.   4:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.0 Energy:   -0.59018 au                                                                      
               (   1 O   2s  :    -0.15) (   1 O   3s  :    -0.33) (   1 O   1p0 :     0.55)                              
               (   1 O   2p0 :     0.32) (   2 H   1s  :     0.21) (   3 H   1s  :     0.21)                              
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.0 Energy:   -0.51096 au                                                                      
               (   1 O   1p+1:    -0.63) (   1 O   2p+1:    -0.45)                                                        
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.0 Energy:    0.03562 au                                                                      
               (   1 O   3s  :     0.23) (   1 O   4s  :     1.58) (   1 O   3p0 :     0.21)                              
               (   2 H   2s  :    -0.49) (   2 H   3s  :    -0.84) (   3 H   2s  :    -0.49)                              
               (   3 H   3s  :    -0.84)                                                                                  
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.0 Energy:    0.05807 au                                                                      
               (   1 O   3p-1:     0.84) (   2 H   2s  :    -0.57) (   2 H   3s  :    -3.04)                              
               (   3 H   2s  :     0.57) (   3 H   3s  :     3.04)                                                        
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.0 Energy:    0.17348 au                                                                      
               (   1 O   3s  :    -0.38) (   1 O   4s  :    -4.12) (   1 O   1p0 :     0.17)                              
               (   1 O   3p0 :    -1.93) (   2 H   2s  :     1.79) (   2 H   3s  :     0.62)                              
               (   2 H   2p-1:    -0.31) (   3 H   2s  :     1.79) (   3 H   3s  :     0.62)                              
               (   3 H   2p-1:     0.31)                                                                                  
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.0 Energy:    0.19602 au                                                                      
               (   1 O   1p+1:     0.17) (   1 O   2p+1:     0.24) (   1 O   3p+1:    -1.45)                              
               (   2 H   2p+1:     0.29) (   3 H   2p+1:     0.29)                                                        
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.0 Energy:    0.22379 au                                                                      
               (   1 O   2s  :     0.19) (   1 O   4s  :    -3.04) (   1 O   2p0 :    -0.18)                              
               (   1 O   3p0 :     0.47) (   2 H   2s  :     0.56) (   2 H   3s  :     0.69)                              
               (   2 H   2p-1:    -0.43) (   2 H   2p0 :    -0.50) (   3 H   2s  :     0.56)                              
               (   3 H   3s  :     0.69) (   3 H   2p-1:     0.43) (   3 H   2p0 :    -0.50)

We sample the SCF tensors on a cubic grid to obtain a volume representation (cube file format) that can be easily plotted in real-space. This is achieved through the VisualizationDriver object and the gen_cubes method.

In this example, we save cube files for the highest occupied (HOMO), lowest unoccupied molecular orbitals (LUMO), and the \(\alpha\)-spin one-particle density. The input dictionary argument cube_dict recognizes the keys cubes and files as comma-separated strings:

  • cubes, specifies which cube files to generate. In our case we ask for 2 molecular orbitals, mo(homo) and mo(lumo), and the \(\alpha\)-spin one-particle density density(alpha).

  • files, specifies the names of the generated cube files.

# initialize visualization driver
visdrv = vlx.VisualizationDriver()

# generate cube files
visdrv.gen_cubes(cube_dict={
    "cubes": "mo(homo),mo(lumo),density(alpha)", 
    "files": "HOMO.cube,LUMO.cube,density_a.cube"}, 
    molecule=h2o, basis=basis, mol_orbs=scfdrv.mol_orbs, density=scfdrv.density)

With the cube files ready, we can use the interactive py3Dmol visualization package to plot them. We initialize a view and include a stick representation of the water molecule. The view is interactive, so you can rotate, translate, zoom in and out.

import py3Dmol as p3d

# generate view
v = p3d.view(width=400, height=400)

v.addModel(h2o_xyz, "xyz")
v.setStyle({'stick':{}})
v.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

We can add additional elements on top of this view, for example volumetric data, as saved in the cube files produced from VeloxChem. For example, the HOMO. We read in the data. Note that this will load the entire file into memory and it could be unfeasible for large systems and/or with low amount of memory available.

with open("HOMO.cube", "r") as f:
    cube = f.read()
    
# negative lobe    
v.addVolumetricData(cube, "cube", {"isoval": -0.02, "color": "blue", "opacity": 0.75})
# positive lobe
v.addVolumetricData(cube, "cube", {"isoval": 0.02, "color": "red", "opacity": 0.75})

v.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Obtaining integrals from VeloxChem

We can access integrals in the atomic orbital (AO) and the molecular orbital (MO) bases. The VeloxChem API returns tensors in its internal datatypes and these can be converted to NumPy arrays as needed.

The AO basis integral drivers

There are multiple AO-basis integral drivers, one for each class of available one-electron operators and for the electron repulsion integrals:

  • vlx.OverlapIntegralsDriver for overlap integrals: \(S_{\mu\nu} = \langle \chi_{\mu} | \chi_{\nu} \rangle\).

  • vlx.KineticEnergyIntegralsDriver for kinetic energy integrals: \(T_{\mu\nu} = \langle \chi_{\mu} | -\frac{1}{2}\nabla^{2}| \chi_{\nu} \rangle\).

  • vlx.NuclearPotentialIntegralsDriver for nuclear attraction integrals: \(V_{\mu\nu} = \left\langle \chi_{\mu} \left| -\sum_{A=1}^{N_\mathrm{at}} \frac{Z_{A}}{|\mathbf{R}_A - \mathbf{r}|} \right| \chi_{\nu} \right\rangle\).

  • vlx.ElectronRepulsionIntegralsDriver for electron repulsion integrals: \((\mu\nu|\kappa\lambda) = \left\langle \chi_{\mu}(\mathbf{r})\chi_{\nu}(\mathbf{r}) \left| \frac{1}{|\mathbf{r} - \mathbf{r}^{\prime}|} \right| \chi_{\kappa}(\mathbf{r}^{\prime})\chi_{\lambda}(\mathbf{r}^{\prime}) \right\rangle\). These are stored in chemist’s notation.

More one-electron integral classes are available for specific properties, such as electric dipole moment integrals.

For one-electron integrals, we start by initializing the driver of interest and then call the compute method with molecule and basis. The to_numpy() method will convert the integrals matrix to a NumPy array.

ovldrv = vlx.OverlapIntegralsDriver()
ovl = ovldrv.compute(h2o, basis)

S = ovl.to_numpy()
print(S.shape)
print(S)

(41, 41)
[[ 1.         -0.21406265  0.19438415 ...  0.          0.
   0.        ]
 [-0.21406265  1.          0.70860733 ...  0.          0.
   0.        ]
 [ 0.19438415  0.70860733  1.         ...  0.          0.
   0.        ]
 ...
 [ 0.          0.          0.         ...  1.          0.
   0.        ]
 [ 0.          0.          0.         ...  0.          1.
   0.51422621]
 [ 0.          0.          0.         ...  0.          0.51422621
   1.        ]]

For electron repulsion integrals, we have two options:

  1. Computing the whole electron repulsion integral 4-index tensor at once and holding it in memory as a NumPy array:

def compute_in_mem(self: ElectronRepulsionIntegralsDriver, molecule: Molecule, basis: MolecularBasis) -> numpy.ndarray[numpy.float64]

  1. Computing the integrals in batches and immediately contracting with a AO density matrix. This is done in two steps:

    • We request integral screening information with the overloaded method compute:

      def compute(self: ElectronRepulsionIntegralsDriver, screening: ericut, threshold: float, molecule: Molecule, ao_basis: MolecularBasis) -> ScreeningContainer

    • We feed the ScreeningContainer and AODensityMatrix objects to the overloaded method compute. The AOFockMatrix object is an input-output parameter:

      def compute(self: ElectronRepulsionIntegralsDriver, ao_fock: AOFockMatrix, ao_density: AODensityMatrix, molecule: Molecule, ao_basis: MolecularBasis, screening: ScreeningContainer) -> None

The second option is of course preferable in terms of memory consumption. In this workshop, we will use the first option as it maps directly to expressions on paper.

eridrv = vlx.ElectronRepulsionIntegralsDriver()
eri = eridrv.compute_in_mem(h2o, basis)
print(eri.shape)

(41, 41, 41, 41)

The MO basis integral driver

For one-electron integrals, the transformation from AO basis to MO basis is:

\[ \mathbf{O}_{\mathrm{MO}} = \mathbf{C}^{t} \mathbf{O}_{\mathrm{AO}} \mathbf{C}, \quad O_{\mathrm{MO}, pq}= \sum_{\mu\nu}C_{\mu p}O_{\mathrm{AO},\mu\nu}C_{\nu q} \]

and we can straightforwardly use NumPy to achieve it:

import numpy as np

C = scfdrv.mol_orbs.alpha_to_numpy()
S_mo = np.einsum("mP,mn,nQ->PQ", C, S, C)
print(S_mo)

[[ 1.00000000e+00 -1.36025021e-17 -3.57964920e-17 ... -1.74283951e-16
  -1.10764804e-16 -2.49608341e-16]
 [-2.74853721e-17  1.00000000e+00  5.06622565e-17 ...  1.04477795e-16
  -6.61797006e-16 -3.19497969e-17]
 [-2.19011408e-17  8.05257573e-17  1.00000000e+00 ... -6.56123751e-18
  -2.76931435e-16 -1.84921523e-15]
 ...
 [-1.74283951e-16  1.04477795e-16 -6.56123751e-18 ...  1.00000000e+00
  -3.08015996e-16  8.33345767e-16]
 [-3.98918637e-16 -2.33970829e-15  3.70130513e-16 ... -3.08015996e-16
   1.00000000e+00  2.24403297e-15]
 [-4.38633966e-16  4.56586120e-16 -1.30971622e-15 ...  8.33345767e-16
  -1.44417777e-16  1.00000000e+00]]

For two electron integrals, the transformation reads:

\[ \langle pq | rs \rangle = \sum_{\mu\nu\kappa\lambda} C_{\mu p}C_{\nu r} (\mu\nu|\kappa\lambda) C_{\kappa q} C_{\lambda s}, \]

with the MO integrals in physicists’ notation. The transformation requires \(O(N^{8})\) operation count.

However, we can perform it more efficiently as a stepwise contraction:

\[ \langle pq | rs \rangle = \sum_{\mu} C_{\mu p} \left(\sum_{\nu} C_{\nu r} \left (\sum_{\kappa} \left(\sum_{\lambda} (\mu\nu|\kappa\lambda) C_{\lambda s} \right) C_{\kappa q} \right)\right). \]

This is implemented in the vlx.MOIntegralsDriver. Similarly to the AO basis ERI evaluation, we will use the compute_in_mem method:

moeridrv = vlx.MOIntegralsDriver()
help(moeridrv.compute_in_mem)

Help on method compute_in_mem in module veloxchem.mointsdriver:

compute_in_mem(molecule, basis, mol_orbs, mints_type) method of veloxchem.mointsdriver.MOIntegralsDriver instance
    Performs in-memory MO integrals calculation for a molecule and a basis
    set.
    
    :param molecule:
        The molecule.
    :param basis:
        The AO basis set.
    :param mol_orbs:
        The molecular orbitals.
    :param mints_type:
        The type of MO integrals to be calculated.
    
    :return:
        The computed MO integrals.

This accepts a mints_type argument to specify which combination of orbitals, occupied or virtual, to use in the transformation:

oooo = moeridrv.compute_in_mem(h2o, basis, scfdrv.mol_orbs, "oooo")
print(oooo.shape)

oovo = moeridrv.compute_in_mem(h2o, basis, scfdrv.mol_orbs, "oovo")
print(oovo.shape)

vvvv = moeridrv.compute_in_mem(h2o, basis, scfdrv.mol_orbs, "vvvv")
print(vvvv.shape)

(5, 5, 5, 5)
(5, 5, 36, 5)
(36, 36, 36, 36)

previous

First steps with VeloxChem

next

The Roothaan-Hall self-consistent field procedure