We also included covalent radii of elements from a
Wikipedia page,
which has data from
Cordero et al (2008), Covalent radii revisited, Dalton Trans. 21, 2832.
>>> gemmi.Element('Zr').covalent_r1.75
van der Waals radii taken from Wikipedia and cctbx:
CIF files that describe small-molecule and inorganic structures
can be read into a SmallStructure object.
Unlike macromolecular Structure,
SmallStructure has no hierarchy.
It is a flat list of atomic sites (SmallStructure::Site)
together with the unit cell and symmetry.
>>> site=SiC.sites[0]>>> site.label'Si1'>>> site.type_symbol'Si4+'>>> site.fract<gemmi.Fractional(0, 0, 0)>>>> site.occ1.0>>> site.u_iso# not specified here0.0>>> site.element# obtained from type_symbol 'Si4+'gemmi.Element('Si')>>> site.charge# obtained from type_symbol 'Si4+'4
The occupancies in small molecules normally represent the actual chemical
occupancy.
This differs from macromolecular crystallography, where models normally store
“crystallographic” occupancy – atoms on special positions have their occupancy
divided by the number of symmetry images in the same place.
This reduction of occupancy simplifies the calculation of structure factors.
and the function determine_and_set_spacegroup("S.H2") is automatically
run to set spacegroup:
>>> st.spacegroup<gemmi.SpaceGroup("P -3 m 1")>
determine_and_set_spacegroup() takes one argument, a string in which characters
specify what to use, and in what order, for space group determination:
S = symmetry operations stored in symops,
H = Hall symbol from spacegroup_hall (we compare symmetry operations
encoded in the Hall symbol, not the strings),
1 = H-M symbol; for space groups such as “P n n n” that have two origin
choices listed in the International Tables, use Origin Choice 1,
2 = H-M symbol, with Origin Choice 2 where applicable,
N = the space group number,
. (after S or H) = if the symmetry operations pass sanity checks,
stop and use them regardless of whether they correspond to one of
the settings tabulated in Gemmi.
If a symbol or operations match one of the 560+ space group settings tabulated
in Gemmi, spacegroup is set to this setting. Otherwise, if . is encountered
and the previous character (S or H) was evaluated to a valid set of symops,
it is assumed that these operations were correct: spacegroup is left null
and cell.images are set from the list of operations.
About 350 (out of 500,000+) entries in the COD use such settings.
Most of them have an unconventional choice of the origin
(e.g. “P 1 21 1 (a,b,c-1/4)”).
To use a different order of items than “S.H2”,
call determine_and_set_spacegroup() again:
>>> st.determine_and_set_spacegroup('H.1')
Errors such as an incorrect format of the symop triplets or of the Hall
symbol are silently ignored, and the consistency between different items
is not checked. That’s because this function is run when reading a file;
throwing an exception at that stage would prevent reading a file.
We have a separate function to check for errors and inconsistencies.
It returns a string, one line – one error:
>>> st.check_spacegroup()''
If the spacegroup setting used in a file is not tabulated in Gemmi,
you can still create a GroupOps object with symmetry operations:
>>> gemmi.GroupOps([gemmi.Op(o)foroinst.symops])<gemmi.GroupOps object at 0x...>>>> # or>>> gemmi.symops_from_hall(st.spacegroup_hall)<gemmi.GroupOps object at 0x...>
In C++ it would be similar, except that the following function
would be used to make gemmi::GroupOps from symops:
If your structure is stored in a macromolecular format (PDB, mmCIF)
you can read it first as macromolecular Structure
and convert it to SmallStructure:
Residues (monomers) and small molecule components of macromolecular models
are called chemical components.
Gemmi can use three sources of knowledge about chemical components:
built-in basic data about 350+ popular components,
the Chemical Component Dictionary (CCD) maintained by the PDB
(25,000+ components),
so-called CIF files compatible with the format of the CCP4 Monomer Library
(more about monomer libraries in the next section).
The built-in data is accessed through the function find_tabulated_residue.
It contains only minimal information about each residue:
assigned category, the “standard” flag (non-standard residues are marked
as HETATM in the PDB, even in polymer), one-letter code,
the number of hydrogens and molecular weight:
>>> gln=gemmi.find_tabulated_residue('GLN')>>> gln.is_amino_acid()True>>> gln.one_letter_code'Q'>>> round(gln.weight,3)146.144>>> gln.hydrogen_count10>>> gemmi.find_tabulated_residue('DOD').is_water()True>>> # PDB marks "non-standard" residues as HETATM.>>> # Pyrrolysine is standard - some microbes have it.>>> gemmi.find_tabulated_residue('PYL').is_standard()True>>> gemmi.find_tabulated_residue('MSE').is_standard()False
One-letter code is an upper case letter if it is a standard residue.
Otherwise, it can be the letter for the parent residue in lower case,
or a space. It is common to use X for non-standard residue –
for this we have helper function fasta_code():
To get more complete information, including atoms and bonds in the monomer,
we need to first read either the CCD
or a monomer library.
The CCD components.cif file describes all the monomers
(residues, ligands, solvent molecules) from the PDB entries.
Importantly, it contains information about bonds.
Note
The absence of bond information in mmCIF files from wwPDB is a
well-known problem.
This information is included in so-called
updated mmCIF files from PDBe,
as well as in BinaryCIF and mmJSON files.
Macromolecular refinement programs need to know more about monomers
than the CCD can tell: they need to know how to restrain the structure.
Therefore, they have their own dictionaries of monomers (a.k.a monomer
libraries), such as the CCP4 Monomer Library (for Refmac),
where each monomer is described by one cif file.
These libraries are often complemented by user’s own cif files.
Gemmi provides a ChemComp class that corresponds to a monomer
from either the CCD or a cif file.
#include<gemmi/read_cif.hpp> // for read_cif_gz#include<gemmi/chemcomp.hpp> // for ChemComp, make_chemcomp_from_blockgemmi::ChemCompmake_chemcomp(constchar*path){gemmi::cif::Documentdoc=gemmi::read_cif_gz(path);// assuming the component description is in the last block of the filereturngemmi::make_chemcomp_from_block(doc.blocks.back());}
Example: ATP
For phosphate motifs of the form R-O-PO3, this rule deprotonates
phosphoryl oxygens (AceDRG-style doubly deprotonated representation).
The alkoxy R-O part is not deprotonated by this rule.
Example: H1T
Oxygen must be neutral, have no hydrogen neighbors, and have exactly one
heavy single-bond neighbor.
This is not general alcohol deprotonation. It normalizes pre-existing
single-bond oxide-like oxygens that are already non-protonated in the graph.
Example: NE5
Shown as charge-separated nitro resonance (R-N+(=O)-O-).
The exact AceDRG representation (formal charges vs valence-only
representation) is under review.
prepare_chemcomp() is the full restraint-preparation pipeline. The order
below is important because later stages depend on graph/charge edits made
earlier.
The function now also has an explicit options object,
PrepareChemcompOptions, which makes the runtime controls part of the API
instead of relying only on environment variables. The legacy overload is still
available.
detect carborane branches before the standard path:
full carborane mode (early return): enabled when the component has at
least one non-hydrogen atom with 4+ boron neighbors, and all
non-hydrogen atoms are only B/C/metal. In this case,
apply_carborane_mode() is used, CCP4 types are assigned, and the
function returns immediately.
Example: 1KW.
otherwise continue with the standard path, but keep a carborane-seed flag
for possible mixed-mode post-processing.
Example: 9UK.
seed missing angles from existing bonds (unless no_angles=false);
run apply_chemical_adjustments();
run add_n_terminal_h3() (may add H3 and corresponding N-centered angles);
add_n_terminal_h3() adds the third proton only for matching
N-terminus-like motifs. sync_n_terminal_h3_angles() is run later (after
table fill) to keep the newly added H3 angles consistent with sibling
N-H/N-H2 restraints.
atom_stereo – maps atom names to pdbx_stereo_config strings used by
chirality generation,
no_angles – skips the initial angle-seeding pass,
strict_mode, compat_mode, trace_mode – tri-state overrides with
values PrepareOverride.Auto, PrepareOverride.Disable, or
PrepareOverride.Enable.
In Auto mode, the corresponding environment variable is still consulted
(GEMMI_ACE_STRICT, GEMMI_ACE_COMPAT, GEMMI_ACE_TRACE). Enable and
Disable override the environment explicitly for that call.
Enables extra validation checks during the pipeline. At present, this
mainly checks for invalid restraint references, duplicate atom ids in
planes, and NaN bond/angle targets at the final stage. In strict mode,
such conditions raise an exception instead of being tolerated or silently
cleaned up.
compat_mode
Enables AceDRG-like behavior in places where Gemmi now has a more
chemistry-first default.
Current effects:
peptide mode is gated by the canonical-AA-backbone test
(N, CA, C, O, OXT and associated bonding / H naming),
torsion candidate selection uses the older AceDRG-like selector instead
of Gemmi’s scored selector,
pyranose chair torsion rewriting is always applied, instead of only when
coordinates are incomplete.
With compat_mode disabled, Gemmi prefers chemically reasonable defaults
over matching AceDRG quirks.
trace_mode
Emits phase-by-phase summaries to stderr. The trace reports how many atoms,
bonds, angles, torsions, chiralities and planes were added or removed in
each pipeline phase, plus net charge changes. This is intended for
debugging and regression analysis rather than normal use.
Structural biologists routinely use prior knowledge about biomolecules
to augment the data obtained in an experiment.
This prior knowledge is what we know about preferred geometries in molecules
(distances between atoms, etc.). This knowledge is extracted primarily from
experimental small molecule databases (COD and CSD) and QM calculations.
One way to store that prior knowledge is in a so-called monomer library.
In addition to describing monomers (chemical components
from the previous section), the monomer library also describes links
between monomers and may contain various other data useful in
macromolecular refinement.
In Gemmi, data from a monomer library is stored in the class MonLib.
Currently, MonLib is modeled after and works only with
the CCP4 monomer library.
This library was introduced in the
early 2000s
to provide restraint templates for Refmac.
There are only two other popular MX refinement programs: PHENIX and BUSTER.
PHENIX provides geostd,
which was forked from CCP4 ML and is still quite similar.
In BUSTER the prior knowledge
is organized differently.
The restraints we use are similar to those used in molecular dynamics
(bond, angle, dihedral, improper dihedral and planarity restraints).
Originally, the monomer library was created because MD potentials
were deemed inadequate for refinement. Since then,
both the restraints in experimental structural biology and MD potentials
have improved, independently of each other. In recent years there have been
a few examples of using AMBER and OpenMM potentials for MX refinement.
Currently, there is no clear advantage to one approach over the other.
MonLib has the following member variables (std:: omitted for readability):
stringmonomer_dir – the top-level directory, in CCP4 it’s $CLIBD_MON.
map<string,ChemComp>monomers – chemical components read from
the monomer library. Usually, we read only the components present in a model.
Chemical components contains restraint templates
map<string,ChemLink>links – link descriptions. Each ChemLink
contains rules for determining to what links it is applicable,
restrain templates for restraining a link, and the names of modifications
to be applied to linked monomers.
map<string,ChemMod>modifications – each modification is a set of rules
for changing ChemComp. The rules can add, remove or modify
atoms (they often remove a hydrogen atom) and restraints.
map<string,ChemComp::Group>cc_groups – groups defined in the library
to classify monomers. Each monomer is assigned a group (such as peptide
or non-polymer). The groups are then used to match link templates against
the links between monomers.
EnerLibener_lib – data from $CLIBD_MON/ener_lib.cif.
read_monomer_lib() takes a monomer directory path, a list of residue names,
and an optional Logger argument.
Example – reading the CCP4 monomer library (Refmac dictionary):
>>> monlib_path=os.environ['CCP4']+'/lib/data/monomers'>>> # Usually, residue names from a model are obtained by calling:>>> #resnames = gemmi.Model.get_all_residue_names()>>> resnames=['LEU','VAL','HIS','SER','ASN','HOH']>>> monlib=gemmi.MonLib()>>> monlib.read_monomer_lib(monlib_path,resnames,logging=sys.stderr)True
EnerLib is the in-memory representation of selected data from
ener_lib.cif (the same schema used in CCP4 monomer-library and
AceDRG-tables layouts).
It is attached to MonLib as MonLib.ener_lib.
_lib_bond -> bond-type entries (atom_type_1, atom_type_2, bond order,
length, esd).
How it is used:
MonLib.read_monomer_lib() reads ener_lib.cif and populates ener_lib,
MonLib.find_ideal_distance() consults it to estimate ideal distances
between atoms, including metal-aware radius handling and fallback tiers
when exact type-pair entries are missing.
Input is typically a CCD-style component definition with atom and bond
information. In practice, gemmidrg must also normalize chemistry before
table lookup, because small differences in protonation state or local bond
annotation can move an atom into a different type bucket.
protonation/deprotonation rules aligned with AceDRG behavior where
feasible,
selected functional-group corrections used to stabilize atom typing and
avoid type drift from equivalent input representations,
graph cleanup steps that make downstream typing and fallback selection
deterministic.
These functional-group corrections are intentionally pragmatic rather than
fully general pKa modeling. Current examples include:
oxoacid normalization (for phosphate/sulfate-like motifs), including
deterministic O-H deprotonation in qualifying local patterns,
carboxylate normalization (side-chain and terminal forms), including
removal of acidic hydrogens where AceDRG would represent a deprotonated
carboxylate,
amine/guanidinium normalization, including AceDRG-style protonation and
hydrogen completion for selected terminal or resonance-stabilized motifs,
selected special cases such as PF6-like phosphorus environments and
metal-adjacent non-metal charge correction used by downstream typing.
Bond/angle assignment is driven by hierarchical matching. Conceptually, it
tries the most specific local environment first, then relaxes constraints
in controlled steps until a statistically supported target is found.
levels 9-11: hash/hybrid/ring summaries from HRS-style hierarchy.
The search does not always start at level 0. It computes a dynamic start
level from key availability, skipping levels that cannot match for the
current atom pair.
For angles, keys are centered on the middle atom hash plus sorted flank
hashes, ring-size/hybridization value key, and then progressively reduced
context (roots, neighbor summaries, atom types).
The lookup ladder is:
1D: full key including types (approx level 0),
2D: no type component (approx 1),
3D: no nb component (approx 2),
4D: roots-only beyond hash/value key (approx 3),
5D: hash + value key (approx 4),
6D: hash-only summary (approx 6).
When center/flank hashes are equal, an additional swapped-center pass is
used to reproduce AceDRG table-orientation behavior. There is also a
wildcard partial-hash fallback used only when the exact hash triple has
no loaded table entries.
Acceptance is thresholded by observation support (defaults are AceDRG-like;
for many paths the default minimum is 3 observations).
In particular:
bond level 0 uses observation count from the matched record,
bond levels 1-8 use the number of contributing entries in the aggregate,
bond levels 9-11 (HRS hierarchy) accept presence of compatible entries,
angle 1D-5D paths are accepted only when count threshold is met,
angle wildcard/HRS paths are used only after specific keyed paths fail.
If type-specific statistical matching fails:
bonds fall back through HRS/element-hybrid routes, then CCP4 ener_lib
compatibility lookup,
angles fall back through HRS and then geometry defaults based on center
hybridization/coordination (with dedicated metal and high-coordination
handling).
This ordering keeps as much chemistry context as possible before using broad
fallbacks. It improves robustness on unusual ligands while still favoring
AceDRG-like targets whenever data are available.
Ring handling is central to output quality and compatibility. It affects
both atom typing and the final restraint lookup.
This area received substantial tuning to match AceDRG conventions:
ring membership and aromaticity are propagated into environment labels,
fused systems are treated as connected ring networks, not isolated
independent rings,
shared atoms in fused systems can keep mixed labels such as [5a,6a],
these labels are used directly in AceDRG signatures and therefore
influence which bond/angle statistics are selected.
The current implementation follows the AceDRG workflow quite closely.
Ring membership is first detected from the molecular graph and then folded
back into per-atom labels:
the ring search is AceDRG-style and currently enumerates rings up to
6 atoms,
each detected ring contributes a size marker such as 5, 6,
5a or 6a to every atom that belongs to it,
if an atom belongs to multiple rings, the markers are combined and
sorted into a bracketed annotation such as [5a,6a],
repeated memberships of the same size are preserved, so a fused atom can
also end up with forms such as [5a,5a] or 3x6a.
These annotations are not cosmetic. They are injected directly into AceDRG
environment signatures for both the center atom and its neighbors. For
example, a plain C can become C[6a], and a neighboring atom description
inside a larger signature can likewise carry [5a,6a]. This means that
fused-ring context changes the COD class, which in turn changes which
statistical bond and angle records are considered compatible.
This is what is meant by treating fused systems as connected ring networks
rather than isolated rings. The code still detects individual rings, but it
does not stop there and type each atom as if it belonged to only one of them.
Instead, every atom accumulates the full set of ring memberships coming from
all rings that pass through it:
an atom on a simple monocycle typically gets a single label such as [6]
or [6a],
an atom shared by two fused aromatic 6-rings can carry [6a,6a],
an atom at a 5/6 fusion can carry [5a,6a],
if the same size occurs three times, the display is compacted to forms such
as 3x6a.
So for a fused bicyclic or polycyclic system, the shared atoms explicitly
remember that they sit at a junction, and the adjacent atoms “see” that
junction when their neighborhood signatures are built. In practice this means:
the center atom type is influenced by the whole local fused network,
not just by one arbitrarily chosen ring,
neighbor descriptors also retain the mixed fused annotation,
two atoms that are both formally aromatic carbons may still get different
AceDRG signatures if one sits in a simple phenyl-like ring and the other
sits at a fused 5/6 or 6/6 junction.
This is important because restraint lookup is keyed by these signatures.
Treating fused systems as disconnected monocycles would collapse distinct
environments together and would tend to select less specific bond and angle
statistics for bridge and fusion atoms.
From Python, ring aromaticity can be inspected directly:
For aromaticity assignment, the implementation uses AceDRG-style
electron counting with a Huckel (4n+2) rule like criterion on ring
pi-electron totals, but with several important details:
aromaticity annotations present in the input dictionary are not used for
this step; AceDRG-style ring aromaticity is recomputed from the molecular
graph and local atom properties,
a ring is only eligible if all of its atoms pass the AceDRG planarity
gate, which in practice means bonding_idx==2 for every atom in the
ring,
this “planarity” is not taken from Cartesian coordinates; it is the
AceDRG internal sp2-like classification derived from bond orders, charges,
and local graph environment,
the strict phase counts pi electrons using the AceDRG non-metal view of
the local environment,
the permissive phase retries with AceDRG’s broader bookkeeping used for
edge cases, while still preserving the same ring list and fused context,
there is also an AceDRG special case for porphyrin-like systems:
when exactly four planar 5-member rings each contain 4 C and 1 N,
they are marked aromatic even if the initial strict pi count fails.
The distinction between the two phases matters:
strict aromaticity is used while building the primary COD classes for
statistical lookup,
if the permissive pass changes any ring labels, the COD classes are
rebuilt for output typing,
however, the strict classification is still retained for the table-lookup
path so that statistics are matched in the same way as AceDRG.
Fallback selection also tries to preserve ring and aromatic context as long
as possible before dropping to broader generic classes. Small differences in
labels such as [6] vs [6a], or [5a] vs [5a,6a], can cascade into
different atom types and therefore different restraints, so matching
AceDRG behavior here is important for practical parity.
Some chemistries need dedicated logic beyond generic rules. For example,
carborane-like systems have specialized handling aimed at reproducing
AceDRG-like typing and restraint targets more closely.
A major goal of this implementation is practical compatibility with
AceDRG output and conventions (not only broad chemical plausibility).
In particular, substantial effort has been invested in matching
AceDRG-like behavior for atom typing, protonation logic, and restraint
selection in edge-case chemistries.
The produced restraint values are refinement targets (ideal values and
sigmas/esds). They are empirical/statistical restraints, not a QM
geometry optimization.
One key concept in restraint generation is the CCP4 “energy type”
(_chem_comp_atom.type_energy). This is a chemistry-aware atom class
used by monomer-library restraint tables.
Gemmi can assign these types directly, without loading AceDRG tables:
The assignment uses the current ChemComp graph as-is. If you want
normalization, hydrogen completion or peptide-specific adjustments first,
run prepare_chemcomp() instead.
These types are not elements. They encode local environment features
such as hydrogen count, local bonding pattern and ring/aromatic context.
The full list of non-element types from ener_lib.cif
(metals and other elements use their element symbol directly):
Carbon sp –
CSP (triple bond or two double bonds),
CSP1 (triple bond + 1H).
Carbon sp2 –
C (no H, carbonyl C),
C1 (1H),
C2 (2H),
CR15 (1H in 5-ring),
CR16 (1H in 6-ring),
CR5 (no H in 5-ring),
CR6 (no H in 6-ring),
CR1 (between two pyrrole units),
CR1H (CR1 + 1H),
CR55 (between 5-5 rings),
CR56 (between 5-6 rings),
CR66 (between 6-6 rings).
Nitrogen sp3, 4 bonds –
NT (no H),
NT1 (1H),
NT2 (2H),
NT3 (3H),
NT4 (4H).
Nitrogen sp3, 3 bonds –
N/N30 (no H),
N31 (1H),
N32 (2H),
N33 (3H).
Oxygen sp2 –
O (no charge, main-chain O),
OC (charged),
OP (charged, bonded to P),
OS (charged, bonded to S),
OB (charged, bonded to B).
Oxygen sp3 –
O2 (2 bonds),
OC2 (2 bonds, charged),
OH1 (alcohol),
OH2 (water).
Phosphorus –
P (4 bonds),
P1 (3 bonds).
Sulfur –
S (no H),
SH1 (1H),
S1 (1 double bond),
S2 (2 bonds),
S3 (3 bonds, sulphoxide).
Silicon –
SI (tetragonal),
SI1 (other).
Hydrogen –
H (generic),
plus specific subtypes: HCH, HCR, HNC1, HNC2,
HNH1, HNH2, HNR5, HNR6, HOH1, HOH2, HSH1.
Historically, energy types were also used to look up ideal bond lengths
and angles. In practice this turned out to be too complex for too little gain.
Currently, their main purpose is to provide hydrogen-bonding type,
van der Waals radii, and ionic radii – i.e. parameters for
non-bonding interactions.
AceDRG environment types (described below)
are far richer, but there are around 800 000 of them.
In principle they could also be used for non-bonding interaction analysis,
but this has not been done yet.
AceDRG environment types are explicit local-neighborhood signatures.
Examples from AceDRG tables include C(C)(H)3, N(CC), O(C)(H),
S(CC)(O)3, P(CC)3(O), and ring/aromatic-aware forms such as
C(C[6a]C[6a]2) and N(C[6a]C[6a]2).
More complex “full” signatures can be substantially richer, for example:
