Geometry

The first step of any quantum chemistry calculation is specification of the system geometry. The InQuanto GeometryMolecular and GeometryPeriodic classes exist to standardize geometry formats and provide convenient functionality for building and manipulating molecular and periodic structures.

Molecular Systems

Initializing Structures

Geometry objects may be constructed in several ways. Atom-by-atom construction is possible with the add_atom() method:

from inquanto.geometries import GeometryMolecular

g = GeometryMolecular()
g.add_atom("H", [0, 0, 0])
g.add_atom("H", [0, 0, 0.735])
g.add_atom("H", [0, 0, 1.])
print(g.df)
   element  x  y      z
id                     
0        H  0  0  0.000
1        H  0  0  0.735
2        H  0  0  1.000

where the dataframe attribute is a pandas.DataFrame of the geometric structure. If the atom position is left empty in add_atom(), random coordinates between ±1 Angstroms are generated. See below:

g = GeometryMolecular()
g.add_atom("H")
print(g.df)
   element         x         y         z
id                                      
0        H  0.816403 -0.913649  0.033713

Note

InQuanto Geometry objects support Angstrom (default) and Bohr as distance units, specified by the distance_units constructor argument. If another unit is used, the user can convert the geometry to Bohrs or Angstroms using the rescale_position_vectors() class method.

Alternatively, one may instantiate a GeometryMolecular object from z-matrices or Cartesian coordinates provided they are appropriately formatted. Examples for both are given below.

water_zmatrix = """h
o 1 1.0
h 2 1.0 1 104.5
"""

g = GeometryMolecular(water_zmatrix)
print("Geometry from z-matrix:\n", g.df)

water_xyz = [
     ["O", [0.0000000, 0.0000000, 0.1271610]],
     ["H", [0.0000000, 0.7580820, -0.5086420]],
     ["H", [0.0000000, -0.7580820, -0.5086420]],
]
g = GeometryMolecular(water_xyz)
print("\nGeometry from xyz:\n", g.df)
Geometry from z-matrix:
    element    x         y        z
id                                
0        h  0.0  0.000000  0.00000
1        o  0.0  0.000000  1.00000
2        h  0.0  0.968148  1.25038

Geometry from xyz:
    element    x         y         z atom
id                                      
0        O  0.0  0.000000  0.127161   O1
1        H  0.0  0.758082 -0.508642   H2
2        H  0.0 -0.758082 -0.508642   H3

It is also possible to instantiate/write from files for a limited number of formats without the use of any extensions. These include:

Calculating and Modifying Properties

A variety of geometrical quantities, including bond lengths, angles and dihedrals, can be calculated with Geometry class methods:

ethane_xyz = [
     ["C", [0.0000000, 0.0000000, 0.7688350]],
     ["C", [0.0000000, 0.0000000, -0.7688350]],
     ["H", [0.0000000, 1.0157000, 1.1533240]],
     ["H", [-0.8796220, -0.5078500, 1.1533240]],
     ["H", [0.8796220, -0.5078500, 1.1533240]],
     ["H", [0.0000000, -1.0157000, -1.1533240]],
     ["H", [-0.8796220, 0.5078500, -1.1533240]],
     ["H", [0.8796220, 0.5078500, -1.1533240]],
 ]
g = GeometryMolecular(ethane_xyz)
print("C=C bond length:", g.bond_length([0, 1]))
print("H-C-C bond angle:", g.bond_angle([2, 0, 1]))
print("2-0-1-5 dihedral angle:", g.dihedral_angle([2, 0, 1, 5]))

g = GeometryMolecular(water_xyz)
print("\n water distance matrix:\n", g.compute_distance_matrix())
C=C bond length: 1.53767
H-C-C bond angle: 110.73395111184878
2-0-1-5 dihedral angle: 180.0

 water distance matrix:
 [[0.         0.98941082 0.98941082]
 [0.98941082 0.         1.516164  ]
 [0.98941082 1.516164   0.        ]]

These quantities can also be modified:

g = GeometryMolecular(water_xyz)
print("Equilibrium structure 0-1 bond length:", g.bond_length([0, 1]))
print("Equilibrium structure 0-1 bond angle:", g.bond_angle([1, 0, 2]))

# modifying a bond length
g.modify_bond_length(atom_ids=[0, 1], new_bond_length=2)
print("Stretched structure 0-1 bond length:", g.bond_length([0, 1]))

# modifying a bond angle
g.modify_bond_angle(atom_ids=[1, 0, 2], theta=90, units="deg")
print("Closed bond angle 2-0-2:", g.bond_angle([1, 0, 2]))
Equilibrium structure 0-1 bond length: 0.989410821414947
Equilibrium structure 0-1 bond angle: 100.0269116572392
Stretched structure 0-1 bond length: 2.0
Closed bond angle 2-0-2: 90.0

Note

Only the last specified atom in the atom_ids has its position changed when modifying geometric properties as above. The position of all other atoms remains unchanged.

Calculating and Modifying Properties By Group

