3.7. Molecule¶

In this chapter the PLAMS module responsible for handling molecular geometries is presented. Information about atomic coordinates can be read from (or written to) files of various types: xyz, pdb, mol or mol2. PLAMS not only extracts relevant data from those files, but also tries to “understand” the structure of the underlying molecule in terms of atoms and bonds between them, allowing you to perform a variety of simple operations like, for example, moving or rotating some parts of the molecule, splitting it into multiple parts, merging two molecules etc. Classes defined in this module are Molecule, Atom and Bond. They interact with each other to provide a basic set of functionalities for geometry handling.

3.7.1. Molecule¶

class Molecule(filename=None, inputformat=None, geometry=1, **other)[source]¶

A class representing basic molecule object.

An instance of this class has the following attributes:

atoms – a list of Atom objects that belong to this molecule
bonds – a list of Bond objects between atoms listed in atoms
lattice – a list of lattice vectors, in case of periodic structures
properties – a Settings instance storing all other information about this molecule

Note

Each Atom in atoms list and each Bond in bonds list has a reference to the parent molecule. Moreover, each atom stores the list of bonds it’s a part of and each bond stores references to atoms it bonds. That creates a complex net of references between objects that are part of a molecule. Consistency of this data is crucial for proper functioning of many methods. Because of that it is advised not to modify contents of atoms and bonds by hand. When you need to alter your molecule, methods add_atom(), delete_atom(), add_bond() and delete_bond() can be used to ensure that all these references are updated properly.

Creating a Molecule object for your calculation can be done in two ways. You can start with an empty molecule and manually add all atoms (and bonds, if needed):

>>> mol = Molecule()
>>> mol.add_atom(Atom(atnum=1, coords=(0,0,0)))
>>> mol.add_atom(Atom(atnum=1, coords=(d,0,0)))

This approach can be useful for building small molecules, especially if you wish to parametrize some of atomic coordinates (like in Simple example), but in general it’s not very practical. Usually one wants to import atomic coordinates from some external file:

>>> mol = Molecule('xyz/Benzene.xyz')

Constructor of a Molecule object accepts three arguments that can be used to supply this information from a file in your filesystem. filename should be a string with a path (absolute or relative) to such a file. inputformat describes the format of the file. Currently, the following formats are supported: xyz, mol, mol2 and pdb. If inputformat argument is not supplied, PLAMS will try to deduce it by examining the extension of the provided file, so in most of cases it is not needed to use inputformat, if only the file has the proper extension. Some formats (xyz and pdb) allow to store more than one geometry of a particular molecule within a single file. In such cases geometry argument can be used to indicate which (in order of appearance in the file) geometry to import. other keyword arguments passed to the constructor are used to populate properties Settings.

If a Molecule is initialized from an external file, the path to this file (filename argument) is stored in properties.source. The base name of the file without extension is kept in properties.name.

It is also possible to write a molecule to a file in one of the formats mentioned above. See write() for details.

lattice attribute is used to store information about lattice vectors in case of periodic structures. Some job types (BANDJob, DFTBJob) will automatically use that data while constructing input files. lattice should be a list of up to 3 vectors (for different types of periodicity: chain, slab or bulk), each of which needs to be a list or a tuple of 3 numbers.

Lattice vectors can be directly read and written to xyz files using the following convention (please mind the fact that this is an unofficial extension to the XYZ format):

3

    H      0.000000      0.765440     -0.008360
    O      0.000000      0.000000      0.593720
    H      0.000000     -0.765440     -0.008360
VEC1       3.000000      0.000000      0.000000
VEC2       0.000000      3.000000      0.000000
VEC3       0.000000      0.000000      3.000000

For 1D (2D) periodicity please supply only VEC1 (VEC1 and VEC2). Writing lattice vectors to xyz files can be disabled by simply reseting the lattice attribute:

>>> mol.lattice = []

