First steps with VeloxChem
Contents
First steps with VeloxChem¶
We want to run a restricted Hartree-Fock test calculation on a water molecule. We will use the following structure.
import py3Dmol as p3d
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
"""
v = p3d.view(width=200, height=200)
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
Running with an input file¶
Let us save the following as water.inp:
@jobs
task: scf
@end
@method settings
basis: aug-cc-pvdz
@end
@molecule
charge: 0
multiplicity: 1
xyz:
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.740848095288 0.582094932012
H 0.000000000000 -0.740848095288 0.582094932012
@end
A VeloxChem input file has:
sections. Each section starts with
@<section-name>and is closed with@end,keywords. You specify them with
<keyword-name>: <keyword-value>.
Documentation for all input options is available at: https://docs.veloxchem.org/inputs/keywords.html
The command-line interface¶
Possibly the most common way of running a quantum chemistry code is by using a command-line interface (CLI). You can access the VeloxChem CLI with:
$ vlx --help
this will print the following help text:
usage:
================= VeloxChem =================
vlx input_file [output_file]
or
python -m veloxchem input_file [output_file]
positional arguments:
input_file Input file
output_file Output file (default: STDOUT)
optional arguments:
-h, --help show this help message and exit
-v, --version show program's version number and exit
Let’s run our water calculation with:
$ vlx water.inp water.out
such that the output will be saved to the plain-text file water.out. Congratulations! You have completed your first quantum chemical calculation with VeloxChem!
The output file¶
Let us look now at the information contained in the output file. VeloxChem logs in the output stream the operations performed during a calculation. You will notice:
The molecular geometry
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 charge : 0
Spin multiplicity : 1
Number of atoms : 3
Number of alpha electrons : 5
Number of beta electrons : 5
The molecular basis and self-consistent field setup
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
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
The report on the iterative SCF procedure
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.041697549819 -0.0000000000 0.00000039 0.00000005 0.00000148
The final summary of SCF energies
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
The molecular orbitals and the electric dipole moment
Last, but not least, the preferred citation for VeloxChem
!========================================================================================================================!
! Rinkevicius, Z.; Li, X.; Vahtras, O.; Ahmadzadeh, K.; Brand, M.; Ringholm, M.; !
! List, N. H.; Scheurer, M.; Scott, M.; Dreuw, A.; Norman, P. !
! VeloxChem: A Python-driven Density-functional Theory Program for Spectroscopy !
! Simulations in High-performance Computing Environments. !
! WIREs Comput Mol Sci 2020, 10 (5), e1457. !
!========================================================================================================================!
Running with a Python script¶
Extracting information from a plain-text output file is not always desirable, for a number of reasons:
The format of the output is not guaranteed to stay the same between different versions of the code. Data-harvesting scripts might fail unexpectedly.
Reading floating-point numbers from their text representation might introduce subtle round-off errors.
VeloxChem offers a powerful Python application programming interface (API) to programmatically access computed data, for example, SCF ground state electric dipole moments, and computational primitives, such as AO basis integrals. Let’s see how we can write a Python script to perform a calculation on the water molecule without using plain-text input and output files. First and foremost, we import the VeloxChem API:
import veloxchem as vlx
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.
The Molecule object¶
The next task is to create a Molecule object. We can use the h2o_xyz string defined above directly.
h2o = vlx.Molecule.from_xyz_string(h2o_xyz)
print(h2o.get_string())
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
The help native method of Python works as well and we can explore the VeloxChem API easily. For example, h2o.get_coordinates() will return the molecular structure as a NumPy array.
help(vlx.Molecule)
h2o.get_coordinates()
Help on class Molecule in module veloxchem.veloxchemlib:
class Molecule(pybind11_builtins.pybind11_object)
| Method resolution order:
| Molecule
| pybind11_builtins.pybind11_object
| builtins.object
|
| Methods defined here:
|
| __eq__(...)
| __eq__(self: veloxchem.veloxchemlib.Molecule, arg0: veloxchem.veloxchemlib.Molecule) -> bool
|
| __init__(...)
| __init__(*args, **kwargs)
| Overloaded function.
|
| 1. __init__(self: veloxchem.veloxchemlib.Molecule) -> None
|
| 2. __init__(self: veloxchem.veloxchemlib.Molecule, arg0: veloxchem.veloxchemlib.Molecule) -> None
|
| 3. __init__(self: veloxchem.veloxchemlib.Molecule, arg0: veloxchem.veloxchemlib.Molecule, arg1: veloxchem.veloxchemlib.Molecule) -> None
|
| 4. __init__(self: veloxchem.veloxchemlib.Molecule, symbols: List[str], coordinates: numpy.ndarray[numpy.float64], units: str = 'angstrom') -> None
|
| 5. __init__(self: veloxchem.veloxchemlib.Molecule, Zs: List[int], coordinates: numpy.ndarray[numpy.float64], units: str = 'angstrom') -> None
|
| broadcast(...)
| broadcast(self: veloxchem.veloxchemlib.Molecule, arg0: int, arg1: object) -> None
|
| center_of_mass = _Molecule_center_of_mass(self)
| Computes center of mass of a molecule.
|
| :return:
| The center of mass.
|
| check_multiplicity(...)
| check_multiplicity(self: veloxchem.veloxchemlib.Molecule) -> None
|
| check_proximity(...)
| check_proximity(self: veloxchem.veloxchemlib.Molecule, arg0: float) -> None
|
| coordination_numbers(...)
| coordination_numbers(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| elem_ids_to_numpy(...)
| elem_ids_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.int32]
|
| get_charge(...)
| get_charge(self: veloxchem.veloxchemlib.Molecule) -> float
|
| get_coordinates = _Molecule_get_coordinates(self)
| Returns atom coordinates in atomic units.
|
| :return:
| A numpy array of atom coordinates (nx3) in atomic units.
|
| get_elemental_composition(...)
| get_elemental_composition(self: veloxchem.veloxchemlib.Molecule) -> list
|
| get_ic_rmsd = _Molecule_get_ic_rmsd(self, ref_mol)
| Gets statistical deviation of bonds, angles and dihedral angles between
| self and reference geometry.
|
| :param ref_mol:
| The reference molecule (or xyz filename).
|
| get_labels = _Molecule_get_labels(self)
| Returns atom labels.
|
| :return:
| A list of atom labels.
|
| get_multiplicity(...)
| get_multiplicity(self: veloxchem.veloxchemlib.Molecule) -> int
|
| get_string(...)
| get_string(self: veloxchem.veloxchemlib.Molecule) -> str
|
| get_sub_molecule(...)
| get_sub_molecule(self: veloxchem.veloxchemlib.Molecule, arg0: int, arg1: int) -> veloxchem.veloxchemlib.Molecule
|
| masses_to_numpy(...)
| masses_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| more_info = _Molecule_more_info(self)
| Returns more information about the molecule.
|
| :return:
| Molecular information in plain text.
|
| nuclear_repulsion_energy(...)
| nuclear_repulsion_energy(self: veloxchem.veloxchemlib.Molecule) -> float
|
| number_of_alpha_electrons(...)
| number_of_alpha_electrons(self: veloxchem.veloxchemlib.Molecule) -> int
|
| number_of_atoms(...)
| number_of_atoms(*args, **kwargs)
| Overloaded function.
|
| 1. number_of_atoms(self: veloxchem.veloxchemlib.Molecule) -> int
|
| 2. number_of_atoms(self: veloxchem.veloxchemlib.Molecule, idElemental: int) -> int
|
| 3. number_of_atoms(self: veloxchem.veloxchemlib.Molecule, iatom: int, natoms: int, idElemental: int) -> int
|
| number_of_beta_electrons(...)
| number_of_beta_electrons(self: veloxchem.veloxchemlib.Molecule) -> int
|
| number_of_electrons(...)
| number_of_electrons(self: veloxchem.veloxchemlib.Molecule) -> int
|
| partial_charges(...)
| partial_charges(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| set_charge(...)
| set_charge(self: veloxchem.veloxchemlib.Molecule, arg0: float) -> None
|
| set_multiplicity(...)
| set_multiplicity(self: veloxchem.veloxchemlib.Molecule, arg0: int) -> None
|
| vdw_radii_to_numpy(...)
| vdw_radii_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| write_xyz = _Molecule_write_xyz(self, xyz_filename)
| Writes molecular geometry to xyz file.
|
| :param xyz_filename:
| The name of the xyz file.
|
| x_to_numpy(...)
| x_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| y_to_numpy(...)
| y_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| z_to_numpy(...)
| z_to_numpy(self: veloxchem.veloxchemlib.Molecule) -> numpy.ndarray[numpy.float64]
|
| ----------------------------------------------------------------------
| Static methods defined here:
|
| from_dict = _Molecule_from_dict(mol_dict)
| Reads molecule from a dictionary.
|
| :param mol_dict:
| The molecule dictionary.
|
| :return:
| The molecule.
|
| from_xyz_string = _Molecule_from_xyz_string(xyz)
| Generate molecule from string in XYZ format.
|
| :param xyz:
| String with XYZ structure.
|
| :return:
| The molecule.
|
| read_str = _Molecule_read_str(xyzstr, units='angstrom')
| Reads molecule from xyz string.
|
| :param xyzstr:
| The xyz string.
| :param units:
| The unit of coordinates.
|
| :return:
| The molecule.
|
| read_xyz = _Molecule_read_xyz(xyzfile)
| Reads molecule from file in XYZ format.
|
| :param xyzfile:
| File with molecular structure in XYZ format.
|
| :return:
| The molecule.
|
| ----------------------------------------------------------------------
| Data and other attributes defined here:
|
| __hash__ = None
|
| ----------------------------------------------------------------------
| Static methods inherited from pybind11_builtins.pybind11_object:
|
| __new__(*args, **kwargs) from pybind11_builtins.pybind11_type
| Create and return a new object. See help(type) for accurate signature.
array([[ 0. , 0. , 0. ],
[ 0. , 1.4, 1.1],
[ 0. , -1.4, 1.1]])
Gaussian molecular basis sets¶
Once we have a molecule, we can define a Gaussian molecular basis set. vlx.MolecularBasis.get_avail_basis shows the list of already available basis sets in VeloxChem for a given chemical element.
vlx.MolecularBasis.get_avail_basis("C")
['6-31++G',
'6-31+G',
'6-31+G(D,P)',
'6-311++G',
'6-311++G(2D,2P)',
'6-311++G(D,P)',
'6-311+G',
'6-311+G(2D,P)',
'6-311G',
'6-311G(2DF,2PD)',
'6-31G',
'6-31G(2DF,P)',
'AUG-CC-PCVDZ',
'AUG-CC-PCVQZ',
'AUG-CC-PCVTZ',
'AUG-CC-PVDZ',
'AUG-CC-PVTZ',
'CC-PCVDZ',
'CC-PVDZ',
'CC-PVDZ-RI',
'CC-PVTZ',
'CC-PVTZ-RI',
'D-AUG-CC-PVDZ',
'D-AUG-CC-PVQZ',
'D-AUG-CC-PVTZ',
'DEF2-QZVP',
'DEF2-QZVPD',
'DEF2-QZVPP',
'DEF2-QZVPPD',
'DEF2-RI-J',
'DEF2-SV(P)',
'DEF2-SVP',
'DEF2-SVPD',
'DEF2-TZVP',
'DEF2-TZVPP',
'DEF2-TZVPPD',
'MIN-CC-PVDZ',
'SADLEJ-PVTZ',
'STO-3G',
'STO-6G']
We generate a basis object for our molecule using one of the bases included in VeloxChem with:
basis = vlx.MolecularBasis.read(h2o, "sto-3g")
print(basis.get_string(h2o))
Molecular Basis (Atomic Basis)
================================
Basis: STO-3G
Atom Contracted GTOs Primitive GTOs
H (1S) (3S)
O (2S,1P) (6S,3P)
Contracted Basis Functions : 7
Primitive Basis Functions : 21
Self-consistent field calculation drivers¶
You launch restricted or unrestricted SCF calculations through their respective driver object: vlx.ScfRestrictedDriver and vlx.ScfUnrestrictedDriver.
scfdrv = vlx.ScfRestrictedDriver()
You can configure the drivers with options passed as a Python dictionary.
settings = {"conv_thresh": 1e-6}
scfdrv.update_settings(settings)
Finally, we run the calculation with the compute method. The output will be exactly identical to the one obtained from running the input file.
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.18 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.09 sec.
* Info * SAD initial guess computed in 0.00 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -74.960337067134 au. Time: 0.48 sec.
* Info * Overlap matrix computed in 0.09 sec.
* Info * Kinetic energy matrix computed in 0.01 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 -74.960337068856 0.0000000000 0.00002391 0.00000577 0.00000000
2 -74.960337069026 -0.0000000002 0.00000785 0.00000178 0.00001947
3 -74.960337069049 -0.0000000000 0.00000050 0.00000012 0.00001040
*** SCF converged in 3 iterations. Time: 0.49 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -74.9603370690 au
Electronic Energy : -84.3039752270 au
Nuclear Repulsion Energy : 9.3436381580 au
------------------------------------
Gradient Norm : 0.0000004994 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.24119 au
( 1 O 1s : -0.99)
Molecular Orbital No. 2:
--------------------------
Occupation: 2.0 Energy: -1.27737 au
( 1 O 1s : 0.23) ( 1 O 2s : -0.83) ( 2 H 1s : -0.16)
( 3 H 1s : -0.16)
Molecular Orbital No. 3:
--------------------------
Occupation: 2.0 Energy: -0.62456 au
( 1 O 1p-1: -0.61) ( 2 H 1s : -0.44) ( 3 H 1s : 0.44)
Molecular Orbital No. 4:
--------------------------
Occupation: 2.0 Energy: -0.45684 au
( 1 O 2s : 0.54) ( 1 O 1p0 : -0.78) ( 2 H 1s : -0.27)
( 3 H 1s : -0.27)
Molecular Orbital No. 5:
--------------------------
Occupation: 2.0 Energy: -0.39296 au
( 1 O 1p+1: 1.00)
Molecular Orbital No. 6:
--------------------------
Occupation: 0.0 Energy: 0.62228 au
( 1 O 2s : -0.91) ( 1 O 1p0 : -0.75) ( 2 H 1s : 0.81)
( 3 H 1s : 0.81)
Molecular Orbital No. 7:
--------------------------
Occupation: 0.0 Energy: 0.75803 au
( 1 O 1p-1: -1.00) ( 2 H 1s : 0.85) ( 3 H 1s : -0.85)
The driver object contains a wealth of useful information:
The final SCF energy:
scfdrv.get_scf_energy()Various tensors, such as the molecular orbital coefficients, AO density and Fock matrices:
scfdrv.scf_tensors.
The scf_tensors property is a Python dictionary with keys:
C, for the molecular orbital coefficients.E, for the orbital energies.S, for the overlap integral.D, for the \(\alpha\) and \(\beta\) density matrices as a tuple of NumPy arrays.F, for the \(\alpha\) and \(\beta\) Fock matrices as a tuple of NumPy arrays.
scfdrv.scf_tensors.keys()
dict_keys(['C', 'E', 'S', 'D', 'F'])
scfdrv.scf_tensors["D"]
(array([[ 1.05303213e+00, -2.22302631e-01, -1.46417671e-02,
-1.46417671e-02, -2.04228488e-16, 5.41902851e-02,
-7.42025150e-18],
[-2.22302631e-01, 9.82076521e-01, -1.64989375e-02,
-1.64989375e-02, -2.16731540e-16, -3.10480924e-01,
6.28860220e-17],
[-1.46417671e-02, -1.64989375e-02, 2.97409605e-01,
-9.77472101e-02, 2.70379139e-01, 2.33638805e-01,
1.91065463e-16],
[-1.46417671e-02, -1.64989375e-02, -9.77472101e-02,
2.97409605e-01, -2.70379139e-01, 2.33638805e-01,
-3.43153996e-16],
[-2.04228488e-16, -2.16731540e-16, 2.70379139e-01,
-2.70379139e-01, 3.70004394e-01, -1.39439910e-15,
1.06935025e-16],
[ 5.41902851e-02, -3.10480924e-01, 2.33638805e-01,
2.33638805e-01, -1.39439910e-15, 6.22944826e-01,
-7.67046565e-17],
[-7.42025150e-18, 6.28860220e-17, 1.91065463e-16,
-3.43153996e-16, 1.06935025e-16, -7.67046565e-17,
1.00000000e+00]]),
array([[ 1.05303213e+00, -2.22302631e-01, -1.46417671e-02,
-1.46417671e-02, -2.04228488e-16, 5.41902851e-02,
-7.42025150e-18],
[-2.22302631e-01, 9.82076521e-01, -1.64989375e-02,
-1.64989375e-02, -2.16731540e-16, -3.10480924e-01,
6.28860220e-17],
[-1.46417671e-02, -1.64989375e-02, 2.97409605e-01,
-9.77472101e-02, 2.70379139e-01, 2.33638805e-01,
1.91065463e-16],
[-1.46417671e-02, -1.64989375e-02, -9.77472101e-02,
2.97409605e-01, -2.70379139e-01, 2.33638805e-01,
-3.43153996e-16],
[-2.04228488e-16, -2.16731540e-16, 2.70379139e-01,
-2.70379139e-01, 3.70004394e-01, -1.39439910e-15,
1.06935025e-16],
[ 5.41902851e-02, -3.10480924e-01, 2.33638805e-01,
2.33638805e-01, -1.39439910e-15, 6.22944826e-01,
-7.67046565e-17],
[-7.42025150e-18, 6.28860220e-17, 1.91065463e-16,
-3.43153996e-16, 1.06935025e-16, -7.67046565e-17,
1.00000000e+00]]))