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?
3 Answers
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')
-
$\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
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
norALLCONNECTED
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. WhenALL
orALLCONNECTED
is specified, all matching criteria in all frames are generated, thus allowing for "animated" measures.
-
$\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
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