10
$\begingroup$

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?

enter image description here

$\endgroup$

3 Answers 3

11
$\begingroup$

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')

bondlenhths

$\endgroup$
1
  • $\begingroup$ +1. That's a beautiful diagram! Thanks for putting that code together for our user who asked the question, and for the rest of the community! $\endgroup$ Dec 28, 2021 at 22:14
10
$\begingroup$

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: enter image description here

$\endgroup$
5
  • $\begingroup$ Oh, thanks. Do you know if the style of the representation can be changed? I have never used Jmol before. can I get a cleaner representation like the one picture I posted? $\endgroup$ Dec 27, 2021 at 21:27
  • $\begingroup$ Jmol is a very versatile visualization program. You can control (almost) everything. Your picture is in 2D representation and mine is in 3D, I didn't found how to plot in 2D in Jmol ;( $\endgroup$
    – Camps
    Dec 27, 2021 at 21:32
  • $\begingroup$ Thanks again. Scuse me one more question. Con you modify the font, or remove only the dotted line? Is it possible to change the atom representation to go from spheres to something else? $\endgroup$ Dec 27, 2021 at 22:20
  • $\begingroup$ @CarlosAndrésdelValle follow-up questions about JMol can be asked in separate questions, where the whole MMSE userbase (~3700 users) can see it. If Camps knows the answer, then perhaps you'll get an answer in this comment chain, but in general we prefer one question per post. Also, if you want more answers (other than just Jmol), it might be a good idea to wait a while before clicking the green checkmark. The green checkmark can give others the impression that you're no longer looking for more answers. $\endgroup$ Dec 27, 2021 at 22:55
  • $\begingroup$ @CarlosAndrésdelValle, yes, you can modify the molecule representation as you wish. Jmol has a lot of options! Go through the Jmol Wiki page and example pages to get more information. $\endgroup$
    – Camps
    Dec 28, 2021 at 16:19
7
$\begingroup$

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

2D molecule image, overlaid with bond lengths

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .