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