Geometry modification

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 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

Manipulation in internal coordinates

class InternalCoordinateManipulationPolicy(value)

Enum representing the options for how atoms move when changes are made in internal coordinates.

Values:

  • FirstPartMoves:

  • SecondPartMoves:

  • LighterPartMoves:

  • HeavierPartMoves:

Distances between atoms

ChemicalSystem.get_distance(a: int, b: int, unit: str = 'angstrom') → float
ChemicalSystem.get_distance(a: Atom, b: Atom, unit: str = 'angstrom') → float

Measures the distance between two atoms.

For periodic systems, the returned distance is the shortest distance under minimum image convention.

ChemicalSystem.set_distance(a: int, b: int, dist: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_distance(a: int, b: int, dist: float, unit: str, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_distance(a: int, b: int, dist: float, moving_atoms: ArrayLike) → None
ChemicalSystem.set_distance(a: int, b: int, dist: float, unit: str, moving_atoms: ArrayLike) → None

Sets the distance between two atoms.

ChemicalSystem.moving_atoms_for_distance_change(a: int, b: int, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → NDArray[np.int_]

Determines which atoms will move when changing the distance between two atoms.

Angles

ChemicalSystem.get_angle(a: int, b: int, c: int, unit: str = 'radians') → float
ChemicalSystem.get_angle(a: Atom, b: Atom, c: Atom, unit: str = 'radians') → float

Measures the angle A-B-C between three atoms.

For periodic systems, the vectors r_BA and r_BC used for the calculation of the angle are taken under minimum image convention.

The returned angle is in the range [0:pi] radians.

ChemicalSystem.set_angle(a: int, b: int, c: int, angle: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_angle(a: int, b: int, c: int, angle: float, unit: str, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_angle(a: int, b: int, c: int, angle: float, rotation_axis: ArrayLike, moving_atoms: ArrayLike) → None
ChemicalSystem.set_angle(a: int, b: int, c: int, angle: float, unit: str, rotation_axis: ArrayLike, moving_atoms: ArrayLike) → None

Sets the angle between three atoms.

ChemicalSystem.moving_atoms_for_angle_change(a: int, b: int, c: int, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → NDArray[np.int_]

Determines which atoms will move when changing the angle between three atoms.

ChemicalSystem.rotation_axis_for_angle_change(a: int, b: int, c: int) → ndarray[Any, dtype[float64]]

Determines the rotation axis for changing the angle between three atoms.

The returned vector is normalized and and oriented such that r_BA, r_BC and the rotation axis follow the right-hand rule.

Dihedral angles

ChemicalSystem.get_dihedral(a: int, b: int, c: int, d: int, unit: str = 'radians') → float
ChemicalSystem.get_dihedral(a: Atom, b: Atom, c: Atom, d: Atom, unit: str = 'radians') → float

Measures the dihedral angle A-B-C-D between four atoms.

For periodic systems, the vectors r_AB, r_BC and r_CD used for the calculation of the dihedral angle are taken under minimum image convention.

The returned angle is in the range ]-pi:pi] radians.

ChemicalSystem.set_dihedral(a: int, b: int, c: int, d: int, angle: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_dihedral(a: int, b: int, c: int, d: int, angle: float, unit: str, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None
ChemicalSystem.set_dihedral(a: int, b: int, c: int, d: int, angle: float, moving_atoms: ArrayLike) → None
ChemicalSystem.set_dihedral(a: int, b: int, c: int, d: int, angle: float, unit: str, moving_atoms: ArrayLike) → None

Sets the dihedral between four atoms.

ChemicalSystem.moving_atoms_for_dihedral_change(a: int, b: int, c: int, d: int, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → NDArray[np.int_]

Determines which atoms will move when changing the dihedral between four atoms.

Other methods

ChemicalSystem.translate(shift: ArrayLike, unit: str = 'angstrom') → None

Translates all atoms in the ChemicalSystem by a vector.

ChemicalSystem.rotate(rot_mat: ArrayLike) → None

Rotate the system according to the rotation matrix.

ChemicalSystem.align_to(other: ChemicalSystem) → 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 ChemicalSystems must have the same number of atoms.

  • The atoms in the two chemical system must be in the same order.

classmethod ChemicalSystem.rotation_matrix_minimizing_rmsd(a: ChemicalSystem, b: ChemicalSystem) → ndarray[Any, dtype[float64]]

Given two chemical systems, returns the rotation matrix that minimizes the RMSD.

classmethod ChemicalSystem.rmsd(a: ChemicalSystem, b: ChemicalSystem, align: bool, unit: str = 'angstrom') → 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 ChemicalSystems must have the same number of atoms.

  • The atoms in the two chemical system must be in the same order.

ChemicalSystem.geometric_center(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]

Position of the geometric center.

ChemicalSystem.center_of_mass(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]

Position of the center of mass.

ChemicalSystem.bounding_box_volume(unit: str = 'angstrom3') → float

Volume of the bounding box.

For periodic systems, the entire extend of the lattice is used in the periodic directions. Example: for a slab (in the xy-plane), the bounding box volume is the cell’s area times the extent of the atomic coordinates in the z-direction.

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

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

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

ChemicalSystem.moments_of_inertia() → Tuple[ndarray[Any, dtype[float64]], ndarray[Any, dtype[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 (units amu*angstrom^2) sorted by magnitude

  • a 3x3 np.array 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.

ChemicalSystem.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.

ChemicalSystem.perturb_coordinates(noise_level: float, unit: str = 'angstrom') → 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).

ChemicalSystem.distance_matrix(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]
ChemicalSystem.distance_matrix(atom_indices: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]
ChemicalSystem.distance_matrix(from_indices: ArrayLike, to_indices: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]

Returns the distance matrix between all, or a subset of atoms in the system.

ChemicalSystem.squared_distance_matrix(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]
ChemicalSystem.squared_distance_matrix(atom_indices: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]
ChemicalSystem.squared_distance_matrix(from_indices: ArrayLike, to_indices: ArrayLike, unit: str = 'angstrom2') → ndarray[Any, dtype[float64]]

Returns the matrix of squared distances between all, or a subset of atoms in the system.

ChemicalSystem.get_distances_to_point(point: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]

Returns the distance of each atom from a given point.

Note that the unit argument determines the unit in which distances are returned. The coordinates of the point are always expected to be in angstrom.

ChemicalSystem.symmetrize_molecule(tolerance: float = 0.05) → T_MoleculeSchoenfliesSymbols

Symmetrizes and reorients a molecule. Returns the Schoenflies symbol of the point group to which the system was symmetrized.

Symmetrizes the atomic coordinates to machine precision and rototranslates the molecule to the AMS standard orientation. This allows optimal use of the system’s symmetry for calculations with the AMS driver and its engines. The standard orientation has the following properties:

  • The origin is a fixed point of the symmetry group.

  • The z-axis is the main rotation axis.

  • The xy-plane is the sigma_h plane (axial groups, C(s)).

  • The x-axis is a C_2 axis (D symmetries).

  • The xz-plane is a sigma_v plane (C_nv symmetries).

  • In T_d and O_h the z-axis is a fourfold axis (S_4 and C_4, respectively) and the (111)-direction is a threefold axis.

See also the Symmetry section in the AMS driver manual’s appendix for more information about usage of molecular symmetry in AMS.

Calling this method may be useful if the system is almost symmetric or to rototranslate a symmetric molecule into the AMS standard orientation.

Note that at the moment it is not possible to symmetrize systems with more than 2000 atoms. Calling this method with a periodic system throws an exception. Calling this method in a molecule with an electrostatic embedding is not supported yet and throws an exception. See the equivalent symmetrize_cell method for working with periodic structures.

ChemicalSystem.symmetrize_molecule_to(label: T_MoleculeSchoenfliesSymbols | str, tolerance: float = 1e-06) → None

Symmetrizes a molecule to a particular symmetry given by a Schoenflies symbol.

Symmetrizes the atomic coordinates to machine precision, but does NOT rototranslate the molecule in any way.

If the symmetrization fails, this routine throws a ChemicalSystemError. Since this method does not rototranslate the molecule, this may happen if the input structure is not in the AMS standard orientation. See the description of the symmetrize_molecule method for information about the AMS standard orientation.

Calling this method with a periodic system throws an exception. Calling this method in a molecule with an electrostatic embedding is not supported yet and throws an exception.

ChemicalSystem.symmetrize() → str

Symmetrizes a system, using either symmetrize_molecule or symmetrize_cell, depending on periodicity. Returns the Schoenflies symbol of the group to which the system was symmetrized.

ChemicalSystem.check_molecule_symmetry(label: T_MoleculeSchoenfliesSymbols | str, tolerance: float = 1e-06) → bool

Checks if a molecule has a particular symmetry given by a Schoenflies symbol.

Also checks that the molecule is in the AMS standard orientation for the specified point group. See the description of the symmetrize_molecule method for information about the AMS standard orientation.

The main use-case for this method is to check if a molecule is symmetric and oriented such that AMS can make full use of its symmetry.

Calling this method with a periodic system throws an exception. Calling this method in a molecule with an electrostatic embedding is not supported yet and throws an exception.