The detailed description of available methods is presented below. Many of these methods require passing atoms belonging to the molecule as arguments. It can by done by using a reference to an Atom object present it atoms list, but not by passing a number of an atom (its position within atoms list). Unlike some other tools, PLAMS does not use integer numbers as primary identifiers of atoms. It is done to prevent problems when atoms within a molecule are reordered or some atoms are deleted. References to Atom or Bond objects can be obtained directly from atoms or bonds lists, or with dictionary-like bracket notation:

>>> mol = Molecule('xyz/Ammonia.xyz')
>>> mol.guess_bonds()
>>> print(mol)
  Atoms:
    1         H      0.942179      0.000000     -0.017370
    2         H     -0.471089      0.815951     -0.017370
    3         N      0.000000      0.000000      0.383210
    4         H     -0.471089     -0.815951     -0.017370
  Bonds:
   (1)--1.0--(3)
   (2)--1.0--(3)
   (3)--1.0--(4)
>>> at = mol[1]
>>> print(at)
         H      0.942179      0.000000     -0.017370
>>> b = mol[(1,3)]
>>> print(b)
(         H      0.942179      0.000000     -0.017370 )--1.0--(         N      0.000000      0.000000      0.383210 )
>>> b = mol[(1,4)]
>>> print(b)
None

Note

Numbering of atoms within a molecule starts with 1.

However, if you feel more familiar with identifying atoms by natural numbers, you can use set_atoms_id() to equip each atom of the molecule with id attribute equal to atom’s position within atoms list. This method can also be helpful to track changes in your molecule during tasks that can reorder atoms.

copy(atoms=None)[source]¶

Return a copy of this molecule. New molecule has atoms, bonds and all other components distinct from original molecule (it is so called “deep copy”).

By default the entire molecule is copied. It is also possible to copy only some part of the molecule, indicated by atoms argument. It should be a list of atoms that belong to this molecule. Only these atoms, together with any bonds between them, are copied and included in the returned molecule.

add_atom(atom, adjacent=None)[source]¶

Add new atom to this molecule.

atom should be an Atom instance that does not belong to the molecule. Bonds between the new atom and other atoms of the molecule can be automatically added based on adjacent argument. It should be a list describing atoms of the molecule that the new atom is connected to. Each element of adjacent list can either be a pair (Atom, order) to indicate new bond’s order (use Bond.AR for aromatic bonds) or an Atom instance (a single bond is inserted in this case).

Example:

>>> mol = Molecule() #create an empty molecule
>>> h1 = Atom(symbol='H', coords=(1.0, 0.0, 0.0))
>>> h2 = Atom(symbol='H', coords=(-1.0, 0.0, 0.0))
>>> o = Atom(symbol='O', coords=(0.0, 1.0, 0.0))
>>> mol.add_atom(h1)
>>> mol.add_atom(h2)
>>> mol.add_atom(o)
>>> mol.add_atom(Atom(symbol='C', coords=(0.0, 0.0, 0.0)), adjacent=[h1, h2, (o,2)])

delete_atom(atom)[source]¶

Delete atom from this molecule.

atom should be an Atom instance that belongs to the molecule. All bonds containing this atom are removed too.

Examples:

>>> #delete all hydrogens
>>> mol = Molecule('protein.pdb')
>>> hydrogens = [atom for atom in mol if atom.atnum == 1]
>>> for i in hydrogens: mol.delete_atom(i)

>>> #delete first two atoms
>>> mol = Molecule('geom.xyz')
>>> mol.delete_atom(mol[1])
>>> mol.delete_atom(mol[1]) #since the second atom of original molecule is now the first

add_bond(arg1, arg2=None, order=1)[source]¶

Add new bond to this molecule.

This method can be used in two different ways. You can call it with just one argument being a Bond instance (other arguments are then ignored):

>>> b = Bond(mol[2], mol[4], order=Bond.AR) #create aromatic bond between 2nd and 4th atom
>>> mol.add_bond(b)

Other way is to pass two atoms (and possibly bond order) and new Bond object will be created automatically:

>>> mol.add_bond(mol[2], mol[4], order=Bond.AR)

In both cases atoms that are to be bond have to belong to the molecule, otherwise an exception is raised.

delete_bond(arg1, arg2=None)[source]¶

Delete bond from this molecule

Just like add_bond(), this method accepts either a single argument that is a Bond instance, or two arguments being instances of Atom. In both cases objects used as arguments have to belong to the molecule.

delete_all_bonds()[source]¶: Delete all bonds from the molecule.

find_bond(atom1, atom2)[source]¶: Find and return a bond between atom1 and atom2. Both atoms have to belong to the molecule. If a bond between chosen atoms does not exist, None is returned.

set_atoms_id()[source]¶: Equip each atom of this molecule with id attribute equal to its position within atoms list.

unset_atoms_id()[source]¶: Delete id attributes of all atoms.

neighbors(atom)[source]¶

Return a list of neighbors of atom within this molecule.

atom has to belong to the molecule. Returned list follows the same order as bonds list of atom.

separate()[source]¶

Separate this molecule into connected components.

Returned is a list of new Molecule objects (all atoms and bonds are disjoint with original molecule). Each element of this list is identical to one connected component of the base molecule. A connected component is a subset of atoms such that there exists a path (along one or more bonds) between any two atoms.

Example:

>>> mol = Molecule('/xyz_dimers/NH3-H2O.xyz')
>>> mol.guess_bonds()
>>> print(mol)
  Atoms:
    1         N     -1.395591     -0.021564      0.000037
    2         H     -1.629811      0.961096     -0.106224
    3         H     -1.862767     -0.512544     -0.755974
    4         H     -1.833547     -0.330770      0.862307
    5         O      1.568501      0.105892      0.000005
    6         H      0.606736     -0.033962     -0.000628
    7         H      1.940519     -0.780005      0.000222
  Bonds:
   (5)--1.0--(7)
   (5)--1.0--(6)
   (1)--1.0--(3)
   (1)--1.0--(4)
   (1)--1.0--(2)
>>> x = mol.separate()
>>> for i in x: print(i)
  Atoms:
    1         N     -1.395591     -0.021564      0.000037
    2         H     -1.629811      0.961096     -0.106224
    3         H     -1.862767     -0.512544     -0.755974
    4         H     -1.833547     -0.330770      0.862307
  Bonds:
   (1)--1.0--(3)
   (1)--1.0--(4)
   (1)--1.0--(2)

  Atoms:
    1         O      1.568501      0.105892      0.000005
    2         H      0.606736     -0.033962     -0.000628
    3         H      1.940519     -0.780005      0.000222
  Bonds:
   (1)--1.0--(3)
   (1)--1.0--(2)

guess_bonds()[source]¶

Try to guess bonds in the molecule based on types and positions of atoms.

All previously existing bonds are removed. New bonds are generated based on interatomic distances and information about maximal number of bonds for each atom type (connectors property, taken from PeriodicTable).

The problem of finding molecular bonds for a given set of atoms in space does not have a general solution, especially considering the fact the chemical bond is itself not a precisely defined concept. For every method, no matter how sophisticated, there will always be corner cases for which the method produces disputable results. Moreover, depending on the context (area of application) the desired solution for a particular geometry may vary. Please do not treat this method as an oracle always providing proper solution. Algorithm used here gives very good results for geometries that are not very far from optimal geometry, especially consisting of lighter atoms. All kinds of organic molecules, including aromatic ones, usually work very well. Problematic results can emerge for transition metal complexes, transition states, incomplete molecules etc.

The algorithm used scales as n log n where n is the number of atoms.

Warning

This method works reliably only for geometries representing complete molecules. If some atoms are missing (for example, a protein without hydrogens) the resulting set of bonds would usually contain more bonds or bonds with higher order than expected.

