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.

See the UnifiedChemicalSystem API docs for a complete overview.

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 coincides with its center of mass
mol.translate(-mol.center_of_mass())

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

Reading and writing, File formats¶

File formats overview:

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.

All text based formats can also be written using Python format strings:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

__format__(format_spec: str) → str¶

Formats a ChemicalSystem into a string representation.

The format_spec string starts with and identifier for the format to write, followed and optional : and a list of space separated key=value pairs configuring options of the format. E.g. the following would produce the AMS System block format in internal units (bohr) with the string “H2O” as a System name in the block header:

cs = ChemicalSystem(...)
print(f"{cs:in:units=internal name=H2O}")

Would produce the output:

System H2O
    Atoms [bohr]
        O 0 0 -0.7262847342654172
        H 0 1.4669943905469853 0.3631423765813393
        H 0 -1.4669943905469853 0.3631423765813393
    End
End

The following formats and options are supported at the moment:

• in for writing the AMS System block. This is the default format, if none is specified. The following options are supported:

  • name=... to put an arbitrary string as the system’s name into the block header.

  • units=[default|internal] to switch between default units and the units used by the ChemicalSystem internally. Internally the ChemicalSystem uses atomic units (e.g. bohr for lengths). Printing the System block in internal units avoids a possible loss in precision in the unit conversion and guarantees an exact ChemicalSystem -> str -> ChemicalSystem round-trip, i.e.:

    cs = ChemicalSystem(...)
    assert ChemicalSystem(f"{cs:in:units=internal}") == cs
    

  • skip=[gui|adf|band|forcefield|dftb|reaxff|qe|...] to avoid printing of a particular property group in the end-of-line string of an atom in the System%Atoms subblock. Multiple groups may be specified as a comma separated list.

  • unused_atomic_properties=[drop|include] to determine whether unused atomic property groups should be written out as Modify%EnableAtomicProperties entries. The default is not to do this, aka drop. This means that unused atomic property groups get lost in the ChemicalSystem -> str -> ChemicalSystem round-trip. If preserving them is imporant, this can be achieved by setting this option to include.

• xyz for writing a plain XYZ file without lattice or atomic properties. This format has no options.

• extended_xyz for writing the AMS extended XYZ format. This format has no options.

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

Elements¶

There are two classes for handling the elements of atoms in ChemicalSystems: UnifiedElements and UnifiedElement.

• UnifiedElements contains read-only static attributes for each element supported in AMS, accessible by their capitalized names. For example:

• UnifiedElements.Hydrogen: Represents the element Hydrogen.

• UnifiedElements.Carbon: Represents the element Carbon.

Each attribute is an instance of UnifiedElement corresponding to that element. Additionally, the UnifiedElements class also provides class methods to retrieve the elements by symbol or atomic number.

Atoms¶

The ChemicalSystem class maintains an array of instances of the UnifiedAtom class 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: UnifiedAtomList¶

A list of all Atom instances that are part of this ChemicalSystem.

property coords: AngstromCoordsArray¶

The coordinates of the atoms in angstrom.

contains_atom(atom: UnifiedAtom) → bool¶

Checks if an Atom instance is part of a ChemicalSystem.

atom_index(atom: UnifiedAtom) → int¶

Given an Atom instance, returns its index in ChemicalSystem.atoms.

add_atom(atom: UnifiedAtom) → None¶

add_atom(atom: UnifiedAtom, coords: ArrayLike, unit: str = 'angstrom') → None

add_atom(Z: int, coords: ArrayLike, unit: str = 'angstrom') → None

add_atom(element: UnifiedElement, coords: ArrayLike, unit: str = 'angstrom') → None

add_atom(symbol: UnifiedAtom.T_AtomSymbol, coords: ArrayLike, unit: str = 'angstrom') → 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: 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 changing the order of atoms within a ChemicalSystem:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

reorder_atoms(atom_order: List[int]) → None¶

Reorders the system’s atoms given a list of atom indices in the desired order.

