Getting Started with the RDKit in C++¶
What is this?¶
This document is intended to provide an overview of how one can use the RDKit functionality from C++. Like the ‘Getting Started with the RDKit in Python’ it is not comprehensive and it’s not a manual. It is modelled very closely on the Python version, and most of the text will be similar if not identical. It is a work-in-progress, and will be added to over the coming months.
Building and Running C++ RDKit Programs¶
Unlike Python scripts, which, once the environment is setup up
correctly, can be run directly from a script or Jupyter notebook
session, the C++ programs must be compiled and linked before they can
be run. This creates an additional step, which varies with operating
system. Using the program cmake makes it easier, but you may need to
experiment and change settings to get this working. You will
need a reasonably modern C++ compiler. On linux systems this will most
likely be the GNU compiler gcc, although it could be Clang. On Macs,
the opposite is true, and on Windows machines it will probably be
Visual C++. As of the March 2018 release, RDKit uses modern C++,
which currently means up to C++17. You will probably need to tell
your compiler that you want to use the new features. On Clang and
gcc, use the flag -std=c++17
. To see how we build our examples
you can look at the “.azure-pipelines” directory in the git
repository for more platform specific build information.
C++ programs generally have quite a lot of extra verbiage in them
compared to, for example, python scripts, that is similar or identical
in all programs. This would produce extraneous clutter in the example
code in this document. Because of this, only the minimum code to
exemplify the point being made is shown in the text. Full programs for
all the examples (generally with multiple examples in the same
program) are given in the $RDBASE/Docs/Book/C++Examples
directory in
the distribution. The particular program will be mentioned in the text.
The functions of the RDKit system are declared in a large number of different header files spread across several directories in the system, and defined across a number of different libraries. The right headers will need to be included in the source code, and libraries linked to during linking. Whilst it’s possible to include all headers and libraries in all executables this will result in slower compile times, especially if you are doing static linking. When linking to the static (.a) libraries rather than the shared-object (.so) ones, the order the libraries appear in the linking list can be important. See the CMakeLists.txt file in C++Examples directory for a good order. In this case, the same library list is used for all examples, so some will be unnecessary for some of the programs. The first 3 programs don’t need the Depictor and SubstructMatch libraries, for instance, although on my Ubuntu 16.04 system, the RDGeometryLib appears to need to be included twice. Working out which libraries need to be linked to and in what order can involve a tedious amount of trial and error.
Should You Use C++ or Python?¶
There is no doubt that it is much easier get started with Python. If you follow the installation instructions, you will be able to start programming and using scripts straightaway. If all you are going to do is use scripts to do relatively simple things, essentially stitching RDKit function calls together, there should be little or no speed issues with using the Python interpreted language, as all the RDKit functions are compiled C++ and well optimised. However, if you are going to do more complicated things, using a lot of your own programming logic and only using the RDKit for peripheral things like I/O, SMARTS matching, preparing 2D images and the like, then it is likely that you will have good performance gains if you write in C++ and compile to a native executable. One reason for faster executables from C++ is that because the code is only compiled once, it can be worth the compiler spending more time optimising the code for speed. In Python, where the compilation is done each time at run-time, this overhead is less acceptable. Writing inefficient code is relatively easy in any language, but the C++ compiler can save you from yourself-at higher optimisation levels, it will re-arrange loops, factorise expressions etc., so that the final executable may be difficult to align with the original source code. For example (Huw!), the gcc will change sqrt(a) * sqrt(b) to sqrt(a*b) removing an expensive square root operation.
Another consideration is the completeness of the API. A lot of the higher level functionality in RDKit is developed in Python, and back-porting to C++ occurs on a demand-driven basis. There are therefore examples of quite useful functionality, such as computing the RMS differences between conformers (GetConformerRMSMatrix) that are not available in C++. Of course, if this affects you you can always implement the C++ version and submit a Pull Request. Indeed, the RMS calculation is on its way.
Memory Management¶
Memory leaks (where objects are created using new
but never
destroyed using delete
) are particularly insidious bugs that can be
difficult to track down
(Valgrind is your friend!) and may be a
surprise to people used to managed-memory languages like python and
java. In RDKit, many of the functions return pointers to molecules
and accept such pointers as arguments. Modern C++ provides tools such
as shared
(std::shared_ptr
) and single use (std::unique_ptr
) (collectively
‘smart’) pointers that can be very helpful in preventing memory leaks.
If an object is created and a pointer to it is stored in either of
these, then when the smart pointer goes out of scope, the object
is automatically deleted. Otherwise, the shared pointer is used
exactly as one would use an ordinary pointer. The unique_ptr is a bit
different in that it cannot be copied (you will get a compiler message
such as
/Users/david/Tmp/rdkit-my-fork/Docs/Book/C++Examples/example1.cpp:16:8: error: call to
implicitly-deleted copy constructor of 'std::__1::unique_ptr<RDKit::RWMol,
std::__1::default_delete<RDKit::RWMol> >'
auto m2 = m1;
^ ~~
/Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin/../include/c++/v1/memory:2494:3: note:
copy constructor is implicitly deleted because 'unique_ptr<RDKit::RWMol,
std::__1::default_delete<RDKit::RWMol> >' has a user-declared move constructor
unique_ptr(unique_ptr&& __u) noexcept
if you try to) although they can be move
d. There is
some anecdotal evidence that VisualStudio doesn’t cope with vectors of
vectors of unique_ptr
(std::vector<std::vector<unique_ptr<RWMol>>>
,
for example).
The Molecule Classes¶
Unlike in the Python libraries, in C++ there are two different
molecule classes, RDKit::ROMol
and RDKit::RWMol
. They are both
declared in GraphMol.h. ROMol (the Read-Only molecule) is used in
most instances. It can’t be edited. On those occasions where you will
need to edit the molecule, you’ll need to use the RWMol (Read-Write).
Reading and Writing Molecules¶
The majority of basic molecular functionality is found in the RDKit namespace, and only a small number of header files will need to be included to cover most use cases:
#include <GraphMol/GraphMol.h>
#include <GraphMol/FileParsers/MolSupplier.h>
#include <GraphMol/FileParsers/MolWriters.h>
Reading Single Molecules¶
Individual molecules can be constructed using a variety of approaches (example1):
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "Cc1ccccc1" );
std::string file_root = getenv( "RDBASE" );
file_root += "/Docs/Book";
std::string mol_file = file_root + "/data/input.mol";
std::shared_ptr<RDKit::ROMol> mol2( RDKit::MolFileToMol( mol_file ) );
std::cout << *mol2 << std::endl;
std::shared_ptr<RDKit::ROMol> mol3( RDKit::SmilesToMol( "Cc1cccc" ) );
All these return a pointer to an ROMol on success, or NULL on failure. Obviously, the object must be deleted when finished with to prevent memory leaks. In the example above, and henceforth in this document, the molecules, apart from mol1, are wrapped in shared pointers so that the objects are deleted as soon as the shared pointer goes out of scope.
There is also an overload of operator""
as a convenient shorthand
method of parsing a SMILES string into a molecule
(example1):
using namespace RDKit;
auto mol = "C[C@H](F)c1ccc(C#N)cc1"_smiles;"
Normally the keyword auto
means you don’t need to know the return
type of the function - the compiler works it out and away you go.
However, in this case it is useful to know that you get a
std::unique_ptr<RWMol>
, which tells you it can’t be copied, the
molecule can be edited, the memory won’t leak and you need to use the
->
operator to use functions on the RWMol object.
If the molecule can’t be sanitized after SMILES parsing, an
RDKit::MolSanitizeException
(derived
from std::exception
) is thrown, and an attempt is made to provide
sensible error messages (example1):
try {
std::shared_ptr<RDKit::ROMol> mol( RDKit::SmilesToMol( "CO(C)C" ) );
} catch( RDKit::MolSanitizeException &e ) {
// empty catch
}
displays something like [15:58:22] Explicit valence of atom # 1 O, 3, is greater than permitted
and (example1)
try {
std::shared_ptr<RDKit::ROMol> mol( RDKit::SmilesToMol( "c1cc1" ) );
} catch( RDKit::MolSanitizeException &e ) {
// empty catch
}
displays something like: [12:20:41] Can't kekulize mol
.
Reading sets of molecules¶
Groups of molecules are read using a Supplier (for example, an
RDKit::SDMolSupplier
or an RDKit::SmilesMolSupplier
)
(example2):
std::unqiue_ptr<RDKit::ROMol> mol;
std::string file_root = getenv( "RDBASE" );
file_root += "/Docs/Book";
std::string sdf_file = file_root + "/data/5ht3ligs.sdf";
bool takeOwnership = true;
RDKit::SDMolSupplier mol_supplier( sdf_file , takeOwnership );
while( !mol_supplier.atEnd() ) {
mol.reset( mol_supplier.next() );
std::cout << mol->getProp<std::string>( "_Name" ) << " has " << mol->getNumAtoms() << " atoms." << std::endl;
}
gives
mol-295 has 20 atoms.
mol-54 has 24 atoms.
mol-15 has 24 atoms.
mol-732 has 26 atoms.
The supplier can be treated as a random-access object (example2):
RDKit::SDMolSupplier mol_supplier( "data/5ht3ligs.sdf" , takeOwnership );
for( int i = int( mol_supplier.length() ) - 1 ; i >= 0 ; --i ) {
unique_ptr<RDKit::ROMol> mol( mol_supplier[i] );
if( mol ) {
std::cout << mol->getProp<std::string>( "_Name" ) << " has "
<< mol->getNumAtoms() << " atoms." << std::endl;
}
}
gives
mol-732 has 26 atoms.
mol-15 has 24 atoms.
mol-54 has 24 atoms.
mol-295 has 20 atoms.
A good practice is to test each molecule to see if it was correctly read before working with it (example2):
bool takeOwnership = true;
RDKit::SDMolSupplier *mol_supplier = new RDKit::SDMolSupplier( "data/5ht3ligs.sdf" , takeOwnership );
for( int i = int( mol_supplier->length() ) - 1 ; i >= 0 ; --i ) {
std::unique_ptr<RDKit::ROMol> mol( (*mol_supplier)[i] );
if( !mol ) {
continue;
}
std::cout << mol->getProp<std::string>( "_Name" ) << " has "
<< mol->getNumAtoms() << " atoms." << std::endl;
}
An alternative type of Supplier, the RDKit::ForwardMolSupplier
can be used to read from file-like objects. This allows the reading
of compressed files, using, for example, the boost::iostreams
objects (example2):
std::string file_root = getenv( "RDBASE" );
file_root += "/Docs/Book";
boost::iostreams::filtering_istream ins;
ins.push( boost::iostreams::gzip_decompressor() );
std::string comp_sdf_file = file_root + "/data/actives_5ht3.sdf.gz";
ins.push( boost::iostreams::file_source( comp_sdf_file ) );
// takeOwnership must be false for this, as we don't want the SDWriter trying
// to delete the boost::iostream
bool takeOwnership = false;
RDKit::ForwardSDMolSupplier forward_supplier( &ins , takeOwnership );
while( !forward_supplier.atEnd() ) {
mol.reset( forward_supplier.next() );
if( mol ) {
std::cout << mol->getProp<std::string>( "_Name" ) << " has " << mol->getNumAtoms() << " atoms." << std::endl;
}
}
Note that the forward suppliers cannot be used in random-access mode, and a compile-time error will result if you attempt to (example2):
error: no match for ‘operator[]’ (operand types are
‘RDKit::ForwardSDMolSupplier’ and ‘int’)
mol = forward_supplier[1];
Writing molecules¶
Single molecules can be converted to text using several functions
present in the RDKit
namespace.
For example, for SMILES (example3):
#include <GraphMol/SmilesParse/SmilesWrite.h>
.
.
std::shared_ptr<RDKit::ROMol> mol( RDKit::MolFromMolFile( "data/chiral.mol" ) );
std::cout << RDKit::MolToSmiles( *mol ) << std::endl;
gives
C[C@H](O)c1ccccc1
and (example3)
bool isomeric = false;
std::cout << RDKit::MolToSmiles( *mol , isomeric ) << std::endl;
produces
CC(O)c1ccccc1
where the isomeric
in the second function call specifies that isomeric
SMILES should not be produced.
Note that the SMILES produced is canonical, so the output should be
the same no matter how a particular molecule is input. For example
(example3)
std::shared_ptr<RDKit::ROMol> mol1( RDKit::SmilesToMol( "C1=CC=CN=C1" ) );
std::cout << RDKit::MolToSmiles( *mol1 ) << std::endl;
std::shared_ptr<RDKit::ROMol> mol2( RDKit::SmilesToMol( "c1cccnc1" ) );
std::cout << RDKit::MolToSmiles( *mol2 ) << std::endl;
std::shared_ptr<RDKit::ROMol> mol3( RDKit::SmilesToMol( "n1ccccc1" ) );
std::cout << RDKit::MolToSmiles( *mol3 ) << std::endl;
all produce c1ccncc1
as output.
If you’d like to have the Kekule form of the SMILES, you need to Kekulize an RWMol copy of the molecule, using the Kekulize function declared in MolOps.h (example3):
#include <GraphMol/MolOps.h>
.
.
std::shared_ptr<RDKit::RWMol> mol4( new RDKit::RWMol( *mol ) );
RDKit::MolOps::Kekulize( *mol4 );
std::cout << RDKit::MolToSmiles( *mol4, true, true ) << std::endl;
gives
C[C@H](O)C1=CC=CC=C1
Note: as of March 2017, the SMILES provided when one requests kekuleSmiles are not canonical. The limitation is not in the SMILES generation, but in the kekulization itself.
MDL Mol blocks are also available (example3):
std::shared_ptr<RDKit::ROMol> mol1( RDKit::SmilesToMol( "C1=CC=CN=C1" ) );
std::cout << RDKit::MolToMolBlock( *mol1 ) << std::endl;
gives
RDKit 2D
6 6 0 0 0 0 0 0 0 0999 V2000
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 N 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 2 2 0
2 3 1 0
3 4 2 0
4 5 1 0
5 6 2 0
6 1 1 0
M END
To include names in the mol blocks, set the molecule’s “_Name” property (example3):
mol1 = RDKit::SmilesToMol( "C1CCC1" );
mol1->setProp( "_Name" , "cyclobutane" );
std::cout << RDKit::MolToMolBlock( *mol1 ) << std::endl;
gives
cyclobutane
RDKit 2D
4 4 0 0 0 0 0 0 0 0999 V2000
1.0607 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.0000 -1.0607 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.0607 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 1.0607 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
3 4 1 0
4 1 1 0
M END
Note that setProp, which is a general function, can be called on an ROMol as well as an RWMol, which came as a surprise to me as I had assumed a read-only molecule would be less changeable than that.
In order for atom or bond stereochemistry to be recognised correctly by most software, it’s essential that the mol block have atomic coordinates. It’s also convenient for many reasons, such as drawing the molecules. Generating a mol block for a molecule that does not have coordinates will, by default, automatically cause coordinates to be generated. These are not, however, added to the molecule, but re-generated each time.
You can either include 2D coordinates (i.e. a depiction), using the function in the RDDepict namespace and declared in RDDepictor.h (example4):
#include <GraphMol/Depictor/RDDepictor.h>
.
.
std::shared_ptr<RDKit::ROMol> mol1( RDKit::SmilesToMol( "C1CCC1" ) );
RDDepict::compute2DCoords( *mol1 );
std::cout << RDKit::MolToMolBlock( *mol1 ) << std::endl;
gives
RDKit 2D
4 4 0 0 0 0 0 0 0 0999 V2000
1.0607 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.0000 -1.0607 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.0607 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0000 1.0607 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
3 4 1 0
4 1 1 0
M END
Or you can add 3D coordinates by embedding the molecule (example4):
#include <GraphMol/DistGeomHelpers/Embedder.h>
#include <GraphMol/ForceFieldHelpers/MMFF/MMFF.h>
.
.
std::shared_ptr<RDKit::ROMol> mol2( RDKit::SmilesToMol( "C1CCC1" ) );
mol2->setProp( "_Name" , "cyclobutane3D" );
RDKit::DGeomHelpers::EmbedMolecule( *mol2 );
RDKit::MMFF::MMFFOptimizeMolecule( *mol2 , 1000 , "MMFF94s" );
std::cout << RDKit::MolToMolBlock( *mol2 ) << std::endl;
gives
cyclobutane3D
RDKit 3D
4 4 0 0 0 0 0 0 0 0999 V2000
-0.8321 0.5405 -0.1981 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.3456 -0.8799 -0.2639 C 0 0 0 0 0 0 0 0 0 0 0 0
0.7190 -0.5613 0.7314 C 0 0 0 0 0 0 0 0 0 0 0 0
0.4587 0.9006 0.5008 C 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
3 4 1 0
4 1 1 0
M END
The optimization step isn’t necessary, but it substantially improves the quality of the conformation.
To get good 3D conformations, it’s almost always a good idea to add hydrogens to the molecule first (example4):
std::shared_ptr<RDKit::ROMol> mol3( RDKit::MolOps::addHs( *mol2 ) );
RDKit::MMFF::MMFFOptimizeMolecule( *mol3 , 1000 , "MMFF94s" );
std::shared_ptr<RDKit::RWMol> mol4( new RDKit::RWMol( *mol3 ) );
RDKit::MolOps::addHs( *mol4 );
Note that there are 2 overloaded versions of addHs. The first takes an ROMol and, because that can’t be edited, returns a pointer to a new ROMol with the result. If you use this version be careful not to leak memory by not deleting mol2 when you are finished with it. The second takes an RWMol which it is able to modify in place. With shared pointers, memory leaks can be avoided (example4):
std::shared_ptr<RDKit::ROMol> mol3sp( RDKit::MolOps::addHs( *mol2 ) );
mol3sp->setProp( "_Name" , "cyclobutaneSP" );
RDKit::MMFF::MMFFOptimizeMolecule( *mol3sp , 1000 , "MMFF94s" );
Once the optimisation is complete, the hydrogens can be removed again (example4):
std::shared_ptr<RDKit::ROMol> mol5( RDKit::MolOps::removeHs( *mol3 );
RDKit::MolOps::removeHs( *mol4 );
Again, there are two versions, one of which has an opportunity for a memory leak.
If you’d like write the molecules to file, use the normal C++ streams (example4):
#include <fstream>
.
.
std::ofstream ofs( "data/foo.mol" );
ofs << RDKit::MolToMolBlock( *mol5 );
Writing sets of molecules¶
Multiple molecules can be written to a file using an object of a
concrete subclass of the MolWriter
class (example5):
#include <GraphMol/FileParsers/MolWriters.h>
.
.
std::string file_root = getenv( "RDBASE" );
file_root += "/Docs/Book";
std::string sdf_file = file_root + "/data/5ht3ligs.sdf";
bool takeOwnership = true;
RDKit::SDMolSupplier mol_supplier( sdf_file , takeOwnership );
std::vector<std::shared_ptr<RDKit::ROMol>> mols;
while( !mol_supplier.atEnd() ) {
std::shared_ptr<RDKit::ROMol> mol( mol_supplier.next() );
if( mol ) {
mols.push_back( mol );
}
}
RDKit::PDBWriter pdb_writer( "data/5ht3ligs.pdb" );
for( std::size_t i = 0 , is = mols.size() ; i < is ; ++i ) {
pdb_writer.write( *mols[i] );
}
A MolWriter can also be initialised to a file-like object, so compressed files can be written or molecules can be written to a string in memory (example5):
#include <sstream>
.
.
std::ostringstream oss;
// takeOwnership must be false for this, as we don't want the SDWriter trying
// to delete the std::ostringstream.
takeOwnership = false;
boost::shared_ptr<RDKit::SDWriter> sdf_writer( new RDKit::SDWriter( &oss , takeOwnership ) );
for( auto mol: mols ) {
sdf_writer->write( *(*it) );
}
std::cout << oss.str() << std::endl;
Other available writers include SmilesWriter and TDTWriter (for those of you with an interest in historical Cheminformatics).
Working with Molecules¶
Looping over Atoms and Bonds¶
Since the adoption of c++11, looping over atoms and bonds has become extremely simple:
std::shared_ptr<RDKit::ROMol> mol( RDKit::SmilesToMol( "C1OC1" ) );
for(auto atom: mol->atoms()) {
std::cout << atom->getAtomicNum() << " ";
}
std::cout << std::endl;
gives
6 8 6
As well as this, there are a method that uses the fact that atoms and bonds can be selected by index number (example6):
for( unsigned int i = 0; i < mol->getNumAtoms() ; ++i ) {
const RDKit::Atom *atom = mol->getAtomWithIdx( i );
std::cout << atom->getAtomicNum() << " ";
}
Likewise with bonds (example6):
for(auto bond: mol->bonds()) {
std::cout << bond->getBondType() << " ";
}
std::cout << std::endl;
for( unsigned int i = 0 , is = mol->getNumBonds() ; i < is ; ++i ) {
const RDKit::Bond *bond = mol->getBondWithIdx( i );
std::cout << bond->getIsAromatic() << " ";
}
std::cout << std::endl;
gives
1 1 1
0 0 0
A bond can be specified by the atoms at its ends, with a NULL pointer being returned if there isn’t one (example6):
std::shared_ptr<RDKit::ROMol> mol2( RDKit::SmilesToMol( "C1OC1Cl" ) );
const RDKit::Bond *bond = mol2->getBondBetweenAtoms( 0 , 1 );
std::cout << bond->getBeginAtomIdx() << " to "
<< bond->getBeginAtomIdx() << " is "
<< bond->getBondType() << std::endl;
if( !mol2->getBondBetweenAtoms( 0 , 3 ) ) {
std::cout << "No bond between 0 and 3" << std::endl;
}
The neighbours of an atom can also be extracted (example6):
const RDKit::Atom *atom = mol2->getAtomWithIdx( 2 );
for(const auto &nbri: make_iterator_range(mol->getAtomBonds(atom))) {
const RDKit::Bond *bond = (*mol)[nbri];
unsigned int nbr_idx = bond->getOtherAtomIdx(atom->getIdx());
int nbr_atnum = bond->getOtherAtom(atom)->getAtomicNum();
std::cout << nbr_idx << " : " << nbr_atnum << std::endl;
}
gives
1 : 8
3 : 17
0 : 6
The auto' keyword here is hiding the fact that
getBondBetweenAtomsreturns the index of the bond, an
unsigned int`.
Ring Information¶
It is relatively easy to obtain ring information for atoms and bonds (example7):
#include <GraphMol/MolOps.h>
.
.
std::shared_ptr<RDKit::ROMol> mol( RDKit::SmilesToMol( "OC1C2C1CC2" ) );
if( !mol->getRingInfo()->isInitialized() ) {
RDKit::MolOps::findSSSR( *mol );
}
for( unsigned int i = 0; i < mol->getNumAtoms() ; ++i ) {
const RDKit::Atom *atom = mol->getAtomWithIdx( i );
std::cout << mol->getRingInfo()->numAtomRings( atom->getIdx() ) << " ";
}
std::cout << std::endl;
for( unsigned int i = 0; i < mol->getNumBonds() ; ++i ) {
const RDKit::Bond *bond = mol->getBondWithIdx( i );
std::cout << mol->getRingInfo()->numBondRings( bond->getIdx() ) << " ";
}
std::cout << std::endl;
gives
0 1 2 1 1 1
0 1 2 1 1 1 1
Obviously, findSSSR only needs to be called once for the molecule. If you only need to know whether the atom or bond is in a ring, just test whether or not the return value is zero (example7):
const RDKit::Bond *bond = mol->getBondWithIdx( 1 );
if( mol->getRingInfo()->numBondRings( bond->getIdx() )) {
std::cout << "Bond " << bond->getIdx() << " is in a ring" << std::endl;;
}
gives
Bond 1 is in a ring
Other information about presence in smallest rings can also be obtained from the RingInfo object of the molecule (example7):
std::cout << "Atom 2 is in ring of size 3 : "
<< mol->getRingInfo()->isAtomInRingOfSize( 2 , 3 ) << std::endl;
std::cout << "Atom 2 is in ring of size 4 : "
<< mol->getRingInfo()->isAtomInRingOfSize( 2 , 4 ) << std::endl;
std::cout << "Atom 2 is in ring of size 5 : "
<< mol->getRingInfo()->isAtomInRingOfSize( 2 , 5 ) << std::endl;
std::cout << "Bond 1 is in ring of size 3 : "
<< mol->getRingInfo()->isBondInRingOfSize( 1 , 3 ) << std::endl;
gives
Atom 2 is in ring of size 3 : 1
Atom 2 is in ring of size 4 : 1
Atom 2 is in ring of size 5 : 0
Bond 1 is in ring of size 3 : 1
More detail about the smallest set of smallest rings (SSSR) is available (example7):
RDKit::VECT_INT_VECT rings;
RDKit::MolOps::symmetrizeSSSR( *mol , rings );
std::cout << "Number of symmetric SSSR rings : " << rings.size() << std::endl;
for( auto ring: rings ) {
for( auto ringat: ring ) {
std::cout << ringat << " ";
}
std::cout << std::endl;
}
gives
Number of symmetric SSSR rings : 2
1 2 3
4 5 2 3
As the name suggests, this is a symmetrized SSSR; if you are
interested in the number of “true” SSSR, use the findSSSR
function
(example7):
std::cout << "Number of SSSR rings : " << RDKit::MolOps::findSSSR( *mol ) << std::endl;
gives
2
The distinction between symmetrized and non-symmetrized SSSR is discussed in more detail below in the section The SSSR Problem.
Modifying molecules¶
Normally molecules are stored in the RDKit with the hydrogen atoms
implicit (i.e. not explicitly present in the molecular graph). When
it is useful to have the hydrogens explicitly present, for example
when generating or optimizing the 3D geometry, the
RDKit::MolOps::addHs
function can be used
(example8).
std::shared_ptr<RDKit::ROMol> mol1( RDKit::SmilesToMol( "CCO" ) );
std::cout << "Number of atoms : " << mol1->getNumAtoms() << std::endl;
std::shared_ptr<RDKit::ROMol> mol2( RDKit::MolOps::addHs( *mol1 ) );
std::cout << "Number of atoms : " << mol2->getNumAtoms() << std::endl;
gives
Number of atoms : 3
Number of atoms : 9
Recall that there are two versions of RDKit::MolOps::addHs
, as
described above.
The Hs can be removed again using the RDKit::MolOps::RemoveHs
function, which again has two forms
(example8):
std::shared_ptr<RDKit::RWMol> mol3( new RDKit::RWMol( *mol2 ) );
RDKit::MolOps::removeHs( *mol3 );
std::cout << "Number of atoms : " << mol3->getNumAtoms() << std::endl;
which returns the atom count to 3.
RDKit molecules are usually stored with the bonds in aromatic rings
having aromatic bond types. This can be changed with the
RDKit::MolOps::Kekulize
function, which must be called with an RWMol
(example9):
std::shared_ptr<RDKit::RWMol> mol( new RDKit::RWMol( *RDKit::SmilesToMol( "c1ccccc1" ) ) );
std::cout << "Order : " << mol->getBondWithIdx( 0 )->getBondType() << std::endl;
std::cout << "Aromatic : " << mol->getBondWithIdx( 0 )->getIsAromatic() << std::endl;
RDKit::MolOps::Kekulize( *mol );
std::cout << "After default Kekulize : Order : " << mol->getBondWithIdx( 0 )->getBondType() << std::endl;
std::cout << "After default Kekulize : Aromatic : " << mol->getBondWithIdx( 0 )->getIsAromatic() << std::endl;
gives
Order : 12
Aromatic : 1
After default Kekulize : Order : 2
After default Kekulize : Aromatic : 0
The bond orders are defined as the enum BondType in Bond.h, and an aromatic bond currently has the value 12. Note that by default the Kekulize function clears the aromatic flags on the atoms and bonds. This is in contrast to the Python version of Kekulize, which preserves the flags by default. The behaviour can be forced explicitly (example9.cpp):
std::shared_ptr<RDKit::RWMol> mol1( new RDKit::RWMol( *RDKit::SmilesToMol( "c1ccccc1" ) ) );
RDKit::MolOps::Kekulize( *mol1 , false );
std::cout << "After Kekulize, markAtomsBonds false : Aromatic : " << mol1->getBondWithIdx( 0 )->getIsAromatic() << std::endl;
std::shared_ptr<RDKit::RWMol> mol2( new RDKit::RWMol( *RDKit::SmilesToMol( "c1ccccc1" ) ) );
RDKit::MolOps::Kekulize( *mol2 , true );
std::cout << "After Kekulize, markAtomsBonds true : Aromatic : " << mol2->getBondWithIdx( 0 )->getIsAromatic() << std::endl;
gives
After Kekulize, markAtomsBonds false : Aromatic : 1
After Kekulize, markAtomsBonds true : Aromatic : 0
Bonds can be restored to the aromatic bond type using the
RDKit::MolOps::sanitizeMol
function:
RDKit::MolOps::sanitizeMol( *mol );
std::cout << "Order : " << mol->getBondWithIdx( 0 )->getBondType() << std::endl;
std::cout << "Aromatic : " << mol->getBondWithIdx( 0 )->getIsAromatic() << std::endl;
gives
Order : 12
Aromatic : 1
once more.
Working with 2D molecules: Generating Depictions¶
The RDKit has 2 different methods for generating depictions (sets of 2D
coordinates) for molecules. They are part of the
RDDepict namespace, is accessed via the RDDepict::Compute2DCoords
function (example10.cpp):
#include <GraphMol/Depictor/RDDepictor.h>
.
.
std::shared_ptr<RDKit::RWMol>xs mol( new RDKit::RWMol( *RDKit::SmilesToMol( "c1nccc2n1ccc2" ) ) );
RDDepict::compute2DCoords( *mol );
The default method is the original RDKit algorithm. The alternative is CoordGen, kindly contributed by Schrodinger. The latter is reckoned to be slower but better at congested structures and macrocycles. CoordGen isn’t always built into the RDKit libraries, so it’s a reasonable idea to use it conditionally (example10.cpp):
#ifdef RDK_BUILD_COORDGEN_SUPPORT
RDDepict::preferCoordGen = true;
#else
std::cout << "CoordGen support not available" << std::endl;
#endif
The CMakeLists.txt with the examples
defines RDK_BUILD_COORDGEN_SUPPORT
.
The 2D conformation is constructed to minimize intramolecular clashes,
i.e. to maximize the clarity of the drawing. Unlike the Python
equivalent, the depiction is not placed in a canonical orientation by
default. This can be forced by passing true
as the third parameter
(example10.cpp):
RDDepict::compute2DCoords( *mol , nullptr , true );
By default, all existing conformations are removed when the 2D coordinates are created. This can be changed by passing false as a 4th parameter. The 2D coordinates are added as another conformation of the molecule so it’s a bit tricky combining them both in the same molecule, and probably best avoided.
The Python API has a convenience function
GenerateDepictionMatching2DStructure
which forces the 2D coordinate
generation to orientate molecules according to a template structure.
A C++ version of the function, generateDepictionMatching2DStructure
was included in late December 2016. If that is later than the version
of RDKit you are using, then the effect can be achieved thus:
(example10.cpp):
#include <Geometry/point.h>
#include <GraphMol/Substruct/SubstructMatch.h>
.
.
RDKit::ROMol *templ = RDKit::SmilesToMol( "c1nccc2n1ccc2" );
RDDepict::compute2DCoords( *templ );
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "c1cccc2ncn3cccc3c21" );
RDKit::MatchVectType matchVect;
if( RDKit::SubstructMatch( *mol1 , *templ , matchVect ) ) {
RDKit::Conformer &conf = templ->getConformer();
RDGeom::INT_POINT2D_MAP coordMap;
for(auto mv: matchVect) {
RDGeom::Point3D pt3 = conf.getAtomPos( mv.first );
RDGeom::Point2D pt2( pt3.x , pt3.y );
coordMap[mv.second] = pt2;
}
RDDepict::compute2DCoords( *mol1 , &coordMap );
}
Here, coordMap
maps the coordinates of atoms in the target
molecule templ onto corresponding atoms in the reference molecule.
It is also possible to produce a 2D picture that attempts to mimic as
closely as possible a 3D conformation. Again, an equivalent of the
Python function
rdkit.Chem.AllChem.GenerateDepictionMatching3DStructure
was
incorporated in December 2016.
Working with 3D Molecules¶
The RDKit can generate conformations for molecules using two different methods. The original method uses distance geometry [1]. The algorithm followed is:
The molecule’s distance bounds matrix is calculated based on the connection table and a set of rules.
The bounds matrix is smoothed using a triangle-bounds smoothing algorithm.
A random distance matrix that satisfies the bounds matrix is generated.
This distance matrix is embedded in 3D dimensions (producing coordinates for each atom).
The resulting coordinates are cleaned up somewhat using a crude force field and the bounds matrix.
Note that the conformations that result from this procedure tend to be fairly ugly. They should be cleaned up using a force field. This can be done within the RDKit using its implementation of the Universal Force Field UFF[2].
More recently, there is an implementation of the method of Riniker and Landrum [3] which uses torsion angle preferences from the Cambridge Structural Database (CSD) to correct the conformers after distance geometry has been used to generate them. With this method, there should be no need to use a minimisation step to clean up the structures; indeed, it is often undesirable as it may move the torsions away from the CSD-based distributions, somewhat negating the point.
The full process of embedding and optimizing a molecule is easier than all the above verbiage makes it sound (example11.cpp):
#include <GraphMol/DistGeomHelpers/Embedder.h>
#include <GraphMol/ForceFieldHelpers/UFF/UFF.h>
.
.
std::shared_ptr<RDKit::ROMol> mol( RDKit::SmilesToMol( "C1CCC1OC" ) );
std::shared_ptr<RDKit::ROMol>xs mol1( RDKit::MolOps::addHs( *mol ) );
// Original distance geometry embedding
RDKit::DGeomHelpers::EmbedMolecule( *mol1 , 0 , 1234 );
RDKit::UFF::UFFOptimizeMolecule( *mol1 );
// new Riniker and Landrum CSD-based method
// using the parameters class
RDKit::DGeomHelpers::EmbedParameters params( RDKit::DGeomHelpers::ETKDG );
params.randomSeed = 1234;
RDKit::DGeomHelpers::EmbedMolecule( *mol2 , params );
The Riniker and Landrum method has a number of parameters that may be altered for various reasons beyond the scope of this document. One that you may want to alter is the random number seed; setting the random number seed to other than the default -1 ensures that the same conformations are produced each time the code is run. This is convenient when testing to ensure reproducibility of results. To make it easier to vary the parameters, there is the EmbedParameters class which is initialised to the default values on construction, and whose individual values can be varied as desired.
The RDKit also has an implementation of the MMFF94 force field available [4], [5], [6], [7], [8]. (example11.cpp):
#include <GraphMol/ForceFieldHelpers/MMFF/MMFF.h>
.
.
RDKit::MMFF::MMFFOptimizeMolecule( *mol2 , 1000 , "MMFF94s" );
Please note that the MMFF atom typing code uses its own aromaticity model,
so the aromaticity flags of the molecule will be modified after calling
MMFF-related methods.
Note the calls to RDKit::MolOps::addHs()
in the examples above. By
default RDKit molecules do not have H atoms explicitly present in the
graph, but they are important for getting realistic geometries, so
they generally should be added. They can always be removed afterwards
if necessary with a call to RDKit::MolOps::removeHs()
With the RDKit, multiple conformers can also be generated using the two
different embedding methods. In both cases this is simply a matter of
running the distance geometry calculation multiple times from
different random start points. The 2nd parameter to
EmbedMultipleConfs
allows the user to
set the number of conformers that should be generated. Otherwise the
procedures are similar to before
(example11.cpp):
RDKit::INT_VECT mol1_cids = RDKit::DGeomHelpers::EmbedMultipleConfs( *mol1 , 10 );
std::cout << "Number of conformations : " << mol1_cids.size() << std::endl;
RDKit::INT_VECT mol2_cids;
int numConfs = 20;
RDKit::DGeomHelpers::EmbedMultipleConfs( *mol2 , mol2_cids , numConfs , params );
std::cout << "Number of conformations : " << mol2_cids.size()
<< std::endl;
The conformer ids are returned in mol1_cids
and mol2_cids
and
there are two overloaded functions with different ways of supplying
the information. As before, the CSD-based method is invoked by
EmbedParameters object, and in the example above the default
number of conformations to be produced has been changed from 10 to 20.
The conformers so generated can be aligned
to each other and the RMS values calculated
(example11.cpp):
#include <GraphMol/MolAlign/AlignMolecules.h>
.
.
std::vector<double> rms_list;
std::vector<unsigned int> m2cids( mol2_cids.begin() , mol2_cids.end() );
RDKit::MolAlign::alignMolConformers( *mol2 ,
static_cast<const std::vector<unsigned int> *>( 0 ) ,
&m2cids ,
static_cast<const RDNumeric::DoubleVector *>( 0 ) ,
false , 50 , &rms_list );
The RMS values for the overlays will be fed into rms_list on return.
Note the somewhat inconvenient issue that EmbedMultipleConfs
returns
a vector of ints
for the conformer ids, but alignMolConformers
requires a vector of unsigned ints
. The reason for this is that
EmbedMultipleConfs
uses -1 to denote a failed embedding. The first
vector of unsigned ints
in the alignMolConformers
declaration is atom ids, and allows
the alignment to be performed on just a subset of atoms which can be
convenient for overlaying a core and seeing how the other bits of the
molecule varied in the different conformations.
There is no C++ equivalent to the Python function
AllChem.GetConformerRMS()
to compute the RMS between two specific
conformers (e.g. 1 and 9) although it is coming.
It is important to remember that unless you specify a random number seed, you will not necessarily get the same conformations each time you run the embedding on the same molecule, especially if you only generate a small number of conformations relative to the number of torsions in the structure. If it’s important in your use-case that you have a good sampling of the conformations including all the low-energy ones, you should be sure to specify a large maximum number of conformations.
Disclaimer/Warning: Conformation generation is a difficult and subtle task. The original, default, 2D->3D conversion provided with the RDKit is not intended to be a replacement for a “real” conformational analysis tool; it merely provides quick 3D structures for cases when they are required. On the other hand, the second method, when a sufficiently large number of conformers are generated, should be adequate for most purposes. It is probably better to ignore the first, historical, method entirely. It is only left as the default method to avoid breaking existing code.
Preserving Molecules¶
Molecules can be preserved, or serialised, or pickled, using the class MolPickler in the namespace of the same name: (example12.cpp):
#include <GraphMol/MolPickler.h>
.
.
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "c1ccncc1" );
std::string pickle;
RDKit::MolPickler::pickleMol( *mol1 , pickle );
RDKit::ROMol mol2;
RDKit::MolPickler::molFromPickle( pickle , mol2 );
std::cout << RDKit::MolToSmiles( mol2 ) << std::endl;
Note that the string is in binary format and will appear as gibberish if printed to a screen. The RDKit pickle format is fairly compact and it is much, much faster to build a molecule from a pickle than from a Mol file or SMILES string, so storing pickles of molecules you will be working with repeatedly can be a good idea: (example12.cpp):
// writing to pickle file
std::string smi_file = getenv("RDBASE");
smi_file += "/Code/GraphMol/test_data/canonSmiles.long.smi";
std::string pkl_name = "canonSmiles.long.pkl";
// tab-delimited file, SMILES in column 0, name in 1, no title line
RDKit::SmilesMolSupplier suppl( smi_file , "\t" , 0 , 1 , false );
std::ofstream pickle_ostream( pkl_name , std::ios_base::binary );
int write_cnt = 0;
while( !suppl.atEnd() ) {
RDKit::ROMol *mol = suppl.next();
RDKit::MolPickler::pickleMol( *mol , pickle_ostream);
delete mol;
++write_cnt;
}
pickle_ostream.close();
std::cout << "Wrote " << write_cnt << " molecules" << std::endl;
// reading from pickle file
std::ifstream pickle_istream( pkl_name , std::ios_base::binary );
int read_cnt = 0;
while( !pickle_istream.eof() ) {
RDKit::ROMol mol3;
try {
RDKit::MolPickler::molFromPickle( pickle_istream , mol3 );
} catch( RDKit::MolPicklerException &e ) {
break;
}
++read_cnt;
}
pickle_istream.close();
std::cout << "Read " << read_cnt << " molecules." << std::endl;
By default, the pickling process does not preserve any properties attached to the molecule, which includes the molecule name (property “_Name”). This can be forced:
RDKit::MolPickler::pickleMol( *mol , pickle_ostream, RDKit::PicklerOps::AllProps );
Drawing Molecules¶
The RDKit has some built-in functionality for drawing molecules, found
in the RDKit namespace, with header files in
$RDBASE/Code/GraphMol/MolDraw2D
. There is an abstract base class
MolDraw2D which defines the interface and does the drawing, with
concrete classes for drawing to SVG or PNG files and Qt and wx
widgets. Only the SVG output is built by default, Cairo support
requires the argument -DRDK_BUILD_CAIRO_SUPPORT=ON
to cmake, and Qt
support -DRDK_BUILD_QT_SUPPORT=ON
. Note that the Qt code hasn’t been
looked at in years, is not tested as part of the release cycle and so
is not known to work. To create an SVG file:
(example13.cpp):
#include <GraphMol/MolDraw2D/MolDraw2DSVG.h>
.
.
RDKit::SDMolSupplier mol_supplier( "data/cdk2.sdf" , true );
RDKit::ROMol *mol1 = mol_supplier.next();
RDDepict::compute2DCoords( *mol1 );
std::ofstream outs("cdk_mol1.svg");
RDKit::MolDraw2DSVG svg_drawer(300, 300, outs);
svg_drawer.drawMolecule( *mol1 );
svg_drawer.finishDrawing();
outs.close();
The procedure for a PNG is slightly different: (example13.cpp):
#include <GraphMol/MolDraw2D/MolDraw2DSVG.h>
.
.
RDKit::MolDraw2DCairo cairo_drawer(300, 300);
cairo_drawer.drawMolecule(*mol1);
cairo_drawer.finishDrawing();
cairo_drawer.writeDrawingText("cdk_mol1.png");
As in Python, you can draw a set of molecules to a grid (example13.cpp):
std::string base_smi("c1ccccc1");
std::vector<std::string> extras = {"F", "Cl", "Br", "OC", "C(=O)O"};
std::vector<ROMol *> mols;
for(auto extra: extras) {
mols.push_back(SmilesToMol(base_smi + extra));
}
MolDraw2DSVG drawer(750, 400, 250, 200);
drawer.drawMolecules(mols);
drawer.finishDrawing();
std::ofstream grids(file_root + "/data/example_13_grid.svg");
grids << drawer.getDrawingText();
grids.flush();
grids.close();
Atoms in a molecule can be highlighted by drawing a coloured solid or open circle around them, and bonds likewise can have a coloured outline applied. An obvious use is to show atoms and bonds that have matched a substructure query (example17.cpp):
drawer.drawMolecule(*mol, &hit1_atoms, &hit1_bonds);
In the example above, hit1_atoms
and hit1_bonds
are vectors of
ints giving the atoms and bonds to be highlighted. They will all be
coloured the same default colour. It is possible to specify the
colours for individual atoms and bonds:
std::map<int, DrawColour> atom_cols;
atom_cols.insert(std::make_pair(13, DrawColour(1.0, 0.0, 0.0)));
atom_cols.insert(std::make_pair(14, DrawColour(0.0, 1.0, 0.0)));
atom_cols.insert(std::make_pair(15, DrawColour(0.0, 0.0, 1.0)));
atom_cols.insert(std::make_pair(16, DrawColour(1.0, 0.55, 0.0)));
std::map<int, DrawColour> bond_cols;
bond_cols.insert(std::make_pair(13, DrawColour(0.8, 0.8, 0.0)));
bond_cols.insert(std::make_pair(14, DrawColour(0.0, 1.0, 1.0)));
bond_cols.insert(std::make_pair(15, DrawColour(1.0, 0.0, 1.0)));
drawer.drawMolecule(*mol, &hit1_atoms, &hit1_bonds, &atom_cols, &bond_cols);
Atoms and bonds can also be highlighted with multiple colours if they fall into multiple sets, for example if they are matched by more than 1 substructure pattern. This is too complicated to show in this simple introduction, but there is an example at the bottom of example17.cpp.
As of version 2020.03, it is possible to add arbitrary small strings
to annotate atoms and bonds in the drawing. The strings are added as
properties common_properties::atomNote
and
common_properties::bondNote
and they will be placed automatically
close to the atom or bond in question in a manner intended to minimise
their clash with the rest of the drawing. There are 3 flags in
MolDraw2DOptions
that will add stereo information (R/S to atoms, E/Z
to bonds) and atom and bond sequence numbers
(example13.cpp):
auto m1 = "Cl[C@H](F)NC\\C=C\\C"_smiles;
MolDraw2DSVG drawer(250, 200);
m1->getAtomWithIdx(2)->setProp(common_properties::atomNote, "foo");
m1->getBondWithIdx(0)->setProp(common_properties::bondNote, "bar");
drawer.drawOptions().addAtomIndices = true;
drawer.drawOptions().addStereoAnnotation = true;
drawer.drawMolecule(*m1);
drawer.finishDrawing();
std::string text = drawer.getDrawingText();
std::ofstream outs(file_root + "/data/example_13_note.svg");
outs << text;
outs.flush();
Substructure Searching¶
Substructure matching can be done using query molecules built from SMARTS. (example14.cpp):
#include <GraphMol/Substruct/SubstructMatch.h>
.
.
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "c1ccccc1O" );
RDKit::RWMol *patt = RDKit::SmartsToMol( "ccO" );
RDKit::MatchVectType res;
if( RDKit::SubstructMatch( *mol1 , *patt , res ) ) {
std::cout << "Pattern matched molecule" << std::endl;
}
for( size_t i = 0 ; i < res.size() ; ++i ) {
std::cout << "(" << res[i].first << "," << res[i].second << ") ";
}
std::cout << std::endl;
SubstructMatch
returns a bool to flag whether there was a match, and
replaces the contents of res
with a mapping of the atom indices in
the pattern and a set of atoms that match in the molecule. In the
above example, the output is:
Pattern matched molecule
(0,0)(1,5)(2,6)
showing that atoms 0, 5 and 6 in the phenol matched the query. If the pattern matches multiple times (as in this case, where 4, 5, 6 is also a match), a single arbitrary set is returned.
All possible matches can also be returned: (example14.cpp):
std::vector<RDKit::MatchVectType> hits_vect;
if( RDKit::SubstructMatch( *mol1 , *patt , hits_vect ) ) {
for( size_t i = 0 ; i < hits_vect.size() ; ++i ) {
std::cout << "Match " << i + 1 << " : ";
for( size_t j = 0 ; j < hits_vect[i].size() ; ++j ) {
std::cout << "(" << hits_vect[i][j].first << ","
<< hits_vect[i][j].second << ")";
}
std::cout << std::endl;
}
}
This gives
Match 1 : (0,0)(1,5)(2,6)
Match 2 : (0,4)(1,5)(2,6)
It is easy to filter lists of molecules: (example14.cpp):
RDKit::SDMolSupplier mol_supplier( "data/actives_5ht3.sdf" , true );
RDKit::RWMol *patt1 = RDKit::SmartsToMol( "c[NH1]" );
std::vector<RDKit::ROMol *> matches;
while( !mol_supplier.atEnd() ) {
RDKit::ROMol *mol3 = mol_supplier.next();
if( mol3 && RDKit::SubstructMatch( *mol3 , *patt1 , res ) ) {
matches.push_back( mol3 );
} else {
delete mol3;
}
}
std::cout << "There were " << matches.size() << " hits in the file." << std::endl;
There should be 22 matches in the file.
Substructure matching can also be done using molecules built from SMILES instead of SMARTS: (example14.cpp):
RDKit::ROMol *mol4 = RDKit::SmilesToMol( "C1=CC=CC=C1OC" );
RDKit::RWMol *smi_mol1 = RDKit::SmilesToMol( "CO" );
if( RDKit::SubstructMatch( *mol4 , *smi_mol1 , res ) ) {
std::cout << "SMILES match" << std::endl;
} else {
std::cout << "Not SMILES match" << std::endl;
}
RDKit::RWMol *smt_mol1 = RDKit::SmartsToMol( "CO" );
if( RDKit::SubstructMatch( *mol4 , *smt_mol1 , res ) ) {
std::cout << "SMARTS match" << std::endl;
} else {
std::cout << "Not SMARTS match" << std::endl;
}
But don’t forget that the semantics of the two languages are not exactly equivalent: (example14.cpp):
RDKit::ROMol *mol4 = RDKit::SmilesToMol( "C1=CC=CC=C1OC" );
RDKit::RWMol *smi_mol2 = RDKit::SmilesToMol( "COC" );
if( RDKit::SubstructMatch( *mol4 , *smi_mol2 , res ) ) {
std::cout << "SMILES match" << std::endl;
} else {
std::cout << "Not SMILES match" << std::endl;
}
RDKit::RWMol *smt_mol2 = RDKit::SmartsToMol( "COC" );
if( RDKit::SubstructMatch( *mol4 , *smt_mol2 , res ) ) {
std::cout << "SMARTS match" << std::endl;
} else {
std::cout << "Not SMARTS match" << std::endl;
}
// Needs aromatic C
RDKit::RWMol *smt_mol3 = RDKit::SmartsToMol( "COc" );
if( RDKit::SubstructMatch( *mol4 , *smt_mol3 , res ) ) {
std::cout << "SMARTS match" << std::endl;
} else {
std::cout << "Not SMARTS match" << std::endl;
}
gives
SMILES match
Not SMARTS match
SMARTS match
Stereochemistry in substructure matches¶
By default, information about stereochemistry is not used in substructure searches: (example15.cpp):
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "CC[C@H](F)Cl" );
RDKit::RWMol *patt1 = RDKit::SmartsToMol( "C[C@H](F)Cl" );
RDKit::MatchVectType res;
if( RDKit::SubstructMatch( *mol1 , *patt1 , res ) ) {
std::cout << "SMARTS 1 match" << std::endl;
} else {
std::cout << "Not SMARTS 1 match" << std::endl;
}
RDKit::RWMol *patt2 = RDKit::SmartsToMol( "C[C@@H](F)Cl" );
if( RDKit::SubstructMatch( *mol1 , *patt2 , res ) ) {
std::cout << "SMARTS 2 match" << std::endl;
} else {
std::cout << "Not SMARTS 2 match" << std::endl;
}
RDKit::RWMol *patt3 = RDKit::SmartsToMol( "CC(F)Cl" );
if( RDKit::SubstructMatch( *mol1 , *patt3 , res ) ) {
std::cout << "SMARTS 3 match" << std::endl;
} else {
std::cout << "Not SMARTS 3 match" << std::endl;
}
All 3 SMARTS patterns match the molecule. To use the chirality
information, you need to pass true
as the optional fourth parameter,
corresponding to useChirality:
(example15.cpp):
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "CC[C@H](F)Cl" );
RDKit::RWMol *patt1 = RDKit::SmartsToMol( "C[C@H](F)Cl" );
RDKit::MatchVectType res;
if( RDKit::SubstructMatch( *mol1 , *patt1 , res , true , true ) ) {
std::cout << "SMARTS 1 chiral match" << std::endl;
} else {
std::cout << "Not SMARTS 1 chiral match" << std::endl;
}
RDKit::RWMol *patt2 = RDKit::SmartsToMol( "C[C@@H](F)Cl" );
if( RDKit::SubstructMatch( *mol1 , *patt2 , res , true , true ) ) {
std::cout << "SMARTS 2 chiral match" << std::endl;
} else {
std::cout << "Not SMARTS 2 chiral match" << std::endl;
}
RDKit::RWMol *patt3 = RDKit::SmartsToMol( "CC(F)Cl" );
if( RDKit::SubstructMatch( *mol1 , *patt3 , res , true , true ) ) {
std::cout << "SMARTS 3 chiral match" << std::endl;
} else {
std::cout << "Not SMARTS 3 chiral match" << std::endl;
}
gives
SMARTS 1 chiral match
Not SMARTS 2 chiral match
SMARTS 3 chiral match
Notice that when useChirality is true, a non-chiral query does match a chiral molecule. The same is not true for a chiral query and a non-chiral molecule: (example15.cpp):
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "CC[C@H](F)Cl" );
RDKit::RWMol *mol2 = RDKit::SmilesToMol( "CC(F)Cl" );
if( RDKit::SubstructMatch( *mol1 , *mol2 , res , true , true ) ) {
std::cout << "Chiral mol, non-chiral query : match" << std::endl;
} else {
std::cout << "Chiral mol, non-chiral query : NO match" << std::endl;
}
RDKit::RWMol *patt5 = RDKit::SmilesToMol( "C[C@H](F)Cl" );
if( RDKit::SubstructMatch( *mol2 , *patt5 , res , true , true ) ) {
std::cout << "Non-chiral mol, chiral query : match" << std::endl;
} else {
std::cout << "Non-chiral mol, chiral query : NO match" << std::endl;
}
gives
Chiral mol, non-chiral query : match
Non-chiral mol, chiral query : NO match
Atom Map Indices in SMARTS¶
It is possible to attach indices to the atoms in the SMARTS
pattern. This is most often done in reaction SMARTS (see Chemical
Reactions), but is more general than that. For example, in the SMARTS
patterns for torsion angle analysis published by Guba et
al. [9] indices are used to define the four atoms of the
torsion of interest. This allows additional atoms to be used to define
the environment of the four torsion atoms, as in
[cH0:1][c:2]([cH0])!@[CX3!r:3]=[NX2!r:4]
for an aromatic C=N
torsion. We might wonder in passing why they didn’t use recursive
SMARTS for this, which would have made life easier, but it is what it
is. The atom lists from GetSubstructureMatches are guaranteed to be in
order of the SMARTS, but in this case we’ll get five atoms so we need
a way of picking out, in the correct order, the four of interest. When
the SMARTS is parsed, the relevant atoms are assigned an atom map
number property that we can easily extract:
(example16.cpp):
RDKit::RWMol *patt1 = RDKit::SmartsToMol( "[cH0:1][c:2]([cH0])!@[CX3!r:3]=[NX2!r:4]" );
std::map<int,unsigned int> ind_map;
for(auto atom: patt1->atoms()) {
int map_num = atom->getAtomMapNum();
if( map_num ) {
ind_map[map_num-1] = atom->getIdx();
}
}
std::vector<unsigned int> map_list;
for(auto im: ind_map) {
map_list.push_back(im.second);
}
for( size_t i = 0 , is = map_list.size() ; i < is ; ++i ) {
std::cout << map_list[i] << " ";
}
std::cout << std::endl;
gives
0 1 3 4
Then, when using the query on a molecule, you can get the indices of the four matching atoms like this: (example16.cpp):
RDKit::ROMol *mol1 = RDKit::SmilesToMol( "Cc1cccc(C)c1C(C)=NC" );
std::vector<RDKit::MatchVectType> hits_vect;
if( RDKit::SubstructMatch( *mol1 , *patt1 , hits_vect ) ) {
for( size_t i = 0 ; i < hits_vect.size() ; ++i ) {
std::cout << "Match " << i + 1 << " : ";
for( size_t j = 0 ; j < map_list.size() ; ++j ) {
std::cout << hits_vect[i][map_list[j]].second << " ";
}
std::cout << std::endl;
}
}
gives
Match 1 : 1 7 8 10
The SSSR Problem¶
As others have ranted about with more energy and eloquence than I intend to, the definition of a molecule’s smallest set of smallest rings is not unique. In some high symmetry molecules, a “true” SSSR will give results that are unappealing. For example, the SSSR for cubane only contains 5 rings, even though there are “obviously” 6. This problem can be fixed by implementing a small (instead of smallest) set of smallest rings algorithm that returns symmetric results. This is the approach that we took with the RDKit.
Because it is sometimes useful to be able to count how many SSSR rings
are present in the molecule, there is a
rdkit.Chem.rdmolops.GetSSSR
function, but this only returns the
SSSR count, not the potentially non-unique set of rings.