Chemical System¶

Overview¶

The ChemicalSystem class serves as a versatile representation of a chemical system. It’s designed to handle various types of chemical structures, such as molecules, surfaces, or crystals.

Here’s how you can initialize a ChemicalSystem object using a System Block string:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem

# Initialize a Chemical System from a 'System Block' string
mol = ChemicalSystem(
"""
System
     Atoms
          O 0.0 0.0 0.0
          H 1.0 0.0 0.0
          H 0.0 1.0 0.0
     End
End
""")

# Guess the bonds in the molecule based on the atomic coordinates
mol.guess_bonds()

# Translate the molecule so that the origin coincide with it's center of mass
mol.translate(-mol.center_of_mass())

# Print the molecule in 'System Block' format
print(mol)

Reading and writing¶

You can create or serialize a ChemicalSystem object using various file formats. Among these, System Block is one of the most significant, offering a versatile text-based way to describe your chemical system. For more information on the syntax and options for the System Block, see the AMS System Block documentation.

When you convert a ChemicalSystem object into a string (either by explicitly calling str(my_chemical_system) or by using a print statement) the output will be in the System Block format

>>> from scm.libbase import UnifiedChemicalSystem as ChemicalSystem

>>> # Read it from a text file in the 'System Block' format:
>>> my_chemical_system = ChemicalSystem.from_in(filename="water.in")

>>> print(my_chemical_system)

System
    Atoms
        O    0.0000000000000000  0.0000000000000000  0.0000000000000000
        H    1.0000000000000000  0.0000000000000000  0.0000000000000000
        H    0.0000000000000000  1.0000000000000000  0.0000000000000000
    End
End

>>> # Write it to a file in the 'System Block' format:
>>> my_chemical_system.write_in(filename="another_water.in")

Note on Serialization

When you read or write a ChemicalSystem using either the System Block format or a kf file, the object is perfectly serialized. In other words, writing the object to a kf file and reading it back will result in an identical ChemicalSystem.

However, be cautious when using the xyz format as it doesn’t offer lossless serialization. Writing and reading back using this format may result in the loss of certain information, such as bonds between atoms.

Available Methods for Serialization

Below are the methods you can use for reading and writing (also known as serializing and deserializing) a ChemicalSystem:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

classmethod from_in(filename: str, name: str = '') → UnifiedChemicalSystem¶: Constructs and returns a new ChemicalSystem from a (possibly named) System block in an AMS input file.

classmethod from_kf(filename: str, section: str = 'Molecule') → UnifiedChemicalSystem¶
classmethod from_kf(kf: libbase.KFFile, section: str = 'Molecule') → UnifiedChemicalSystem: Constructs and returns a new ChemicalSystem from a section on a KF file.

classmethod from_xyz(filename: str) → UnifiedChemicalSystem¶: Constructs and returns a new ChemicalSystem from an extended XYZ file.

classmethod from_input(input_file: libbase.InputFile, prefix: str) → UnifiedChemicalSystem¶: Constructs and returns a new ChemicalSystem from an InputFile instance, given a prefix, e.g. ‘System[1]%’.

write_in(filename: str) → None¶: Writes a ChemicalSystem to file as an AMS System block.

write_kf(filename: str, section: str = 'Molecule', omode: libbase.KFFile.OpenMode = KFFile.OpenMode.Any) → None¶
write_kf(kf: libbase.KFFile, section: str = 'Molecule') → None: Writes a ChemicalSystem to a section on a KF file.

write_xyz(filename: str, extended_xyz_format: bool = True) → None¶

Writes a ChemicalSystem to an XYZ file.

By default the file is written in the AMS extended XYZ file format, which includes lattice vectors, atomic properties, regions and some system properties such as the total charge. In order to write a plain standard XYZ files without any of these extensions, set extended_xyz_format to False.

Atoms¶

The ChemicalSystem class maintains an array of instances of the Atom class (UnifiedAtom) to represent the atoms in the system.

Adding and Modifying Atoms

The following code provides examples for adding atoms to a ChemicalSystem and modifying their properties.

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem
from scm.libbase import UnifiedAtom as Atom

mol = ChemicalSystem()

# Three different ways of adding atoms
mol.add_atom("O", coords=[0, 0, 0])
mol.add_atom(Atom("H"), coords=[2, 0, 0])
mol.add_atom(1, coords=[0, 2, 0], unit="Angstrom")

# Change the coordinates of the Oxygen atom (units: Bohr)
mol.atoms[0].coords = [1, 2, 3]

# This is equivalent to the line above:
mol.coords[:, 0] = [1, 2, 3]

# Change the first hydrogen atom to deuterium by adjusting the mass (units: Dalton):
mol.atoms[1].mass = 2.014

Available Methods

Here are methods to interact with the atoms in your ChemicalSystem:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property atoms¶: A list of all Atom instances that are part of this ChemicalSystem.

property coords¶: The coordinates of the atoms in bohr.

contains_atom(atom: libbase.UnifiedAtom) → bool¶: Checks if an Atom instance is part of a ChemicalSystem.

atom_index(atom: libbase.UnifiedAtom) → int¶: Given an Atom instance, returns its index in ChemicalSystem.atoms.

add_atom(atom: libbase.UnifiedAtom) → None¶
add_atom(atom: libbase.UnifiedAtom, coords: ArrayLike, unit: str = 'bohr') → None
add_atom(Z: int, coords: ArrayLike, unit: str = 'bohr') → None
add_atom(element: libbase.UnifiedElement, coords: ArrayLike, unit: str = 'bohr') → None
add_atom(symbol: UnifiedAtom.T_AtomSymbol, coords: ArrayLike, unit: str = 'bohr') → None: Add a new Atom to the ChemicalSystem.

