Simple properties#
Simple properties do not require a DFT calculation. They are typically only functions of the atom types and geometries.
Getting cartesian positions#
If you want the \((x,y,z)\) coordinates of the atoms, use the func:ase.Atoms.getpositions. If you are interested in the fractional coordinates, use func:ase.Atoms.getscaledpositions.
from ase.build import molecule
atoms = molecule('C6H6') # benzene
# access properties on each atom
print(' # sym p_x p_y p_z')
print('------------------------------')
for i, atom in enumerate(atoms):
print('{0:3d}{1:^4s}{2:-8.2f}{3:-8.2f}{4:-8.2f}'.format(i,
atom.symbol,
atom.x,
atom.y,
atom.z))
# get all properties in arrays
sym = atoms.get_chemical_symbols()
pos = atoms.get_positions()
num = atoms.get_atomic_numbers()
atom_indices = range(len(atoms))
print()
print(' # sym at# p_x p_y p_z')
print('-------------------------------------')
for i, s, n, p in zip(atom_indices, sym, num, pos):
px, py, pz = p
print('{0:3d}{1:>3s}{2:8d}{3:-8.2f}{4:-8.2f}{5:-8.2f}'.format(i, s, n,
px, py, pz))
# sym p_x p_y p_z
------------------------------
0 C 0.00 1.40 0.00
1 C 1.21 0.70 0.00
2 C 1.21 -0.70 0.00
3 C 0.00 -1.40 0.00
4 C -1.21 -0.70 0.00
5 C -1.21 0.70 0.00
6 H 0.00 2.48 0.00
7 H 2.15 1.24 0.00
8 H 2.15 -1.24 0.00
9 H 0.00 -2.48 0.00
10 H -2.15 -1.24 0.00
11 H -2.15 1.24 0.00
# sym at# p_x p_y p_z
-------------------------------------
0 C 6 0.00 1.40 0.00
1 C 6 1.21 0.70 0.00
2 C 6 1.21 -0.70 0.00
3 C 6 0.00 -1.40 0.00
4 C 6 -1.21 -0.70 0.00
5 C 6 -1.21 0.70 0.00
6 H 1 0.00 2.48 0.00
7 H 1 2.15 1.24 0.00
8 H 1 2.15 -1.24 0.00
9 H 1 0.00 -2.48 0.00
10 H 1 -2.15 -1.24 0.00
11 H 1 -2.15 1.24 0.00
Molecular weight and molecular formula#
[BROKEN LINK: index:molecular weight]
We can quickly compute the molecular weight of a molecule with this recipe. We use func:ase.Atoms.getmasses to get an array of the atomic masses of each atom in the Atoms object, and then just sum them up.
from ase.build import molecule
atoms = molecule('C6H6')
masses = atoms.get_masses()
molecular_weight = masses.sum()
molecular_formula = atoms.get_chemical_formula(mode='reduce')
# note use of two lines to keep length of line reasonable
s = 'The molecular weight of {0} is {1:1.2f} gm/mol'
print(s.format(molecular_formula, molecular_weight))
The molecular weight of C6H6 is 78.11 gm/mol
Note that the argument reduce=True for func:ase.Atoms.getchemicalformula collects all the symbols to provide a molecular formula.
Center of mass#
[BROKEN LINK: index:center of mass] The center of mass (COM) is defined as:
COM = \(\frac{\sum m_i \cdot r_i}{\sum m_i}\)
The center of mass is essentially the average position of the atoms, weighted by the mass of each atom. Here is an example of getting the center of mass from an Atoms object using func:ase.Atoms.getcenterofmass.
from ase.build import molecule
import numpy as np
# ammonia
atoms = molecule('NH3')
# cartesian coordinates
print('COM1 = {0}'.format(atoms.get_center_of_mass()))
# compute the center of mass by hand
pos = atoms.positions
masses = atoms.get_masses()
COM = np.array([0., 0., 0.])
for m, p in zip(masses, pos):
COM += m*p
COM /= masses.sum()
print('COM2 = {0}'.format(COM))
# one-line linear algebra definition of COM
print('COM3 = {0}'.format(np.dot(masses, pos) / np.sum(masses)))
COM1 = [1.29821427e-19 5.91861899e-08 4.75435401e-02]
COM2 = [0.00000000e+00 5.91861899e-08 4.75435401e-02]
COM3 = [1.29821427e-19 5.91861899e-08 4.75435401e-02]
You can see see that these centers of mass, which are calculated by different methods, are the same.
Moments of inertia#
[BROKEN LINK: index:moment of inertia] The moment of inertia is a measure of resistance to changes in rotation. It is defined by \(I = \sum_{i=1}^N m_i r_i^2\) where \(r_i\) is the distance to an axis of rotation. There are typically three moments of inertia, although some may be zero depending on symmetry, and others may be degenerate. There is a convenient function to get the moments of inertia: func:ase.Atoms.getmomentsofinertia. Here are several examples of molecules with different types of symmetry.:
from ase.build import molecule
print('linear rotors: I = [0 Ia Ia]')
atoms = molecule('CO2')
print(' CO2 moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
print('')
print('symmetric rotors (Ia = Ib) < Ic')
atoms = molecule('NH3')
print(' NH3 moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
atoms = molecule('C6H6')
print(' C6H6 moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
print('')
print('symmetric rotors Ia < (Ib = Ic)')
atoms = molecule('CH3Cl')
print('CH3Cl moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
print('')
print('spherical rotors Ia = Ib = Ic')
atoms = molecule('CH4')
print(' CH4 moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
print('')
print('unsymmetric rotors Ia != Ib != Ic')
atoms = molecule('C3H7Cl')
print(' C3H7Cl moments of inertia: {}'.format(atoms.get_moments_of_inertia()))
linear rotors: I = [0 Ia Ia]
CO2 moments of inertia: [ 0. 44.45273132 44.45273132]
symmetric rotors (Ia = Ib) < Ic
NH3 moments of inertia: [1.71022353 1.71022474 2.67047664]
C6H6 moments of inertia: [ 88.78025559 88.78027717 177.56053276]
symmetric rotors Ia < (Ib = Ic)
CH3Cl moments of inertia: [ 3.2039126 37.969823 37.96982492]
spherical rotors Ia = Ib = Ic
CH4 moments of inertia: [3.19164619 3.19164619 3.19164619]
unsymmetric rotors Ia != Ib != Ic
C3H7Cl moments of inertia: [ 19.41420447 213.18480664 223.1578698 ]
If you want to know the principle axes of rotation, we simply pass vectors=True to the function, and it returns the moments of inertia and the principle axes.
from ase.build import molecule
atoms = molecule('CH3Cl')
moments, axes = atoms.get_moments_of_inertia(vectors=True)
print('Moments = {0}'.format(moments))
print('axes = {0}'.format(axes))
Moments = [ 3.2039126 37.969823 37.96982492]
axes = [[0. 0. 1.]
[0. 1. 0.]
[1. 0. 0.]]
This shows the first moment is about the z-axis, the second moment is about the y-axis, and the third moment is about the x-axis.
Computing bond lengths and angles#
A typical question we might ask is, “What is the structure of a molecule?” In other words, what are the bond lengths, angles between bonds, and similar properties. The Atoms object contains an func:ase.Atoms.getdistance method to make this easy. To calculate the distance between two atoms, you have to specify their indices, remembering that the index starts at 0.
from ase.build import molecule
# ammonia
atoms = molecule('NH3')
print('atom symbol')
print('===========')
for i, atom in enumerate(atoms):
print('{0:2d} {1:3s}' .format(i, atom.symbol))
# N-H bond length
s = 'The N-H distance is {0:1.3f} angstroms'
print(s.format(atoms.get_distance(0, 1)))
atom symbol
===========
0 N
1 H
2 H
3 H
The N-H distance is 1.017 angstroms
Bond angles are a little trickier. If we had vectors describing the directions between two atoms, we could use some simple trigonometry to compute the angle between the vectors: \(\vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}| \cos(\theta)\). So we can calculate the angle as \(\theta = \arccos\left(\frac{\vec{a} \cdot \vec{b}}{|\vec{a}||\vec{b}|}\right)\), we just have to define our two vectors \(\vec{a}\) and \(\vec{b}\). We compute these vectors as the difference in positions of two atoms. For example, here we compute the angle H-N-H in an ammonia molecule. This is the angle between N-H\(_1\) and N-H\(_2\). In the next example, we utilize functions in mod:numpy to perform the calculations, specifically the func:numpy.arccos function, the func:numpy.dot function, and func:numpy.linalg.norm functions.
from ase.build import molecule
# ammonia
atoms = molecule('NH3')
print('atom symbol')
print('===========')
for i, atom in enumerate(atoms):
print('{0:2d} {1:3s}'.format(i, atom.symbol))
a = atoms.positions[0] - atoms.positions[1]
b = atoms.positions[0] - atoms.positions[2]
from numpy import arccos, dot, pi
from numpy.linalg import norm
theta_rad = arccos(dot(a, b) / (norm(a) * norm(b))) # in radians
print('theta = {0:1.1f} degrees'.format(theta_rad * 180./pi))
atom symbol
===========
0 N
1 H
2 H
3 H
theta = 106.3 degrees
from ase.build import molecule
from numpy import pi
# ammonia
atoms = molecule('NH3')
print('theta = {0} degrees'.format(atoms.get_angle(1, 0, 2)))
theta = 106.33462423179175 degrees
Dihedral angles#
There is support in ase for computing dihedral angles. Let us illustrate that for ethane. We will compute the dihedral angle between atoms 5, 1, 0, and 4. That is a H-C-C-H dihedral angle, and one can visually see (although not here) that these atoms have a dihedral angle of 60° (Figure ref:fig:ethane-dihedral).
# calculate an ethane dihedral angle
from ase.build import molecule
import numpy as np
atoms = molecule('C2H6')
print('atom symbol')
print('===========')
for i, atom in enumerate(atoms):
print('{0} {1}'.format(i, atom.symbol))
# compute the dihedral angle
dihedral = atoms.get_dihedral(4, 0, 1, 5)
print('')
print('Dihedral angle (4, 0, 1, 5) = {0:1.0f} degrees'.format(dihedral))
atom symbol
===========
0 C
1 C
2 H
3 H
4 H
5 H
6 H
7 H
Dihedral angle (4, 0, 1, 5) = 60 degrees
In this section we covered properties that require simple calculations, but not DFT calculations, to compute.