The atom_order shall be a list of integers. The list shall be of length num_atoms and its values a permutation of range(0,num_atoms). It shall contain the indices of the atoms in the original system in the order in which the atoms should occur after reordering.

Example: for a 6 atom system the atom_order = [4, 5, 0, 1, 2, 3] would have the effect of moving the last two atoms to the front of the system. (One might also look at atom_order as a mapping from the atom index in the reordered system to the atom index in the original system. Using the above example atom_order[1] == 5, which is the index that the new atom at index 1 had in the original system.)

sorted_atom_order(comp: atom_comparator_func_t) → List[int]¶

Returns the atom order based on a user defined comparator function as a list of atom indices.

The returned list can be passed into the reorder_atoms method to actually sort the atoms according to comp. See also documentation of the reorder_atoms and sort_atoms method.

sort_atoms(comp: atom_comparator_func_t) → None¶

Sorts the system’s atoms according to a user defined comparator function.

Atoms are sorted in ascending order where the user defined comparator function implements the < operator.

Example for sorting atoms descending in their atomic number, i.e. heavy atoms first:

chemsys.sort_atoms(lambda A, B: A.Z > B.Z)

Note that the atom sorting is stable: two atoms for which comp(A,B) returns False will keep their relative order in the sorted system.

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 a compatible lattice or no lattice at all. Lattices compatibility is checked with the lattice.is_close() method. If both systems have a lattice, and the lattice is compatible but not exactly the same, the lattice of the original system (self) is kept. 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: ArrayLike) → 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.

copy() → UnifiedAtom

Creates a deep copy of the Atom.

has_identical_properties(other: UnifiedAtom) → bool

has_identical_properties(other: UnifiedAtom, properties: List[T_UnifiedPropertyGroup]) → bool

Helper for @overload to raise when called.

in_chemicalsystem() → bool

Check if this Atom instance is part of a ChemicalSystem.

transmute(Z: int)

transmute(element: UnifiedElement)

transmute(symbol: T_AtomSymbol)

Helper for @overload to raise when called.

property Z: int

The atomic number of the atom.

property adf: UnifiedADFProperties | None

Atomic properties used by the ADF engine.

property band: UnifiedBANDProperties | None

Atomic properties used by the AMS BAND engine.

property coords: AngstromCoordsArray

The coordinates of the atom in angstrom.

property dftb: UnifiedDFTBProperties | None

Atomic properties used by the AMS DFTB engine.

property element: UnifiedElement

Returns the atom’s element.

property forcefield: UnifiedForcefieldProperties | None

Atomic properties used by the AMS Forcefield engine.

property gui: UnifiedGUIProperties | None

Atomic properties used for visualization in the GUI.

property is_ghost: bool

Whether the atom is a Ghost atom.

Ghost atoms are a special construct in AMS for calculating basis set superposition errors. Please refer to the ADF and BAND manuals for details.

property mass: float

The mass of the atom in dalton.

property qe: UnifiedQEProperties | None

Atomic properties used by the AMS QuantumESPRESSO engine.

property reaxff: UnifiedReaxFFProperties | None

Atomic properties used by the AMS ReaxFF engine.

property symbol: T_AtomSymbol

The atom’s symbol (e.g. ‘C’ or ‘Au’). May be ‘Gh.C’ for ghost atoms.

Atomic species

One often needs to group a system’s atoms into sets of indistinguishable atoms, essentially answering the question: which atoms are considered the same?

In AMS these groups of equivalent atoms are often called “atomic species”. Note that the assignment of atoms into species is application dependent: an electronic structure code may only care about the elements of the atoms, while a molecular dynamics driver will also need to distinguish the different isotopes of the same element. Engines implementing classical forcefields will even need to distinguish atoms based on their forcefield atom types and formal charges.

The ChemicalSystem class allows users to provide a binary predicate function that implements the comparison of two atoms, to allow for any definition of atomic species:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

determine_species(comp: atom_comparator_func_t) → Tuple[int, List[int], List[int]]¶

Determines the present atomic species, based on a user defined comparator function.