translate(vector, unit='angstrom')[source]¶

Move this molecule in space by vector, expressed in unit.

vector should be an iterable container of length 3 (usually tuple, list or numpy array). unit describes unit of values stored in vector.

rotate_lattice(matrix)[source]¶

Rotate only lattice vectors of this molecule with given rotation matrix.

matrix should be a container with 9 numerical values. It can be a list (tuple, numpy array etc.) listing matrix elements row-wise, either flat ([1,2,3,4,5,6,7,8,9]) or in two-level fashion ([[1,2,3],[4,5,6],[7,8,9]]).

Note

This method does not check if supplied matrix is a proper rotation matrix.

rotate(matrix, lattice=False)[source]¶

Rotate this molecule with given rotation matrix. If lattice is True, rotate also the lattice vectors.

matrix should be a container with 9 numerical values. It can be a list (tuple, numpy array etc.) listing matrix elements row-wise, either flat ([1,2,3,4,5,6,7,8,9]) or in two-level fashion ([[1,2,3],[4,5,6],[7,8,9]]).

Note

This method does not check if supplied matrix is a proper rotation matrix.

align_lattice(convention='x', zero=1e-10)[source]¶

Rotate this molecule in such a way that lattice vectors are aligned with coordinate system.

This method is meant to be used with periodic systems only. Using it on a Molecule instance with empty lattice attribute has no effect.

Possible values of convention argument are:

x (default) – first lattice vector aligned with X axis. Second vector (if present) aligned with XY plane.
z (convention used by ReaxFF) – second lattice vector (if present) aligned with YZ plane. Third vector (if present) aligned with Z axis.

zero argument can be used to specify the numerical tolerance for zero (used to determine if some vector is already aligned with a particular axis or plane).

The returned value is True if any rotation happened, False otherwise.

rotate_bond(bond, atom, angle, unit='radian')[source]¶

Rotate given bond by an angle expressed in unit.

bond should be chosen in such a way, that it divides the molecule into two parts (using a bond being part of a ring results in an error). atom has to belong to bond and is used to pick which “half” of the molecule is rotated. Positive angle denotes counterclockwise rotation (looking along the bond, from the stationary part of the molecule).

closest_atom(point, unit='angstrom')[source]¶

Return the atom of this molecule that is the closest one to some point in space.

point should be an iterable container of length 3 (for example: tuple, Atom, list, numpy array). unit describes unit of values stored in point.

distance_to_point(point, unit='angstrom', result_unit='angstrom')[source]¶

Calculate the distance between this molecule and some point in space (distance between point and closest_atom()).

point should be an iterable container of length 3 (for example: tuple, Atom, list, numpy array). unit describes unit of values stored in point. Returned value is expressed in result_unit.

distance_to_mol(other, result_unit='angstrom', return_atoms=False)[source]¶

Calculate the distance between this molecule and some other molecule.

The distance is measured as the smallest distance between a pair of atoms, one belonging to each of the molecules. Returned distance is expressed in result_unit.

If return_atoms is False, only a single number is returned. If return_atoms is True, this method returns a tuple (distance, atom1, atom2) where atom1 and atom2 are atoms fulfilling the minimal distance, with atom1 belonging to this molecule and atom2 to other.

wrap(self, length, angle=2*pi, length_unit='angstrom', angle_unit='radian')[source]¶

Transform the molecule wrapping its x-axis around z-axis. This method is useful for building nanotubes or molecular wedding rings.

Atomic coordinates are transformed in the following way:

z coordinates remain untouched
x axis gets wrapped around the circle centered in the origin of new coordinate system. Each segment of x axis of length length ends up as an arc of a circle subtended by an angle angle. The radius of this circle is R = length/angle.
part of the plane between the x axis and the line y=R is transformed into the interior of the circle, with line y=R being squashed into a single point - the center of the circle.
part of the plane above line y=R is dropped
part of the plane below x axis is transformed into outside of the circle
transformation is done in such a way that distances along y axis are preserved

