How to plot molecules with angles and bond lengths

Question

There are many programs out there where you pass an XYZ coordinate file and the program will give you a visualization of the molecule. I want to make a plot like the one in the picture, with bond lengths on top and a second one with bond angles. Again, plotting the molecule, calculating the bond angles and lengths is easy. The problem is to put them on the same image in the right place. I have tried many molecular visualization programs such as gdis, xcrysden, and Python libraries such as RDKit and sisl. None of these allow me to do what I want. I guess I could photoshop them on top but that will be very tedious. My molecule has around 70 angles and 45 bonds. Does someone know a good way to do this?

Answer 1

Only a half answer (bond lengths), but maybe interesting for people working with the RDKit.

The code was written in a Jupyter Noteboook.

Calculate the 3D shape of a molecule for the bond lengths and the 2D shape for the depiction.

Get the bond lengths and add them to the 2D structure.

from rdkit.Chem import AllChem as Chem
from rdkit.Chem import rdMolTransforms
from rdkit.Chem.Draw import rdMolDraw2D, rdDepictor, IPythonConsole
rdDepictor.SetPreferCoordGen(True)
from IPython.display import Image

smiles = 'Nc1cc2[nH]cc(C(=O)O)c(=O)c2cc1F'

# 3D for lengths
mol3d = Chem.AddHs(Chem.MolFromSmiles(smiles))
Chem.EmbedMolecule(mol3d, randomSeed=42)
Chem.MMFFOptimizeMolecule(mol3d)

# 2D for depiction
mol2d = Chem.AddHs(Chem.MolFromSmiles(smiles))
rdDepictor. Compute2DCoords(mol2d)

bonds = [(x.GetBeginAtomIdx(), x.GetEndAtomIdx()) for x in mol3d.GetBonds()]

conf = mol3d.GetConformer()

b_lengths = []

for b in bonds:
    bl = rdMolTransforms.GetBondLength(conf,b[0],b[1])
    b_lengths.append(str(round(bl,2)))

d = rdMolDraw2D.MolDraw2DCairo(600, 400)

for n in range(len(b_lengths)):
    mol2d.GetBondWithIdx(n).SetProp('bondNote', b_lengths[n])

d.DrawMolecule(mol2d)
d.FinishDrawing()
d.WriteDrawingText('bonds.png')
Image('bonds.png')

Answer 2

Using Jmol you can label all the bonds, angle, etc. but in 3D using the command measure:

Renders a measurement between the specified atoms. See also set (measure). Two general syntaxes are available. In the older syntax, a series of two to four atom numbers are given, and the appropriate measure (distance, angle, or dihedral angle) is then displayed. The newer, more general syntax is as follows:

measure RANGE <minValue> <maxValue> ALL|ALLCONNECTED|DELETE (\<atom expression>) (\<atom expression>)

Using this syntax one can specify a set of measurements to define all at once. Note that these sets are embedded atom expressions that must be enclosed in parentheses. If neither ALL nor ALLCONNECTED is present, only the first matching atom in the entire model set (all frames, so probably the first frame) is matched in each atom expression. When ALL or ALLCONNECTED is specified, all matching criteria in all frames are generated, thus allowing for "animated" measures.

The image below is an example:

Answer 3

I also manage to do the bond lengths in a slightly different way. I also used RDKit and my code is based on this discussion. I wanted to upload my answer earlier but the holidays are busy days.

First of all, as I was working with .mol files instead of the smiles I noticed that GetBondLength was giving me the wrong values. Then, I calculated the distances directly from the coordinates of the .mol file. To load the data I used

coords = []
new_coords = []
with open(file) as f:
    for _ in range(3):
        next(f)
    for line in f:
        coords.append(line.split())

N = int(coords[0][0])

for i in range(1,int(N)+1):
    new_coords.append([float(coords[i][0]),float(coords[i][1]),float(coords[i][2])])

c = np.asarray(new_coords)

Then, to add the lengths I used

for bnd in molecule.GetBonds():
        bl = np.linalg.norm(c[bnd.GetBeginAtomIdx()]- c[bnd.GetEndAtomIdx()])
        molecule.GetBondWithIdx(bnd.GetIdx()).SetProp('bondNote',f'{bl:.2f}')

After that, I used the style proposed in the discussion to draw the molecule. Also, because I had some Hydrogen bonds That weren't showing I created a Chem.RWMol molecule, added manually the bonds as bonds of type zero, and then applied the bond lengths to draw

molecule = Chem.MolFromMolFile(File, removeHs=False)
new_mol = Chem.RWMol(molecule)
new_mol.AddBond(8,25,Chem.rdchem.BondType('21'))
new_mol.AddBond(10,27,Chem.rdchem.BondType('21'))
molecule2 = Chem.Mol(new_mol)

Draw_with_bond_lengths(molecule2)

The final result is