The comparator shall return True for two atoms if they are to be considered of the same species.

Returns a 3-tuple consisting of:

1. num_species is the total number of detected distinct species.

2. species_prototypes a list of size num_species with the index of the first atom in the system that was recognized as being of a distinct species. This atom is considered the species’ prototype atom and all atoms of the same species compare equal to it by the user defined comparator function comp.

3. atom_to_species a list of size num_atoms mapping each atom to its species.

See also the related get_species function, which returns the same data as a list of atom indices for each species.

In practice is usually easiest to implement the comparison as a lambda function right at the call site:

# Distinguish isotopes as separate species:
... = chemsys.determine_species(lambda A, B: A.Z == B.Z and A.mass == B.mass)

Note that you can access the ChemicalSystem instance from within the lambda, allowing access to things like regions or bonds, that are not part of the Atom class itself.

# Elements as species, but separately for region "substrate" and the rest:
... = chemsys.determine_species(
     lambda A, B: A.Z == B.Z and
     chemsys.is_atom_in_region(A, "substrate") == chemsys.is_atom_in_region(B, "substrate")
)

# Distinguish species based on the number of bonds
# e.g. all 4-bonded atoms are equivalent:
... = chemsys.determine_species(
   lambda A, B: chemsys.bonds.num_bonds(chemsys.atom_index(A)) == chemsys.bonds.num_bonds(chemsys.atom_index(B))
)

If you primarily want to iterate over all atoms within a species, the get_species function may return the same information in a more convenient format:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

get_species(comp: atom_comparator_func_t) → List[List[int]]¶

Determines the present atomic species, based on a user defined comparator function.

The comparator shall return True for two atoms if they are to be considered of the same species.

Returns a list of atom indices for each species, allowing for easy iteration over all atoms belonging to a species. Atom indices within each species are in ascending order.

for species in chemsys.get_species():
    for iat in species:
        # do something with chemsys.atoms[iat]

See also the related determine_species function, which returns the same information in a different format.

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: UnifiedADFPropertiesList | None¶

An optional list of all ADFProperties for all atoms.

property bandprops: UnifiedBANDPropertiesList | None¶

An optional list of all BANDProperties for all atoms.

property qeprops: UnifiedQEPropertiesList | None¶

An optional list of all QEProperties for all atoms.

property dftbprops: UnifiedDFTBPropertiesList | None¶

An optional list of all DFTBProperties for all atoms.

property reaxffprops: UnifiedReaxFFPropertiesList | None¶

An optional list of all ReaxFFProperties for all atoms.

property forcefieldprops: UnifiedForcefieldPropertiesList | None¶

An optional list of all ForcefieldProperties for all atoms.

property guiprops: UnifiedGUIPropertiesList | None¶

An optional list of all GUIProperties for all atoms.

enable_atomic_properties(group_prefix: Literal['gui', 'adf', 'band', 'forcefield', 'dftb', 'reaxff', 'qe']) → None¶

Enables the use of a group of atomic properties.

disable_atomic_properties(group_prefix: Literal['gui', 'adf', 'band', 'forcefield', 'dftb', 'reaxff', 'qe']) → None¶

Disables the use of a group of atomic properties. Any set properties within the group will be discarded.

atomic_properties_enabled(group_prefix: Literal['gui', 'adf', 'band', 'forcefield', 'dftb', 'reaxff', 'qe']) → bool¶

Checks if a group of atomic properties is enabled or not.

has_ghost_atoms() → bool¶

Checks whether the ChemicalSystem contains any Ghost atoms.

Ghost atoms are a special construct in AMS for calculating basis set superposition errors. Please refer to the ADF and BAND manuals for details.

Properties within 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: float¶

The mass of the atom in dalton.

property is_ghost: bool¶

Whether the atom is a Ghost atom.

Ghost atoms are a special construct in AMS for calculating basis set superposition errors. Please refer to the ADF and BAND manuals for details.

property adf: UnifiedADFProperties | None¶

Atomic properties used by the ADF engine.

property band: UnifiedBANDProperties | None¶