To transform the geometry by moving more than one atom at a time, one may define a grouping scheme within a structure. A chemically relevant example is switching between the staggered and eclipsed geometries of ethane. Below we define a "methyl" grouping scheme which contains two sub-groups:

g = GeometryMolecular(ethane_xyz)
g.set_groups("methyl", {"ch3_1": [0, 2, 3, 4], "ch3_2": [1, 5, 6, 7]})
print(g.df)
   element         x        y         z atom methyl
id                                                 
0        C  0.000000  0.00000  0.768835   C1  ch3_1
1        C  0.000000  0.00000 -0.768835   C2  ch3_2
2        H  0.000000  1.01570  1.153324   H3  ch3_1
3        H -0.879622 -0.50785  1.153324   H4  ch3_1
4        H  0.879622 -0.50785  1.153324   H5  ch3_1
5        H  0.000000 -1.01570 -1.153324   H6  ch3_2
6        H -0.879622  0.50785 -1.153324   H7  ch3_2
7        H  0.879622  0.50785 -1.153324   H8  ch3_2

Transformations may then be performed “by group”, so that all atoms within a sub-group are modified simultaneously, preserving internal structure. For example, below we stretch the C-C bond, then rotate the atoms in "ch3_2" such that the ethane is in an eclipsed configuration:

g.modify_bond_length_by_group(atom_ids=[0, 1], bond_length=1.8, group="methyl")
g.modify_dihedral_angle_by_group(atom_ids=[2, 0, 1, 7], theta=0.0, group="methyl")
print(g.df)
   element             x        y         z atom methyl
id                                                     
0        C  0.000000e+00  0.00000  0.768835   C1  ch3_1
1        C  0.000000e+00  0.00000 -1.031165   C2  ch3_2
2        H  0.000000e+00  1.01570  1.153324   H3  ch3_1
3        H -8.796220e-01 -0.50785  1.153324   H4  ch3_1
4        H  8.796220e-01 -0.50785  1.153324   H5  ch3_1
5        H  8.796220e-01 -0.50785 -1.415654   H6  ch3_2
6        H -8.796220e-01 -0.50785 -1.415654   H7  ch3_2
7        H  3.349405e-16  1.01570 -1.415654   H8  ch3_2

Similarly to single atom modifications, the sub-group that gets moved when transformations are applied is that which contains the last specified atom in atom_ids.

Generating Many Structures

Convenience methods for generating a series of structures with varying geometric properties are also available. For example, below we generate a range of GeometryMolecular objects for water with different H-O-H bond angles:

water_geom = GeometryMolecular(water_xyz)
geometries = water_geom.scan_bond_angle([1, 0, 2], [100, 101, 103, 104, 105, 106, 107, 108])
print([g.bond_angle([1,0,2]) for g in geometries])
[100.0, 101.0, 103.0, 103.99999999999999, 105.0, 106.0, 107.0, 108.0]

Similar functions are implemented for scanning over bond lengths and dihedrals (scan_bond_length() and scan_dihedral_angle()), and equivalently for varying properties by group (scan_bond_angle_by_group() for example).

Periodic systems

For periodic systems there exists the GeometryPeriodic class, which holds unit cell lattice vectors in addition to atomic positions:

from inquanto.geometries import GeometryPeriodic
import numpy as np

# Diamond fcc
d = 3.576

a, b, c = [
   np.array([0, d/2, d/2]),
   np.array([d/2, 0, d/2]),
   np.array([d/2, d/2, 0])
]
atoms = [
   ["C", [0, 0, 0]],
   ["C", a/4 + b/4 + c/4]
]

cfcc_geom = GeometryPeriodic(geometry=atoms, unit_cell=[a, b, c])
print("Diamond fcc lattice vectors:\n", cfcc_geom.unit_cell)
print("\nCartesian atomic positions:\n", cfcc_geom.df)
Diamond fcc lattice vectors:
 [[0.    1.788 1.788]
 [1.788 0.    1.788]
 [1.788 1.788 0.   ]]

Cartesian atomic positions:
    element      x      y      z atom
id                                  
0        C  0.000  0.000  0.000   C1
1        C  0.894  0.894  0.894   C2

GeometryPeriodic supports all of the methods discussed above, excluding instantiation by z-matrices. In addition, one may easily construct supercells with the build_supercell() method, which edits the geometry in-place:

cfcc_geom.build_supercell(dimensions=[2, 2, 1])
print("Supercell atomic positions:\n", cfcc_geom.df)
Supercell atomic positions:
    element      x      y      z atom
id                                  
0        C  0.000  0.000  0.000   C1
1        C  0.894  0.894  0.894   C2
2        C  1.788  0.000  1.788   C3
3        C  2.682  0.894  2.682   C4
4        C  0.000  1.788  1.788   C5
5        C  0.894  2.682  2.682   C6
6        C  1.788  1.788  3.576   C7
7        C  2.682  2.682  4.470   C8

To visualize both molecular and periodic structures, see the inquanto-nglview extension.