remove_atom(atom_index: int) → None¶: Removes a single atom from the system, given its index.

remove_atoms(atom_indices: ArrayLike) → None¶: Removes multiple atoms from the system, given their atom indices.

set_atom(atom_index: int, atom: libbase.UnifiedAtom) → None¶

Safe assignment of an Atom instance to atoms[atom_index].

Is essentially the same as self.atoms[atom_index] = atom, but makes sure that all the atomic properties carried by the Atom instance are enabled on the target ChemicalSystem: the direct assignment will throw a ChemicalSystemError if the needed atomic properties are not enabled, while calling this method will enable them automatically and will not throw.

If Atom is part of a ChemicalSystem, its will overwrite whatever coordinates atoms[atom_index] had before. If Atom is not part of a ChemicalSystem, the coordinates of atoms[atom_index] will remain unchanged.

Methods for combining two ChemicalSystems into one:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

add_other(other: UnifiedChemicalSystem) → None¶

Merges another ChemicalSystem into this one.

The number of atoms of self increases by the number of atoms of other. All atomic coordinates, properties and bonds between atoms will be kept. The total charge of other will be added to the total charge of self.

The regions of each atom do not change in the process. Regions with the same name in self and other are merged.

Systems with a lattice can only be merged with systems having the same lattice or no lattice at all. If only one side has a lattice, that side determines the lattice of the result.

If both sides have a compatible lattice and bonds, merging them is only supported if either both or none of the two have the lattice displacements of the bonds set. If one side has lattice displacements, and the other does not, a ChemicalSystemError is raised.

__iadd__(arg0: UnifiedChemicalSystem) → UnifiedChemicalSystem¶

Merges another ChemicalSystem into this one.

Note that lhs += rhs is just the operator version of lhs.add_other(rhs). See ChemicalSystem.add_other for details about merging systems.

__add__(arg0: UnifiedChemicalSystem) → UnifiedChemicalSystem¶

Creates a new ChemicalSystem by merging two others.

Note that C = A + B is equivalent to C = copy(A); C.add_other(B). See ChemicalSystem.add_other for details about merging systems.

Splitting of ChemicalSystems into parts:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

extract_atoms(atom_indices: ArrayLike) → UnifiedChemicalSystem¶

Returns a new system build from a subset of atoms.

Bonds within the subset will be preserved, but bonds to atoms not extracted will disappear. The returned system will have the same lattice as the original and a total charge of zero. Atomic properties and regions of the extracted atoms are preserved, but the returned system does not have any selected atoms.

split(part_indices: Sequence[numpy.int64]) → List[UnifiedChemicalSystem]¶

Splits the system into parts and returns a list of these parts as separate systems.

Accepts a num_atoms long sequence, assigning the atoms of the system to the parts. The length of the returned list of parts is max(part_indices)+1.

Example: for a 6 atom system and a part_indices of [0, 0, 0, 1, 1, 1] a list of two systems will be returned. The system will contain the first three atoms of the original, and the second system the other three atoms.

The returned systems will have the same lattice as the original and a total charge of zero. Atomic properties and regions are preserved, but the returned systems do not have any selected atoms.

The UnifiedAtom class

The UnifiedAtom class contains all the relevant data for an individual atom:

class UnifiedAtom

class UnifiedAtom(Z: int)

class UnifiedAtom(element: UnifiedElement)

class UnifiedAtom(symbol: T_AtomSymbol)

Class representing a single atom and its properties.

property Z¶: The atomic number of the atom.

property coords¶: The coordinates of the atom in bohr.

property symbol¶: The element symbol, e.g. ‘C’ or ‘Au’.

in_chemicalsystem() → bool¶: Check if this Atom instance is part of a ChemicalSystem.

Molecular properties¶

These are properties of the Chemical System:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property charge¶: The total charge of the system in atomic units (i.e. in units of elementary charge).

total_mass() → float¶: Total mass in atomic mass units (dalton).

formula() → str¶

Chemical formula in Hill notation (e.g. H2O, CH4):

If Carbon is present, the order is: C, H, alphabetical.
If Carbon is not present all the elements (including Hydrogen) are listed alphabetically.
Single-letter elements come before 2-letters elements (e.g. ‘K’ comes before ‘Kr’, ‘B’ before ‘Be’)

Geometry and manipulation¶

The Chemical System provides a range of methods for manipulating the geometry and retrieving geometrical properties of your molecular system.

Here is a simple example showcasing some of the methods:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem

# Initialize a Chemical System from an XYZ file
mol = ChemicalSystem.from_xyz("some_mol.xyz")
print("Original system:")
print(mol)

# Translate the system to center its geometric center at the origin
mol.translate(-mol.geometric_center())

# Define a rotation matrix for a 90-degree rotation around the z-axis
rot_mat = [[0,-1, 0],
           [1, 0, 0],
           [0, 0, 1]]

# Apply the rotation
mol.rotate(rot_mat)
print("Roto-translated system:")
print(mol)

This will output:

Original system:
System
    Atoms
        H    0.0000000000000000  0.0000000000000000  0.0000000000000000
        F    1.0000000000000000  0.0000000000000000  0.0000000000000000
    End
End

Roto-translated system:
System
    Atoms
        H    0.0000000000000000 -0.5000000000000000  0.0000000000000000
        F    0.0000000000000000  0.5000000000000000  0.0000000000000000
    End