Atomic properties used by the AMS BAND engine.

property qe: UnifiedQEProperties | None¶

Atomic properties used by the AMS QuantumESPRESSO engine.

property dftb: UnifiedDFTBProperties | None¶

Atomic properties used by the AMS DFTB engine.

property reaxff: UnifiedReaxFFProperties | None¶

Atomic properties used by the AMS ReaxFF engine.

property forcefield: UnifiedForcefieldProperties | None¶

Atomic properties used by the AMS Forcefield engine.

property gui: UnifiedGUIProperties | None¶

Atomic properties used for visualization in the GUI.

Classes for the Engine/Module-Specific Properties:

class UnifiedADFProperties¶

Atomic properties used by the ADF engine.

clear() → None¶

Unsets any ADF atomic properties.

copy() → UnifiedADFProperties¶

Creates a deep copy of the ADFProperties.

empty() → bool¶

Returns whether any of the ADF atomic properties is set.

property ChgU: float | None¶

Optional atom centered point charge for ADF’s 3D-RISM method.

property EpsU: float | None¶

Lennard-Jones parameter for ADF’s 3D-RISM method.

property R: str | None¶

The radius of the atom for constructing the COSMO surface.

If the string represents a real number that is used as the radius, otherwise the string is used for looking up the radius in the Solvation%Radii block.

property SigU: float | None¶

Lennard-Jones parameter for ADF’s 3D-RISM method.

property f: str | None¶

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

property fg: str | None¶

TODO

property fs: str | None¶

TODO

property imol: str | None¶

TODO

property nuclear_charge: float | None¶

Alternative nuclear charge for an atom.

property type: str | None¶

The atom sub-type to be used by the ADF engine.

This can be used to make the ADF engine treat two atoms of the same elements as separate atom types:

System
    Atoms
        C  x y z   adf.type=1
        C  x y z   adf.type=2
    End
End

This is used as a replacements for the “dot-notation” from ADF<2025, which is no longer allowed in ADF>=2025:

System
    Atoms
        C.1  x y z
        C.2  x y z
    End
End

property x: ndarray[Any, dtype[float64]] | None¶

Custom atomic x-axis to be used by the ADF engine.

property z: ndarray[Any, dtype[float64]] | None¶

Custom atomic z-axis to be used by the ADF engine.

class UnifiedBANDProperties¶

Atomic properties used by the AMS BAND engine.

clear() → None¶

Unsets any BAND atomic properties.

copy() → UnifiedBANDProperties¶

Creates a deep copy of the BANDProperties.

empty() → bool¶

Returns whether any of the BAND atomic properties is set.

property nuclear_charge: float | None¶

Alternative nuclear charge for an atom.

class UnifiedDFTBProperties¶

Atomic properties used by the AMS DFTB engine.

clear() → None¶

Unsets any DFTB atomic properties.

copy() → UnifiedDFTBProperties¶

Creates a deep copy of the DFTBProperties.

empty() → bool¶

Returns whether any of the DFTB atomic properties is set.

property Vext: float | None¶

The external potential at the atomic nucleus used by the DFTB engine.

class UnifiedReaxFFProperties¶

Atomic properties used by the AMS ReaxFF engine.

clear() → None¶

Unsets any ReaxFF atomic properties.

copy() → UnifiedReaxFFProperties¶

Creates a deep copy of the ReaxFFProperties.

empty() → bool¶

Returns whether any of the ReaxFF atomic properties is set.

property charge: float | None¶

The charge of an atom used for electrostatics by the ReaxFF engine.

class UnifiedForcefieldProperties¶

Atomic properties used by the AMS Forcefield engine.

clear() → None¶

Unsets any Forcefield atomic properties.

copy() → UnifiedForcefieldProperties¶

Creates a deep copy of the ForcefieldProperties.

empty() → bool¶

Returns whether any of the Forcefield atomic properties is set.

property charge: float | None¶

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

property type: str | None¶

The atom type to be used by the Forcefield engine.

class UnifiedGUIProperties¶

Atomic properties used for visualization in the GUI.

clear() → None¶

Unsets any GUI atomic properties.

copy() → UnifiedGUIProperties¶

Creates a deep copy of the GUIProperties.

empty() → bool¶

Returns whether any of the GUI atomic properties is set.

property color: Tuple[int, int, int] | None¶

The color of an atom.

property radius: float | None¶

The radius of an atom for visualization purposes.

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: float¶

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

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:

Methods for changing distances between atoms:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

get_distance(a: int, b: int, unit: str = 'angstrom') → float¶

get_distance(a: UnifiedAtom, b: UnifiedAtom, unit: str = 'angstrom') → float

Measures the distance between two atoms.

set_distance(a: int, b: int, dist: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.LighterPartMoves) → None¶

set_distance(a: int, b: int, dist: float, unit: str, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.LighterPartMoves) → None

set_distance(a: int, b: int, dist: float, moving_atoms: ArrayLike) → None

set_distance(a: int, b: int, dist: float, unit: str, moving_atoms: ArrayLike) → None

Sets the distance between two atoms.

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

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

Methods for changing angles:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

set_angle(a: int, b: int, c: int, angle: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None¶

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

set_angle(a: int, b: int, c: int, angle: float, rotation_axis: ArrayLike, moving_atoms: ArrayLike) → None

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.

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.

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.

Methods for changing dihedral angles:

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

set_dihedral(a: int, b: int, c: int, d: int, angle: float, policy: InternalCoordinateManipulationPolicy = InternalCoordinateManipulationPolicy.SecondPartMoves) → None¶

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

set_dihedral(a: int, b: int, c: int, d: int, angle: float, moving_atoms: ArrayLike) → None

set_dihedral(a: int, b: int, c: int, d: int, angle: float, unit: str, moving_atoms: ArrayLike) → None

Sets the dihedral between four atoms.

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

class UnifiedChemicalSystem

class UnifiedChemicalSystem(system_block: str)

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

translate(shift: ArrayLike, unit: str = 'angstrom') → 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) → ndarray[Any, dtype[float64]]¶

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

