User Guide#
MolSystem is a Python package for storing molecular and periodic (crystalline) systems. It provides and object oriented interface over an underly SQL database that handles systems with multiple configurations, plus templates and subsets. You can loosely think of a system as being an abstract molecular or periodic system and the configurations as the physical instances with the coordinates, bonds, etc. that describe different structures, conformers, or frames in a trajectory.
Perceiving bonds#
Structures read from formats without connectivity – an extended XYZ trajectory frame,
for example – have atoms but no bonds, so anything that works with molecules (finding
molecules, keeping them whole across a periodic boundary, extracting clusters) has
nothing to go on. perceive_bonds() finds the bonds from the geometry and adds them to
the configuration:
n_bonds = configuration.perceive_bonds()
Two atoms are bonded when their distance is less than tolerance (default 1.2) times
the sum of their covalent radii (Pyykkö radii from mendeleev). Hydrogen is limited
to one bond, the shortest, so a hydrogen bond is never mistaken for a covalent one, and
the alkali and alkaline-earth metals are treated as ions and left unbonded. In a
periodic configuration every periodic image within reach is considered and each bond
records the cell offset of its partner, so molecules straddling the cell boundary are
bonded correctly and covalent crystals get all their bonds – primitive diamond, for
example, has four bonds between its two atoms through different images. All bonds are
single bonds: anything that derives a SMILES or Lewis structure from them would read
ethylene as [CH2+][CH2-]. For a molecule, assign_bond_orders() then gives the
bonds their orders and the atoms their formal charges:
n_assigned, n_failed = configuration.assign_bond_orders(charge=0)
Each molecule (connected fragment) gets the bond orders of RDKit’s
DetermineBondOrders, trying the charges 0, -1, +1, -2 and +2 in turn and keeping the
first that gives a valid structure; single alkali, alkaline-earth and halogen atoms get
their ionic charges. A molecule RDKit cannot handle, such as hypervalent PF6-, keeps its
single bonds, which are already its Lewis structure. Given the net charge, the molecules’
charges are checked against it and a mismatch is logged. Open-shell (a radical is
indistinguishable from an ion by geometry alone, so set the multiplicity first) and
periodic configurations are left alone. Read Structure does this for extended XYZ files
when it perceives the bonds.
The defaults can be adjusted per call:
configuration.perceive_bonds(
replace=True, # discard existing bonds first
tolerance=1.25, # looser criterion, e.g. hot MD frames
radii={"H": 0.35}, # override a covalent radius (Å)
max_bonds={"H": 1, "O": 2}, # per-element valence limits
exclude=(), # bond the metals too, e.g. an ionic crystal
)
Only the covalent radii method is available at present; the method argument is
there so that others, such as a Voronoi tessellation, can be added later.
Charges and spin in structure files#
A configuration carries its charge and spin multiplicity as properties of the structure as a whole:
configuration.charge = -1
configuration.spin_multiplicity = 1
The atoms may also carry formal charges of their own, in the optional
formal_charge attribute – the charge on the carboxylate oxygen of an acetate ion,
for instance, rather than on the ion as a whole:
configuration.atoms.get_column_data("formal_charge") # e.g. [0, 0, 0, -1]
The two are distinct and are written to a structure file as such: the molecular charge
and multiplicity as properties of the structure, and the formal charges on the
individual atoms that carry them – the M CHG line of an SDF file, say. Both come
back unchanged on reading the file, so the charges stay where they belong and a neutral,
closed-shell structure is not mistaken for a radical.
Radicals and ions in SMILES and RDKit#
to_smiles() and to_RDKMol() carry the charges and unpaired electrons into RDKit,
so that radicals and ions give the right SMILES:
configuration.from_smiles("[CH2]C") # charge 0, multiplicity 2
configuration.to_smiles(canonical=True) # '[CH2]C', the ethyl radical
from_smiles() sets the charge and spin multiplicity from the SMILES, and the formal
charges of the atoms. Since the hydrogens are all explicit, an atom short of its valence
is a radical, an anion or a cation: the formal charges are used when they add up to the
configuration’s charge, and otherwise, as for a structure read from a file without them,
charges and radical electrons are chosen that match the charge and multiplicity, with
the fewest charges, negative charges on electronegative atoms, positive ones on metals
and electropositive atoms, and charge-separated groups such as the nitro group
C[N+](=O)[O-] kept neutral. Hypervalent atoms need no charge of their own: sulfuric acid
is O=S(=O)(O)O. The same atoms and bonds are the radical [CH2]C, the cation [CH2+]C or the
anion [CH2-]C, depending on the configuration’s charge and multiplicity.
Properties and their units#
Every property has units, held as a string that Pint can interpret, and a property that is genuinely dimensionless – a statistical inefficiency, a count, a ratio – has an empty unit string rather than no units at all:
configuration.properties.units("temperature, inefficiency#LAMMPS#oplsaa+") # ""
The distinction matters because the units are what a value is converted to when a
property is stored or read, so a property whose units are undefined cannot be converted.
Anything that creates properties – including reading them back from a structure file –
records dimensionless units as "".
User tables#
Besides the structures, the database holds the tables a SEAMM flowchart builds – with
the Table step, or by storing results into a table – so that a job’s tables are in its
seamm.db with everything else. They are reached by name through
SystemDB.user_tables:
tables = db.user_tables
table = tables.create(
"energies",
columns=[("SMILES", "string", None), ("E (kJ/mol)", "float", None)],
index_column="SMILES",
)
row = table.append_row(SMILES="CCO")
table.set_cell(row, "E (kJ/mol)", -235.1)
df = table.to_dataframe() # a pandas copy, typed by the declared types
Each column has a declared type – boolean, integer, float, string or
json – and a default, which fills the column in rows that do not give a value. The
values are stored exactly as written (the text "1.0960" stays text) and read back
according to the declared type; values written to text columns are stored as text. A
table may have an index column whose values identify the rows; rows also have an
internal id, which is not meant for users, and a current row that SEAMM’s steps write
to (None means the next write appends a row).
The tables are stored under internal names, table_1, table_2, …, with
internal column names, so any display name is allowed, including names that differ only
in case. A registry table, _tables, records each table’s name, columns, index column,
current row and metadata, and a journal, _table_changes, records every change
(creating and dropping tables, adding columns, appending rows and setting values). Row
ids are never reused. Nothing is committed by these calls; the caller commits.
A database opened read-only (file:seamm.db?mode=ro) can be read but not written:
creating or changing a table raises PermissionError.