End

Available Methods

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

translate(shift: ArrayLike) → None¶: Translates all atoms in the ChemicalSystem by a vector.

rotate(rot_mat: ArrayLike) → None¶: Rotate the system according to the rotation matrix.

align_to(other: UnifiedChemicalSystem) → None¶

Translate and rotate the system to maximally align it with ‘other’. It will first translate the system to the center of mass of ‘other’, then it will use the Kabsch algorithm to rotate the system in order to minimize the RMSD.

Notes:

The two UnifiedChemicalSystems must have the same number of atoms.
The atoms in the two chemical system must be in the same order.

classmethod rotation_matrix_minimizing_rmsd(a: UnifiedChemicalSystem, b: UnifiedChemicalSystem) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Given two chemical systems, returns the rotation matrix that minimizes the RMSD.

classmethod rmsd(a: UnifiedChemicalSystem, b: UnifiedChemicalSystem, align: bool) → float¶

Computes the RMSD between two systems. ‘align’: whether or not the systems should be roto-translated as to minimize the RMSD.

Notes:

The two UnifiedChemicalSystems must have the same number of atoms.
The atoms in the two chemical system must be in the same order.

geometric_center() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Position of the geometric center in bohr.

center_of_mass() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Position of the center of mass in bohr.

inertia_tensor() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶

Returns the system’s inertia tensor as a 3x3 matrix in atomic units (amu*bohr^2).

The inertia tensor is not defined for periodic systems. Calling this method on a periodic systems throws a ChemicalSystemError.

moments_of_inertia() → Tuple[numpy.ndarray[Any, numpy.dtype[numpy.float64]], numpy.ndarray[Any, numpy.dtype[numpy.float64]]]¶

Calculates the system’s moments of inertia and the corresponding principal axes.

Returns a tuple:

a 3 component np.array of the inertial moments in atomic units (amu*bohr^2) sorted by magnitude
a 3x3 np.array with with the corresponding axes as the column vectors

The moments of inertia and the principal axes are not defined for periodic systems. Calling this method on a periodic systems throws a ChemicalSystemError.

is_linear() → bool¶

Checks if a ChemicalSystem is linear.

A molecule is considered linear if all of its atoms are located on a line. A 1D periodic system is considered linear if all atoms are located on a line parallel to the only lattice vector. Slab and bulk periodic systems can never be linear.

perturb_coordinates(noise_level: float, unit: str = 'bohr') → None¶

Perturb the atomic coordinates by adding random numbers between [-noise_level,noise_level] to each Cartesian component.

This can be useful if you want to break the symmetry of your system (e.g. for a geometry optimization).

Bonds¶

The Chemical System contains a bonds property, which is an instance of the UnifiedBonds class. This contains the (possibly empty) bonding information on the system.

Creating and adding bonds

There are several ways of getting bonds into your chemical system. The following snippet will show some of the most common approaches:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem
from scm.libbase import UnifiedBond as Bond

# Here we initialize a Chemical System from an XYZ file.
# The chemical system will not have any bonds at the moment, because
# bonds are not defined in the .xyz file.
water = ChemicalSystem.from_xyz("water.xyz")

# You can let the bond guessing algorithm simply guess the bonds:
water.guess_bonds()


# Or you can add them manually.
# First we clear the bonds we just guessed:
water.bonds.clear_bonds()

# and add bonds manually:
water.bonds.add_bond(0, 1, Bond(1.0)) # From the oxygen to the first Hydrogen
water.bonds.add_bond(0, 2, Bond(1.0)) # From the oxygen to the second Hydrogen


# You can also explicitly define the bonds in the System Block...
system_block = """
System
    Atoms
        O    0.0000000000000000  0.0000000000000000  0.3695041700000000
        H    0.0000000000000000  0.7838367199999999 -0.1847520900000000
        H    0.0000000000000000 -0.7838367199999999 -0.1847520900000000
    End
    BondOrders
        1 2 1.00
        1 3 1.00
    End
End
"""

# ...and then load the system:
water = ChemicalSystem(system_block)

# This is how you can iterate over the bonds and get the indices of the bonded atoms:
for i, j, bond in water.bonds:
    print(f"Atoms: {i}-{j} Bond order: {bond.order}")

Note on atom indexing: as you can also see from the example above, atom indexing from within python and in the System Block differ: In python the first atom has index 0, while in the System Block the first atom has index 1.

The Bond Class

This contains the properties of a single bond.

class UnifiedBond¶

class UnifiedBond(order: float)

class UnifiedBond(order: float, lattice_displacements: ArrayLike)

A class representing a bond between two atoms.

property lattice_displacements¶

Lattice displacements for the second atom that is part of the bond.

If atom A and B are bonded with lattice displacements [1, 0, 0], this means that atom A is bonded to the image of atom B that is displaced along the first lattice vector. Note that A > B in terms of atom indices.

property order¶: The order of the bond: 1 = single bond, 1.5 = aromatic bond, 2 = double bond, 3 = triple bond.

The Bonds Class

This contains the information all the bonds, and methods to retrieve and manipulate bonds:

class UnifiedBonds

A class representing a set of bonds between atoms in a ChemicalSystem.

add_bond(from_atom: int, to_atom: int, bond: libbase.UnifiedBond) → None¶: Adds a bond between two atoms, given their indices.

add_bonds(from_atoms: ArrayLike, to_atoms: ArrayLike, bonds: List[libbase.UnifiedBond]) → None¶
add_bonds(from_atoms: ArrayLike, to_atoms: ArrayLike, orders: ArrayLike) → None: Adds bonds between multiple atom pairs at the same time. All arguments should have the same length.