classmethod rmsd(a: UnifiedChemicalSystem, b: UnifiedChemicalSystem, 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 UnifiedChemicalSystems must have the same number of atoms.

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

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

Position of the geometric center.

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

Position of the center of mass.

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.

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.

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.

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 = '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).

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. See the equivalent symmetrize_cell method for working with periodic structures.

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.

check_molecule_symmetry(label: T_MoleculeSchoenfliesSymbols | str, tolerance: float = 1e-07) → 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.

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.

get_inverted() → UnifiedBond¶

Returns a copy of the bond with inverted lattice displacements.

property lattice_displacements: ndarray[Any, dtype[int64]]¶

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: float¶

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

The Bonds Class

Methods to inspect the bonds:

class UnifiedBonds(num_atoms: int, num_lattice_vectors: int)

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

num_bonds() → int¶

num_bonds(atom: int) → int

num_bonds(from_atom: int, to_atom: int) → int

If called without an atom index, returns the total number of bonds between all atoms. If called with an atom index, returns the number of bonds of that atom. if called with two atom indices, returns the number of bonds between the two atoms.

atoms_are_bonded(from_atom: int, to_atom: int) → bool¶

Checks whether two atoms are bonded or not.

any_bond_between(from_atoms: List[int], to_atoms: List[int]) → bool¶

Checks if there is any direct bond between two groups of atoms.

get_bonded_atoms(atom: int) → List[int]¶

Returns a list of unique atom indices that are bonded to the given atom index.

property num_atoms: int¶

The total number of atoms in the system.

property num_lattice_vectors: int¶

The number of lattice vectors of the system.

Methods to iterate over bonds:

class UnifiedBonds(num_atoms: int, num_lattice_vectors: int)

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

__len__() → int¶

Returns the total number of bonds between all atoms.

__iter__() → Iterator[Tuple[int, int, UnifiedBond]]¶

Returns an iterator over all bonds.

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

get_bonds_for_atom(from_atom: int) → Iterator[Tuple[int, int, 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.

get_bonds_between_atoms(from_atom: int, to_atom: int) → Iterator[Tuple[int, int, UnifiedBond]]¶

Returns an iterator over all bonds from between two atoms given by their indices.

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

Methods to modify the bonds:

class UnifiedBonds(num_atoms: int, num_lattice_vectors: int)

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

add_bond(from_atom: int, to_atom: int, bond: UnifiedBond) → None¶

Adds a bond between two atoms, given their indices.

add_bonds(from_atoms: ArrayLike, to_atoms: ArrayLike, bonds: List[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.

remove_bond(bidx: int | UnifiedBond) → None¶

Removes a bond, given either its index in the .bonds attribute or an instance of a Bond.

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.

remove_bonds_between_atoms(from_atom: int, to_atom: int) → None¶

Removes the bonds between two atoms, gived their indices.

clear_bonds() → None¶

Removes all bonds.

Direct read-only access to the modified Compressed Sparse Row format:

class UnifiedBonds(num_atoms: int, num_lattice_vectors: int)

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

property row_offset: ndarray[Any, dtype[int64]]¶

Bonds are stored as sparse matrix in a modified Compressed Spare Row format: The index in row_offset represents from_atom, and the value is an offset in column_index and bonds. The values in column_index represent to_atom. Because there can be duplicate indices (periodic bonds), it is not strictly a CSR matrix.

property column_index: ndarray[Any, dtype[int64]]¶

Bonds are stored as sparse matrix in a modified Compressed Spare Row format: The index in row_offset represents from_atom, and the value is an offset in column_index and bonds. The values in column_index represent to_atom. Because there can be duplicate indices (periodic bonds), it is not strictly a CSR matrix.

property bond_index: ndarray[Any, dtype[int64]]¶

Bonds are stored as sparse matrix in a modified Compressed Spare Row format: The index in row_offset represents from_atom, and the value is an offset in column_index and bonds. The values in column_index represent to_atom. Because there can be duplicate indices (periodic bonds), it is not strictly a CSR matrix.

property bonds: UnifiedBondList¶

Bonds are stored as sparse matrix in a modified Compressed Spare Row format: The index in row_offset represents from_atom, and the value is an offset in column_index and bonds. The values in column_index represent to_atom. Because there can be duplicate indices (periodic bonds), it is not strictly a CSR matrix.

Supported file formats and conversions to other formats:

class UnifiedBonds(num_atoms: int, num_lattice_vectors: int)

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

classmethod from_kf(kf: KFFile, section: str) → UnifiedBonds¶

Constructs and returns new Bonds from a section on a KF file.

write_kf(kf: KFFile, section: str, write_list_format: bool = True) → None¶

Writes Bonds to a section on a KF file.

The write_list_format argument determines whether the retrocompatible list format (with variables fromAtoms, toAtoms, bondOrders) is written. Otherwise the Compressed Sparse Row format is written directly.

classmethod from_list(num_atoms: int, num_lattice_vectors: int, from_atoms: ArrayLike, to_atoms: ArrayLike, bond_orders: ArrayLike, lattice_displacements: ArrayLike = [[]], indexing: int = 0) → UnifiedBonds¶

Creates a set of bonds from lists of atom indices, bond orders and optional lattice displacements.

The indexing argument can be used to switch to one-based atom indices, as they are used in the AMS System block format.

classmethod from_sparse(num_lattice_vectors: int, row_offset: ArrayLike, column_index: ArrayLike, bond_orders: ArrayLike, lattice_displacements: ArrayLike = [[]], indexing: int = 0) → UnifiedBonds¶

Low level construction of a set of bonds directly from the (modified) CRS storage.

The indexing argument can be used to switch to one-based atom indices, as they are used in the AMS System block format.

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: UnifiedBonds¶

The bonds of the system.

has_bonds() → bool¶

Checks if the system has any bonding information.

guess_bonds() → None¶

Guesses bonds based on the atomic elements and the geometry. Keeps existing bonds.

set_lattice_displacements_from_minimum_image_convention() → None¶

Applies the minimum image convention which sets the lattice_displacements for all bonds to be the shortest possible.

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() → ndarray[Any, dtype[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: UnifiedAtom, to_atom: UnifiedAtom) → bool

Checks if removing the bonds between two atoms would cut the graph into two disjoint subgraphs.

atom_is_in_ring(atom: int | 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.

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:

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

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

3. 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: int | UnifiedAtom) → List[float]¶

shortest_path_lengths_from(from_atom: int | 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) → List[int] | None¶

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

shortest_path_between(from_atom: UnifiedAtom, to_atom: UnifiedAtom, dist_func: molgraph_edgeweight_func_t) → List[int] | None

shortest_path_between(from_atom: UnifiedAtom, to_atom: UnifiedAtom) → List[int] | None

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: UnifiedAtom, to_atom: UnifiedAtom, dist_func: molgraph_edgeweight_func_t) → float

shortest_path_length_between(from_atom: UnifiedAtom, to_atom: 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: UnifiedLattice¶

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, unit: str = 'angstrom')

A class representing the lattice of a ChemicalSystem.

cartesian_to_fractional(cartesian_coord: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Convert Cartesian coordinates to fractional coordinates for a single point. No conversion will take place in non-periodic directions.

cartesians_to_fractionals(cartesian_coords: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Convert Cartesian coordinates to fractional coordinates for a set of points. No conversion will take place in non-periodic directions.

copy() → UnifiedLattice¶

Creates a deep copy of the Lattice.

fractional_to_cartesian(fractional_coord: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Convert fractional coordinates to Cartesian coordinates for a single point. No conversion will take place in non-periodic directions.

fractionals_to_cartesians(fractional_coords: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Convert fractional coordinates to Cartesian coordinates for a set of points. No conversion will take place in non-periodic directions.

classmethod from_kf(kf: KFFile, section: str) → UnifiedLattice¶

Constructs and returns a new Lattice from a section on a KF file.

get_angles(unit: str = 'rad') → ndarray[Any, dtype[float64]]¶

Get the angles between the lattice vectors (note: always returns 3 numbers).

get_cell_depth(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Return unit cell depth, or thickness, in each dimension.

The cell depth is the distance between the two opposite bounding planes of the parallelepiped spanned by the lattice vectors. For three lattice vectors v_i, v_j, and v_k, the i-th component of the cell depth is the height of lattice vector i over the plane spanned by v_j and v_k.

Always returns a 3-component vector. For < 3D periodic systems, the remaining components for the non-periodic directions at the end are set to 0.

get_lengths(unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Get the length of the lattice vectors (note: always returns 3 numbers).

get_reciprocal_lattice_vectors(unit: str = 'angstrom-1') → ndarray[Any, dtype[float64]]¶

Get the reciprocal lattice vectors as a numpy array.

get_volume(unit: str = 'angstrom') → float¶

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

is_close(other: UnifiedLattice, tol: float = 0.001, unit: str = 'angstrom') → bool¶

is_close(other: UnifiedLattice, length_tol: float, angle_tol: float, length_unit: str = 'angstrom', angle_unit: str = 'rad') → bool

Checks if two sets of lattice vectors are close to each other.

is_orthogonal() → bool¶

Is the lattice orthogonal?

map_vector_to_central_cell(vec: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Maps a vector into the central cell, such that its fractional coordinates are in range [-0.5,+0.5). No mapping will take place in non-periodic directions.

map_vectors_to_central_cell(vecs: ArrayLike, unit: str = 'angstrom') → ndarray[Any, dtype[float64]]¶

Maps a set of vectors into the central cell, such that their fractional coordinates are in range [-0.5,+0.5). No mapping will take place in non-periodic directions.

write_kf(kf: KFFile, section: str) → None¶

Writes a Lattice to a section on a KF file.

property num_vectors: int¶

Returns the number of lattice vectors (i.e. the number of PBC).

property vectors: AngstromCoordsArray¶

The lattice vectors in angstrom.