Before:

After:

get_center_of_mass(unit='angstrom')[source]¶: Return the center of mass of this molecule (as a tuple). Returned coordinates are expressed in unit.

get_mass()[source]¶: Return mass of the molecule, expressed in atomic units.

get_formula()[source]¶

Calculate the molecular formula for this molecule.

Returned value is a single string. It contains simple molecular formula (it only includes atom types and total number of atoms of each type).

apply_strain(strain)[source]¶

Apply a strain deformation to a periodic system.

This method can be used only for periodic systems (the ones with non-empty lattice attribute). strain should be a container with n*n numerical values, where n is the size of the lattice. It can be a list (tuple, numpy array etc.) listing matrix elements row-wise, either flat ([1,2,3,4,5,6,7,8,9]) or in two-level fashion ([[1,2,3],[4,5,6],[7,8,9]]).

__len__()[source]¶: Length of a molecule is the number of atoms.

__str__()[source]¶

Return string representation of this molecule.

Information about atoms are printed in xyz format fashion – each atom in a separate, enumerated line. Then, if the molecule contains any bonds, they are printed. Each bond is printed in a separate line, with information about both atoms and bond order. Example:

Atoms:
  1         N       0.00000       0.00000       0.38321
  2         H       0.94218       0.00000      -0.01737
  3         H      -0.47109       0.81595      -0.01737
  4         H      -0.47109      -0.81595      -0.01737
Bonds:
  (1)----1----(2)
  (1)----1----(3)
  (1)----1----(4)

__iter__()[source]¶: Iterate over atoms.

__getitem__(key)[source]¶

Bracket notation can be used to access atoms or bonds directly.

If key is a single int (mymol[i]), return i-th atom of this molecule. If key is a pair of ints (mymol[(i,j)]), return bond between i-th and j-th atom (None if such a bond does not exist).

This method is read only (things like mymol[3] = Atom(...) are forbidden). Numbering of atoms withing a molecule starts with 1.

__add__(other)[source]¶

Create a new molecule that is a sum of this molecule and other:

>>> newmol = mol1 + mol2

The new molecule has atoms, bonds and all other elements distinct from both components. properties of newmol are properties of mol1 soft_updated with properties of mol2.

__iadd__(other)[source]¶

Add other molecule to this one:

>>> protein += water

All atoms and bonds present in other are copied and copies are added to this molecule. properties of this molecule are soft_updated with properties of other.

read(filename, inputformat=None, frame=1)[source]¶

Read molecular coordinates from file.

filename should be a string with a path to the file. If inputformat is not None, it should be one of supported formats (keys occurring in class attribute _readformat). Otherwise, format of the file is deduced from file’s extension (for files without extension xyz format is assumed).

If chosen format allows multiple geometries in a single file, frame can be used to pick one of them.

write(filename, outputformat=None)[source]¶

Write molecular coordinates to a file.

filename should be a string with a path to the file. If outputformat is not None, it should be one of supported formats (keys occurring in class attribute _writeformat). Otherwise, format of the file is deduced from file’s extension (for files without extension xyz format is assumed).

as_dict()[source]¶

Store all the information about this Molecule in a dictionary.

Returned dictionary is, in principle, identical to self.__dict__ of the current instance, apart from the fact that all Atom and Bond instances in atoms and bonds lists are replaced with dictionaries storing corresponing information.

This method is a counterpart of from_dict().

classmethod from_dict(dictionary)[source]¶

Generate a new Molecule instance based on the information stored in dictionary.

This method is a counterpart of as_dict().

3.7.2. Atom¶

class Atom(atnum=0, symbol=None, coords=None, unit='angstrom', bonds=None, mol=None, **other)[source]¶

A class representing a single atom in three dimensional space.

An instance of this class has the following attributes:

atnum – atomic number (zero for “dummy atoms”)
coords – tuple of length 3 storing spatial coordinates
bonds – list of bonds (see Bond) this atom is a part of
mol – a Molecule this atom belongs to
properties – a Settings instance storing all other information about this atom (initially it is populated with **other keyword arguments passed to the constructor)

All the above attributes can be accessed either directly or using one of the following properties:

x, y, z – allow to read or modify each coordinate separately
symbol – allows to read or write atomic symbol directly. Atomic symbol is not stored as an attribute, instead of that atomic number (atnum) indicates the type of atom. In fact, symbol this is just a wrapper around atnum that uses PeriodicTable as a translator:
```
>>> a = Atom(atnum=8)
>>> print(a.symbol)
O
>>> a.symbol = 'Ca'
>>> print(a.atnum)
20
```
mass – atomic mass, obtained from PeriodicTable, read only
radius – atomic radius, obtained from PeriodicTable, read only
connectors – number of connectors, obtained from PeriodicTable, read only

Note

When creating a new atom, its type can be chosen either by setting an atomic number or a symbol (atnum and symbol constructor arguments). Symbol takes precedence – if it is supplied, atnum argument is ignored.

Values stored in coords tuple do not necessarily have to be numeric, you can also store any string there. This might come handy for programs that allow parametrization of coordinates in the input file (to enforce some geometry constraints for example):

>>> a = Atom(symbol='C', coords=(1,2,3))
>>> print(a)
         C       1.00000       2.00000       3.00000
>>> a.y = 'param1'
>>> print(a)
         C       1.00000        param1       3.00000

However, non-numerical coordinates cannot be used together with some methods (for example distance_to() or translate()). Trying to do this will raise an exception.

Internally, atomic coordinates are always expressed in angstroms. Most of methods that read or modify atomic coordinates accept keyword argument unit allowing to choose unit in which results and/or arguments are expressed (see Units for details). Throughout the entire code angstrom is the default length unit. If you don’t specify unit parameter in any place of your script, all automatic unit handling described above boils down to occasional multiplication/division by 1.0.

str(symbol=True, suffix='', suffix_dict={}, unit='angstrom', space=14, decimal=6)[source]¶

Return a string representation of this atom.

Returned string is a single line (no newline characters) that always contains atomic coordinates (and maybe more). Each atomic coordinate is printed using space characters, with decimal characters reserved for decimal digits. Coordinates values are expressed in unit.

If symbol is True, atomic symbol is added at the beginning of the line. If symbol is a string, this exact string is printed there.

suffix is an arbitrary string that is appended at the end of returned line. It can contain identifiers in curly brackets (like for example f={fragment}) that will be replaced by values of corresponding keys from suffix_dict dictionary. See Format String Syntax for details.

Example:

>>> a = Atom(atnum=6, coords=(1,1.5,2))
>>> print(a.str())
         C      1.000000      1.500000      2.000000
>>> print(a.str(unit='bohr'))
         C      1.889726      2.834589      3.779452
>>> print(a.str(symbol=False))
      1.000000      1.500000      2.000000
>>> print(a.str(symbol='C2.13'))
     C2.13      1.000000      1.500000      2.000000
>>> print(a.str(suffix='protein1'))
         C      1.000000      1.500000      2.000000 protein1
>>> a.properties.info = 'membrane'
>>> print(a.str(suffix='subsystem={info}', suffix_dict=a.properties))
         C      1.000000      1.500000      2.000000 subsystem=membrane

__str__()[source]¶: Return a string representation of this atom. Simplified version of str() to work as a magic method.

__iter__()[source]¶: Iteration through atom yields coordinates. Thanks to that instances of Atom can be passed to any method requiring point or vector as an argument.

translate(vector, unit='angstrom')[source]¶

Move this atom in space by vector, expressed in unit.

vector should be an iterable container of length 3 (usually tuple, list or numpy array). unit describes unit of values stored in vector.