July 24, 2025
I was chatting with Pat Walters at the CADD GRC and he brought up point 7 from his famous 14 points LinkedIn post. The usual justification given for using mol2 files is that they are the only file format we have that allows including partial charges. This is a widely held opinion, but it’s not really true. It’s actually quite easy to store arbitrary atom (or bond) properties in SD files and the RDKit has a simple mechanism for doing this. This post is a short tutorial on how to use this functionality.
The handling of atomic and bond properties in SD files is described in more detail in the docs.
Rather than doing the usual thing and generating Gasteiger-Marsilli charges, here I’m going to use the DASH-props tree that we developed in the Riniker lab to rapidly generate approximate AM1BCC charges. DASH uses a hierarchical tree of substructures to compute high-quality estimates of AM1BCC charges for a molecule. Because it’s based on substructures, the method is fast and conformation independent. If want to learn more about the algorithm, take a look at the open-access paper.
If you want to run the code in this notebook, you’ll need to install DASH props. This requires adding some packages to your conda environment. There’s an environment file in the DASH-tree repo that you can use, but if you already have a working RDKit environment, I think this is sufficient:
% conda install -c conda-forge pytables tqdmand then pip installing the DASH tree package:
% python -m pip install git+https://github.com/rinikerlab/DASH-treeAn aside: the way these properties are stored in the SD file is very easy to parse (The format is actually the result of some discussions I had with a couple of cheminformatics software vendors back in 2019; I think I ended up being the only one to implement what we had discussed), so adding support for this format to other software or toolkits that can already read SD files wouldn’t be much of a lift. If you’re a developer and have questions, please just ask.
from rdkit import Chem
from rdkit.Chem import Draw
import rdkit
print(rdkit.__version__)
# import the required stuff from DASH:
from serenityff.charge.tree.dash_tree import DASHTree, TreeType
from serenityff.charge.data import dash_props_tree_path2025.03.4# Load the property tree.
# Note, that the files will be automatically downloaded the first time the tree is loaded from the ETHZ Research Collection.
tree = DASHTree(tree_folder_path=dash_props_tree_path, tree_type=TreeType.FULL)Loading DASH tree data
Loaded 122 trees and dataConstruct a sample molecule and compute the charges:
Set the charges as atomic properties:
# this is the one extra call you have to make in order to
# convert the atom properties to a molecule property which will
# be written to the SD file:
Chem.CreateAtomDoublePropertyList(esomep_h,'DASH_AM1BCC_CHARGES')
# write the SD data to a string so that we can look at it:
from io import StringIO
sio = StringIO()
with Chem.SDWriter(sio) as w:
w.write(esomep_h)
sdf = sio.getvalue()
print(sdf)
RDKit 2D
25 25 0 0 0 0 0 0 0 0999 V2000
3.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
1.5000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.7500 -1.2990 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.7500 -1.2990 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.5000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.7500 1.2990 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.7500 1.2990 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.0000 0.0000 0.0000 N 0 0 0 0 0 0 0 0 0 0 0 0
-1.5000 -2.5981 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.1155 -3.5040 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
-3.0000 -2.5981 0.0000 N 0 0 0 0 0 0 0 0 0 0 0 0
1.5000 -2.5981 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
3.0000 -2.5981 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0
0.7500 -3.8971 0.0000 N 0 0 0 0 0 0 0 0 0 0 0 0
4.5000 0.0000 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
3.0000 -1.5000 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
3.0000 1.5000 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.5000 2.5981 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
1.5000 2.5981 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.7500 1.2990 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.5923 -0.7860 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.7500 -3.8971 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-3.5923 -1.8121 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
1.5000 -5.1962 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.2268 -4.0171 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 2 0
3 4 1 0
4 5 2 0
5 6 1 0
6 7 2 0
5 8 1 0
4 9 1 0
9 10 2 0
9 11 1 0
3 12 1 0
12 13 2 0
12 14 1 0
7 2 1 0
1 15 1 0
1 16 1 0
1 17 1 0
6 18 1 0
7 19 1 0
8 20 1 0
8 21 1 0
11 22 1 0
11 23 1 0
14 24 1 0
14 25 1 0
M END
> <atom.dprop.DASH_AM1BCC_CHARGES> (1)
-0.057020346264131966 -0.05498106388730991 -0.1417576117737365 -0.24649839652124084 0.19183992335343653 -0.14084077400735293 -0.13325346021838766 -0.85373148869180315 0.69066109221699112
-0.61948158338060078 -0.63810034626413192 0.6751982937456793 -0.58508759253159914 -0.67224362594143339 0.054944129070337355 0.054944129070337355 0.054944129070337355 0.1595973470895479
0.1595973470895479 0.41622965373586801 0.41622965373586801 0.31563710234748005 0.31563710234748005 0.31876819330440859 0.31876819330440859
$$$$
You can see that the property atom.dprop.DASH_AM1BCC_CHARGES has been added to the SD string. Also that we’re writing way too many sig figs… it’s probably worth adding an option to control that.
When that SD string is parsed, the atom properties are extracted automatically and assigned to the atoms:
suppl = Chem.SDMolSupplier()
suppl.SetData(sdf,removeHs=False)
m = next(suppl)
for atom in m.GetAtoms():
print(atom.GetIdx(), atom.GetDoubleProp("DASH_AM1BCC_CHARGES"))0 -0.057020346264131966
1 -0.05498106388730991
2 -0.1417576117737365
3 -0.24649839652124084
4 0.19183992335343653
5 -0.14084077400735293
6 -0.13325346021838766
7 -0.8537314886918032
8 0.6906610922169911
9 -0.6194815833806008
10 -0.6381003462641319
11 0.6751982937456793
12 -0.5850875925315991
13 -0.6722436259414334
14 0.054944129070337355
15 0.054944129070337355
16 0.054944129070337355
17 0.1595973470895479
18 0.1595973470895479
19 0.416229653735868
20 0.416229653735868
21 0.31563710234748005
22 0.31563710234748005
23 0.3187681933044086
24 0.3187681933044086That’s it… hopefully it is useful