num_bonds() → int¶
num_bonds(atom: int) → int: If called with an atom index, returns the number of bonds of that atom. If called without an atom index, returns the total number of bonds between all atoms.

atoms_are_bonded(from_atom: int, to_atom: int) → bool¶: Checks whether two atoms are bonded or not.

get_bonds_for_atom(from_atom: int) → Iterator[Tuple[int, int, libbase.UnifiedBond]]¶

Returns an iterator over all bonds from a particular atom given by its index.

Unpacking the iterator yields a 3-tuple of two atom indices A and B, and a Bond instance. Note that A <= B.

remove_bond(bidx: Union[int, libbase.UnifiedBond]) → None¶: Removes a bond, given either its index in the .bonds attribute or an instance of a Bond.

remove_bonds_between_atoms(from_atom: int, to_atom: int) → None¶: Removes the bonds between two atoms, gived their indices.

remove_bonds_to_atom(atom: int) → None¶: Removes all bonds to a particular atom, given its index.

remove_bonds_to_atoms(remove_atoms: ArrayLike) → None¶: Removes all bonds for a list of atoms given by their indices.

Methods in the Chemical System

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property bonds¶: The bonds of the system.

guess_bonds() → None¶: Guesses bonds based on the atomic elements and the geometry. Keeps existing bonds.

num_molecules() → int¶

Counts the number of connected molecules that are part of this system.

What is considered a molecule is based on the bonds of the system: anything connected by any bonds is considered to be a molecule. For a fully connected system, this method returns 1. If the system has no bonding information, each atom is its own molecule and this method returns the total number of atoms in the system.

molecule_indices() → numpy.ndarray[Any, numpy.dtype[numpy.int64]]¶

Returns a num_atoms sized array mapping the atoms to connected molecules of the system.

What is considered a molecule is based on the bonds of the system: anything connected by any bonds is considered to be a molecule. For a system consisting of two water molecules (where the first three atoms make up one of the molecules), this method returns [0, 0, 0, 1, 1, 1]. If the system has no bonding information, each atoms is considered its own molecule and this method returns [0, 1, 2, ..., num_atoms].

split_into_molecules() → List[UnifiedChemicalSystem]¶

Splits the system into individual molecules based on the connectivity between atoms.

The length of the returned list is num_molecules. Which set of atoms ends up in which element of the list is determined by the indices returned by the molecule_indices method.

The returned systems will have the same lattice as the original and a total charge of zero. Atomic properties and regions are preserved, but the returned systems do not have any selected atoms.

This method is just a shortcut for: cs.split(cs.molecule_indices()).

bond_cuts_molecule(from_atom: int, to_atom: int) → bool¶
bond_cuts_molecule(from_atom: libbase.UnifiedAtom, to_atom: libbase.UnifiedAtom) → bool: Checks if removing the bonds between two atoms would cut the graph into two disjoint subgraphs.

atom_is_in_ring(atom: Union[int, libbase.UnifiedAtom]) → bool¶

Checks if an atom is part of any ring.

An atom is in a ring, if removing the bond to (at least) one of its neighbors does not cut the graph into two disjoint subgraphs.

property molgraph_edgeweight_func_t¶: A method (Callable[[int, int, UnifiedBond], float]) taking two atom indices and a Bond instance between the two atoms. It shall return the non-negative distance measurement between the two atoms. This is an (often optional) input to methods that determine the shortest path between atoms in the molecular graph.

molgraph_dijkstra(from_atidx: int, dist_func: molgraph_edgeweight_func_t, to_atidx: int = - 1) → Tuple[List[float], List[int], List[bool]]¶

General method implementing the Dijkstra algorithm on the molecular graph. The Dijktra algorithm solves the single-source shortest path problem in weighted graphs with non-negative weights.

from_atidx is the index of source atom, i.e. the atom from which we start searching for paths.
dist_func is a method (of type molgraph_edgeweight_func_t = Callable[[int, int, UnifiedBond], float]) taking two atomic indices and a Bond instance as input. It returns the weight of the graph edge, aka the distance between the two atoms. The returned distance may not be negative, but you can return float("inf") to ignore a particular edge in the search.
to_atidx marks an optional target atom and can be used to solve the single-pair shortest path problem by prematurely terminating the Dijkstra algorithm as soon as the shortest path to the target atom has been found.

Returns a 3-tuple consisting of:

path_lengths is a list of length num_atoms that upon return of the method will be populated with the lengths of the shortest paths from the source atom. Unreachable atoms, aka atoms in other molecules, will have an infinite path length.
path_through is a list of length num_atoms that upon return will for each atom contain the preceding atoms index in the shortest path from source atom. By stepping backwards through this array all shortest paths from the source atom can be recovered. If no path has been found to a particular atom, the corresponding list element will be set to -1.
visited is a list of length num_atoms marking which atoms have been visited by the Dijkstra algorithm prior to its termination. An atom counts as visited if it has been popped of the top of the priority queue to have its neighbors added to the queue. Note that in case of premature termination of the algorithm due to using the to_atidx argument, shortest paths have only been determined for all visited atoms. (Other atoms may still have a path length set, but it is not guaranteed to be the shortest path from the source.)

shortest_path_lengths_from(from_atom: Union[int, libbase.UnifiedAtom]) → List[float]¶
shortest_path_lengths_from(from_atom: Union[int, libbase.UnifiedAtom], dist_func: molgraph_edgeweight_func_t) → List[float]: Returns the lengths of all shortest paths from from_atom to all other atoms.

shortest_path_between(from_atom: int, to_atom: int, dist_func: molgraph_edgeweight_func_t) → Optional[List[int]]¶

shortest_path_between(from_atom: int, to_atom: int) → Optional[List[int]]

shortest_path_between(from_atom: libbase.UnifiedAtom, to_atom: libbase.UnifiedAtom, dist_func: molgraph_edgeweight_func_t) → Optional[List[int]]

shortest_path_between(from_atom: libbase.UnifiedAtom, to_atom: libbase.UnifiedAtom) → Optional[List[int]]

Returns the shortest path between two atoms as a list of atom indices.

Note that the list contains both from_atom and to_atom as the first and last element. The length of the path (measured in number of hops) is therefore equal to the length of the returned list minus one.

shortest_path_length_between(from_atom: int, to_atom: int, dist_func: molgraph_edgeweight_func_t) → float¶
shortest_path_length_between(from_atom: int, to_atom: int) → float
shortest_path_length_between(from_atom: libbase.UnifiedAtom, to_atom: libbase.UnifiedAtom, dist_func: molgraph_edgeweight_func_t) → float
shortest_path_length_between(from_atom: libbase.UnifiedAtom, to_atom: libbase.UnifiedAtom) → float: Returns the length of the shortest path between two atoms.

Lattice and Periodic Systems¶

The Chemical System can handle periodic system with arbitrary numbers of Periodic Boundaries Conditions (i.e. 0,1,2,3).

The lattice property of the Chemical System contains the lattice vectors information.

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property lattice¶: The lattice of the system.

has_lattice() → bool¶

Checks if the system has a lattice.

Note that the lattice attribute of a ChemicalSystem is never None. A system is considered to not have a lattice if the number of lattice vectors is 0. Using this method is the same a checking lattice.num_vectors == 0.

Here are some examples:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem
from scm.libbase import UnifiedLattice as Lattice

# Create an empty chemical system
chain = ChemicalSystem()

# Add an Helium atom in the origin
chain.add_atom("He", [0, 0, 0])

# Make it a 1D periodic chain by adding a lattice vector along X (units: Bohr)
chain.lattice = Lattice([[2],
                         [0],
                         [0]])
print(chain)

# You can also initialize it from a System Block:
graphene = ChemicalSystem("""
System
     Atoms
          C    0.0000000000000000  0.0000000000000000  0.0000000000000000
          C    1.2300000000000000  0.7101408311032394  0.0000000000000000
     End
     Lattice
          2.4600000000000000   0.0000000000000000   0.0000000000000000
          -1.2300000000000000   2.1304224933097191   0.0000000000000000
     End
End
""")

print(graphene)

These are the methods in the UnifiedLattice class.

class UnifiedLattice¶

class UnifiedLattice(vectors: ArrayLike)

A class representing the lattice of a ChemicalSystem.

cartesian_to_fractional(cartesian_coord: numpy.ndarray[Any, numpy.dtype[numpy.float64]]) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Convert Cartesian coordinates to fractional coordinates for a single point.

cartesians_to_fractionals(cartesian_coords: numpy.ndarray[Any, numpy.dtype[numpy.float64]]) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Convert Cartesian coordinates to fractional coordinates for a set of points.

fractional_to_cartesian(fractional_coord: numpy.ndarray[Any, numpy.dtype[numpy.float64]]) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Convert fractional coordinates to Cartesian coordinates for a single point.

fractionals_to_cartesians(fractional_coords: numpy.ndarray[Any, numpy.dtype[numpy.float64]]) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Convert fractional coordinates to Cartesian coordinates for a set of points.

get_angles() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Get the angles between the lattice vectors (note: always returns 3 numbers).

get_lengths() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Get the length of the lattice vectors (note: always returns 3 numbers).

get_vector_with_index(i: int) → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Get the lattice vector of a given index.

get_volume() → float¶: Get the volume of the unit cell (for 2D: area, for 1D length)

is_orthogonal() → bool¶: Is the lattice orthogonal?

property num_vectors¶: Returns the number of lattice vectors (i.e. the number of PBC).

property vectors¶: The lattice vectors in bohr.

Here are methods from the ChemicalSystem that relate to periodic systems:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

set_fractional_coordinates(frac_coords: ArrayLike) → None¶: Sets the fractional coordinates of all atoms. In non-periodic directions the plain coordinates should be passed in bohr.

get_fractional_coordinates() → numpy.ndarray[Any, numpy.dtype[numpy.float64]]¶: Returns the fractional coordinates of all atoms with respect to the cell. In non-periodic directions the plain coordinates are returned in bohr.

map_atoms(start_range: Union[float, ArrayLike]) → Tuple[bool, numpy.ndarray[Any, numpy.dtype[numpy.int64]]]¶

Maps all atoms into a unit cell from [start_range:start_range+1] in fractional coordinates.

No mapping will take place in non-periodic directions. The returned tuple signifies if any mapping has taken place and the shifts in units of the lattice vectors that was applied to map the atoms into the cell.

map_atoms_around_atom(atom: libbase.UnifiedAtom) → bool¶

Map all atoms around the chosen atom that will be in the center of the new unit cell.

No mapping will take place in non-periodic directions. The return value signifies if any mapping has taken place.

apply_strain(strain_matrix: ArrayLike) → None¶

Apply a strain deformation to a periodic system.

The atoms in the unit cell will be strained accordingly, keeping the fractional atomic coordinates constant.

The strain_matrix argument should be (or be convertable to) an NxN numpy array, where N is the number of lattice vectors of the system.

apply_strain_voigt(strain_voigt: Sequence[float]) → None¶

Apply a strain deformation to a periodic system.

The atoms in the unit cell will be strained accordingly, keeping the fractional atomic coordinates constant.

The strain_voigt argument should be passed in Voigt form:

for 3D periodic systems: [e_xx, e_yy, e_zz, gamma_yz, gamma_xz, gamma_xy]
for 2D periodic systems: [e_xx, e_yy, gamma_xy]
for 1D periodic systems: [e_xx]

with e.g. e_xy = gamma_xy/2.

perturb_lattice(noise_level: float) → None¶

Perturb the lattice vectors by applying random strain with matrix elements between [-noise_level,noise_level].

This can be useful if you want to deviate from an ideal symmetric geometry, for example if you look for a phase change due to high pressure.

supercell(supercell: ArrayLike) → None¶

Create a supercell by scaling the lattice vectors.

Copies of atoms will have the same properties as in the initial unit cell. Bonds are replicated between copies and across the new unit cell boundary.

supercell_trafo(supercell: ArrayLike) → None¶

Create a supercell by creating the matrix product of the lattice vectors and the supercell matrix.

Copies of atoms will have the same properties as in the initial unit cell. Bonds are replicated between copies and across the new unit cell boundary.

slice_thickness(ref_atom: int, top: float, bottom: float, miller: ArrayLike, translate: float = 0.0) → None¶: Create a 2D slab by slicing a 3D system by specifying the upper and lower bound of the resulting slab.

slice_layers(ref_atom: int, num_layers: int, miller: ArrayLike, translate: float = 0.0) → None¶: Create a 2D slab by slicing a 3D system by specifying the amount of layers in the resulting slab.

to_primitive_cell(precision: float = 0.1) → None¶: Convert a 2D or 3D cell to its primitive representation.

to_conventional_cell(precision: float = 0.1) → None¶: Convert a 2D or 3D cell to its conventional representation.

density(unit: str = 'dalton/bohr3') → float¶: Returns the density of the system in the specified unit. Only valid for 3D periodic systems.

set_density(target_density: float, unit: str = 'dalton/bohr3') → None¶: Applies a uniform strain to match the specified target density. Only valid for 3D periodic systems.

Regions¶

Regions are ‘groups’ of atoms within a ChemicalSystem. Regions are used for some of the AMS driver and its engines features, but they can also be useful tools for bookkeeping and manipulations (see also the AMS documentation on Regions).

Here is a simple example showing how to create and use regions:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem
import numpy as np

# Create a ChemicalSystem and define two regions (one for each water molecule)
mol = ChemicalSystem(
"""
System
     Atoms
          O -2.676678  0.759085  0.370636 region=water_1
          H -3.468900  0.415690  0.814339 region=water_1
          H -3.005004  1.433129 -0.246320 region=water_1
          O  0.039085  0.303872  1.265871 region=water_2
          H -0.874303  0.483757  0.975166 region=water_2
          H  0.293563 -0.534617  0.849654 region=water_2
     End
End
"""
)

# Translate all the atoms in the region "water_2"
for i_atom in mol.get_atoms_in_region("water_2"):
     mol.atoms[i_atom].coords += np.array([1, 0, 0])

Region names are not case sensitive and may not include certain characters. If a string is a valid name for a region can be checked with the is_valid_region_name method:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

classmethod is_valid_region_name(name: str) → bool¶

Checks if a string is a valid region name.

Valid region names do not contain any of the following characters:

, + - * / \ | $ & ^ < > ( ) [ ] { } " '

Valid region names do not have leading or trailing whitespace. Note that region names are case-insensitive, so “A” and “a” correspond to the same region.

There are various methods for querying region information:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property num_regions¶

Returns the number of regions used in the system.

Note that there is no such thing as an empty region. Every region contains at least one atom.

property region_names¶

Returns a list of the names of all regions in the system.

Note that there is no such thing as an empty region. Every region contains at least one atom.

has_region(region: str) → bool¶

Checks if a region exists in the system, given its region name.

Note that region should be a valid region name, not a region expression. If you want to check if a region expression contains any atoms, use num_atoms_in_region(...) > 0 instead.

do_regions_intersect(regionA: str, regionB: str) → bool¶: Checks if two regions or region expressions intersect, i.e. have at least one atom in common.

num_atoms_in_region(region: str) → int¶

Returns the number of atoms in a region or region expression.

Returns zero if the region does not exist at all, or if the region expression evaluates to an empty set.

num_atoms_outside_region(region: str) → int¶

Returns the number of atoms outside of a region or region expression.

Returns the total number of atoms in the system if the region does not exist at all, or if the region expression evaluates to an empty set.

get_atoms_in_region(region: str) → numpy.ndarray[Any, numpy.dtype[numpy.int64]]¶: Returns a sorted array of atom indices of all atoms in a region or region expression.

get_atoms_outside_region(region: str) → numpy.ndarray[Any, numpy.dtype[numpy.int64]]¶: Returns a sorted array of atom indices of all atoms outside of a region or region expression.

is_atom_in_region(atom: Union[int, libbase.UnifiedAtom], region: str) → bool¶: Checks if an atom is in a region or region expression.

is_atom_outside_region(atom: Union[int, libbase.UnifiedAtom], region: str) → bool¶: Checks if an atom is outside of a region or region expression.

get_regions_of_atom(atom: Union[int, libbase.UnifiedAtom]) → List[str]¶: Returns an alphabetically sorted list of the names of all regions an atom is part of.

Almost all methods for querying region information also support so called region expressions in addition to plain region names. Region expressions may include region names, the parentheses (), as well as the operators + for union, - for set difference and & for set intersection. Additionally the wildcard * means all atoms, while $ stands for the set of selected atoms.

# number of atoms that are either in region A or B (or both)
num_at_AB = mol.num_atoms_in_region("A + B")

# loop over all atoms that are not in the intersection of region A and B
for atidx in mol.get_atoms_in_region("* - (A & B)"):
     ...

The following methods can be used to change the assignment of atoms to regions:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

set_atoms_in_region(atom_indices: ArrayLike, region: str) → None¶

Creates or sets an entire region given the indices of the atoms to be part of the region.

If atom_indices is empty and the region previously existed, it will effectively be deleted.

add_atoms_to_region(atom_indices: ArrayLike, region: str) → None¶: Adds multiple atoms to a region, given their atom indices.

remove_atoms_from_region(atom_indices: ArrayLike, region: str) → None¶

Removes multiple atoms from a region, given their atom indices.

No error is raised if any of the atoms is not part of the region.

remove_region(region: str) → None¶: Removes a region. Atoms in that region will keep existing, but will not longer be part of that region.

add_atom_to_region(atom: Union[int, libbase.UnifiedAtom], region: str) → None¶: Adds an atom to a region.

remove_atom_from_region(atom: Union[int, libbase.UnifiedAtom], region: str) → None¶

Removes an atom from a region.

Throws a ChemicalSystemError if the atom was not part of the region it should be removed from.

remove_atom_from_all_regions(atom: Union[int, libbase.UnifiedAtom]) → None¶: Removes an atom from all regions.

Atom selection¶

The ChemicalSystem keeps track of a set of selected atoms. In the GUI an atom is selected by simply clicking on it, and the current selection is highlighted by a cyan outline and shading.

The following methods allow inspecting the current selection:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

num_selected_atoms() → int¶: The number of currently selected atoms.

is_atom_selected(atom: Union[int, libbase.UnifiedAtom]) → bool¶: Checks if an atom is currently selected.

get_selected_atoms() → numpy.ndarray[Any, numpy.dtype[numpy.int64]]¶: Returns an array of indices of all selected atoms.

You can also use the $ symbol to refer to the set of currently selected atoms in region expressions. The following bit of code loops over all selected atoms in the region named myregion:

for atidx in mol.get_atoms_in_region("$ & myregion"):
     ...

There are many methods to change the current selection. Almost all of them work by adding or removing atoms from the current selection, e.g. the select_atom method adds an atom to the current selection and is equivalent to clicking it in the GUI. The exception to this is the set_selected_atoms method, which completely replaces the current selection.

Note that the order in which atoms are selected is tracked for small selections. Small selections of up to 4 atoms are used in the GUI for interactive manipulation in internal coordinates using the sliders at the bottom of the molecule view. The selection order is relevant for manipulations of e.g. dihedral angles, as the dihedral angle between atoms (3,1,2,4) is different than between atoms (1,2,3,4). For large selections >10 atoms there are no use cases in which the selection order is relevant and for performance reasons get_selected_atoms always returns the indices of the selected atoms in ascending order.

Basic methods to change the current selection:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

select_atom(atom: Union[int, libbase.UnifiedAtom]) → None¶: Selects an atom, i.e. adds it to the current selection.

deselect_atom(atom: Union[int, libbase.UnifiedAtom]) → None¶: Deselects an atom, i.e. removes it from the current selection.

select_atoms(atom_indices: ArrayLike) → None¶: Selects multiple atoms at once, given their atom indices.

deselect_atoms(atom_indices: ArrayLike) → None¶: Deselects multiple atoms at once, given their atom indices.

set_selected_atoms(atom_indices: ArrayLike) → None¶: Sets the selection to the given atom indices. Any previous selection is cleared.

GUI style methods from AMSinput:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

select_all() → None¶: Selects all atoms.

deselect_all() → None¶: Deselects all atoms, or in other words: clears the current selection.

invert_selection() → None¶

Inverts the set of selected atoms.

All previously unselected atoms will become selected. All previously selected atoms will become unselected.

select_connected() → None¶: Selects all atoms bonded to the currently selected atoms.

select_molecule() → None¶

Using the bonds, selects entire molecules that include any currently selected atom.

This is equivalent to repeatedly calling select_connected until the selection stops growing.

select_region(region: str) → None¶: Selects atoms in a region or region expression.

select_atom_close_to_origin() → None¶: Selects the atom that is closest to the origin of the coordinate system.

select_within_radius(radius: float, unit: str = 'bohr') → None¶: Selects all atoms within a given radius of any of the currently selected atoms.

make_selection_cappable() → None¶

Extends the current selection but does not cross single bonds, unless they are to hydrogen atoms.

The intended use of this method is to select a suitable QM region for QM/MM calculations where one wants the QM region to be separated by single bonds from the rest of the molecule.

select_atoms_of_same_type() → None¶: Selects all atoms whose element is the same as of a currently selected atom.

Methods taking predicate functions:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

select_atoms_if(pred: Callable[[libbase.UnifiedAtom], bool]) → None¶: Selects atoms based on a predicate function.

deselect_atoms_if(pred: Callable[[libbase.UnifiedAtom], bool]) → None¶: Deselects atoms based on a predicate function.

select_connected_if(pred: Callable[[libbase.UnifiedAtom, libbase.UnifiedAtom, libbase.UnifiedBond], bool]) → None¶

Selects atoms bonded to the currently selected atoms based on a predicate function.

The predicate function is called on a pair of atoms and their connecting bond. The unselected atom attached to the bond only becomes selected of the predicate function returns True.

select_molecule_if(pred: Callable[[libbase.UnifiedAtom, libbase.UnifiedAtom, libbase.UnifiedBond], bool]) → None¶

Using the bonds and a predicate function, selects all molecules that include a currently selected atom.

This is equivalent to repeatedly calling select_connected_if with the same predicate until the selection stops growing.

Atomic Properties¶

Atomic properties provide additional, customizable information for each atom. They are organized into categories, each relevant to different computation engines or modules.

Generic Properties: These are common to all atoms, like mass.
Engine/Module-Specific Properties: Specific to certain computation engines or modules. For example, properties relevant only to the ADF engine would be in the adfprops group and prefixed with ‘adf’ in the System Block.

You can toggle these property groups on or off. Here’s an example, where we initialize some atomic properties in the System Block:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem

mol = ChemicalSystem(
"""
System
     Atoms
          O 0.0 0.0 0.0 forcefield.type=O_water
          H 1.0 0.0 0.0 forcefield.type=H_water mass=2.014
          H 0.0 1.0 0.0 forcefield.type=H_water
     End
End
"""
)

# The Generic Properties are always available (and have default values if not specified)
print(mol.atoms[1].mass) # Outputs 2.014
print(mol.atoms[2].mass) # Outputs 1.00798

# The forcefield properties are enabled because there were defined in the system block
print(mol.atomic_properties_enabled('forcefield')) # Outputs True

# Accessing the forcefield atomic properties
print(mol.atoms[0].forcefield.type) # Outputs O_water

# The adf properties are not enabled
print(mol.atomic_properties_enabled('adf')) # Outputs False

# Enable the 'adf' property group before using it:
mol.enable_atomic_properties("adf")
mol.atoms[0].adf.f = "fragment_1"

List of properties and utility methods in the ChemicalSystem:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

A class representing a chemical system in the Amsterdam Modeling Suite.

property adfprops¶: An optional list of all ADFProperties for all atoms.

property forcefieldprops¶: An optional list of all ForcefieldProperties for all atoms.

property guiprops¶: An optional list of all GUIProperties for all atoms.

enable_atomic_properties(group_prefix: Union[Literal[gui], Literal[adf], Literal[forcefield]]) → None¶: Enables the use of a group of atomic properties.

disable_atomic_properties(group_prefix: Union[Literal[gui], Literal[adf], Literal[forcefield]]) → None¶: Disables the use of a group of atomic properties. Any set properties within the group will be discarded.

atomic_properties_enabled(group_prefix: Union[Literal[gui], Literal[adf], Literal[forcefield]]) → bool¶: Checks if a group of atomic properties is enabled or not.

has_ghost_atoms() → bool¶: Checks whether the ChemicalSystem contains any Ghost atoms.

Properties with in the atom class:

class UnifiedAtom

class UnifiedAtom(Z: int)

class UnifiedAtom(element: UnifiedElement)

class UnifiedAtom(symbol: T_AtomSymbol)

Class representing a single atom and its properties.

property mass¶: The mass of the atom in dalton.

property is_ghost¶: Whether the atom is a Ghost atom.

property adf¶: Atomic properties used by the ADF engine.

property forcefield¶: Atomic properties used by the AMS Forcefield engine.

property gui¶: Atomic properties used for visualization in the GUI.

Classes for the Engine/Module-Specific Properties:

class UnifiedGUIProperties¶

Atomic properties used for visualization in the GUI.

clear() → None¶: Unsets any GUI atomic properties.

property color¶: The color of an atom.

empty() → bool¶: Returns whether any of the GUI atomic properties is set.

property radius¶: The radius of an atom for visualization purposes.

class UnifiedADFProperties¶

Atomic properties used by the ADF engine.

property ChgU¶: Optional atom centered point charge for ADF’s 3D-RISM method.

property EpsU¶: Lennard-Jones parameter for ADF’s 3D-RISM method.

property SigU¶: Lennard-Jones parameter for ADF’s 3D-RISM method.

clear() → None¶: Unsets any ADF atomic properties.

empty() → bool¶: Returns whether any of the ADF atomic properties is set.

property f¶: Specifies which fragment an atom belongs to; see the Fragments block in the ADF engine input.

class UnifiedForcefieldProperties¶

Atomic properties used by the AMS Forcefield engine.

property charge¶: The charge of an atom used for electrostatics by the Forcefield engine.

clear() → None¶: Unsets any Forcefield atomic properties.

empty() → bool¶: Returns whether any of the Forcefield atomic properties is set.

property type¶: The atom type to be used by the Forcefield engine.

Converting to and from PLAMS Molecules¶

ChemicalSystems can PLAMS’ molecules can be converted into each other using the following conversion functions: scm.utils.conversions.plams_molecule_to_chemsys and scm.utils.conversions.chemsys_to_plams_molecule.

Example:

from scm.libbase import UnifiedChemicalSystem as ChemicalSystem
from scm.plams import *
from scm.utils.conversions import plams_molecule_to_chemsys, chemsys_to_plams_molecule

# Create a PLAMS molecule from a SMILES string:
plams_molecule = from_smiles('C-C')

# Convert from a PLAMS molecule to a ChemicalSystem:
chemical_system = plams_molecule_to_chemsys(plams_molecule)

# Convert back from ChemicalSystem to a PLAMS molecule:
plams_molecule = chemsys_to_plams_molecule(chemical_system)