This guide introduces methods for visualizing Bader charges, using a glycine dipeptide as a case study. The initial structure was created from a SMILES string and optimized to a linear backbone geometry.
from pfcc_extras.structure.ase_rdkit_converter import smiles_to_atoms
diglycine = smiles_to_atoms('NCC(=O)NCC(=O)O')
Table of Contents:
Preparation
Set up the calculator and retrieve the atomic charges.
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator from pfp_api_client.pfp.estimator import Estimator calc_mode, model_version = "PBE_PLUS_D3", "v8.0.0" estimator = Estimator(calc_mode=calc_mode, model_version=model_version) calculator = ASECalculator(estimator)
diglycine.calc = calculator
charge = diglycine.get_charges()
Visualization
(a) 3D Color Mapping of Atomic Charges
To visualize the distribution, we map the atomic charges to a color gradient.
First, we determine the scale by calculating abs_charge_max , which ensures the color map is symmetrically balanced between the maximum positive and negative values.
# Calculate the maximum absolute charge to define the color scale limits abs_max_charge = max(abs(charge))
We generate a color code ( c_hex ) based on these charges. We use the 'bwr' (blue-white-red) diverging colormap: blue for negative, white for zero, and red for positive charges. The intensity of the color corresponds to the magnitude of the charge.
For more colormap options, refer to the Matplotlib Colormap Documentation. https://matplotlib.org/stable/users/explain/colors/colormaps.html#diverging
from matplotlib import cm, colors import numpy as np cmap=cm.bwr_r # Normalize charge to the range [0, 1] for the colormap rgba = cmap((diglycine.get_charges()/abs_max_charge+1)/2) c_hex = [colors.rgb2hex(i) for i in rgba]
Now, prepare to display the molecule in the viewer (v).
from pfcc_extras import view_ngl v = view_ngl(diglycine) v
At this point, the molecule is rendered as follows:
By running the following code in a separate cell, the atomic coloring will be overlaid on the viewer. We use the Van der Waals radii to scale the atoms for better visibility.
from ase.data import vdw_radii
r_scale = 0.3 # Adjust vdW radius scale as preferred
radi = [float(vdw_radii[a.number]) * r_scale for a in diglycine]
# radi = np.full(num_atoms, 0.5) # Alternative: Display all atoms with the same size
v.view.clear_representations()
for i, color, radius in zip(diglycine, c_hex, radi):
v.view.add_representation("spacefill", selection=[i.index], color=color, radius=radius)
v.view.add_unitcell()
v.view.add_representation('licorice', radius=0.20, color='lightgray')
The viewer above will be updated as follows:
Finally, we will create a colorbar.
import matplotlib.pyplot as plt
import numpy as np
# Colorbar dimensions: Increase 4 for width, 0.15 for thickness
fig, ax = plt.subplots(figsize=(4, 0.15))
# Range and color settings
norm = plt.Normalize(vmin=-abs_max_charge, vmax=abs_max_charge)
sm = plt.cm.ScalarMappable(cmap='bwr', norm=norm)
# Set 5 tick marks
ticks = np.linspace(-abs_max_charge, abs_max_charge, 5)
cb = fig.colorbar(sm, cax=ax, orientation='horizontal', ticks=ticks)
# Tick and label settings
cb.set_ticklabels([f'{t:+.2f}' for t in ticks]) # Clean numerical labels
cb.set_label('Partial charge ($e$)', labelpad=5)
plt.show()
# plt.savefig('colorbar.png', dpi=300, bbox_inches='tight', transparent=True) # To save
By setting transparent=True in plt.savefig , you can create a colorbar with a transparent background as shown below. (Uncomment the line to save the file.)
Final Result:
Supplementary: To display atom indices and element species, run the following in a separate cell.
v.view.add_label(color="black",
labelType="text",
labelText=[diglycine[i].symbol + str(i) for i in range(diglycine.get_global_number_of_atoms())],
zOffset=1.5,
attachment='middle_center',
radius=1.0)
(b) 2D Structural Mapping with Charge Values
We can also use RDKit to display the charges on a 2D representation of the molecule.
The charge variable defined in the preparation step is used for this visualization.
from IPython.display import Image, display
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit import Chem
from rdkit.Chem import AllChem
# Generate 2D structure from SMILES
smiles = "NCC(=O)NCC(=O)O"
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol) # Add hydrogen atoms
AllChem.Compute2DCoords(mol) # Generate 2D coordinates
# Map charge values onto each atom as notes
for i, i_charge in enumerate(charge):
lbl = f"{i_charge:+.2f}"
mol.GetAtomWithIdx(i).SetProp("atomNote", lbl)
# Create a canvas
d2d = rdMolDraw2D.MolDraw2DCairo(600, 400)
opts = d2d.drawOptions()
# Optional: Display atom indices
# opts.addAtomIndices = False
# Reduce index font size (to distinguish them from charge values)
# opts.atomIndexFontSize = 10
# Draw the molecule
d2d.DrawMolecule(mol)
# To add a title:
# d2d.DrawMolecules([mol], legends=["Diglycine Bader Charge Distribution"])
d2d.FinishDrawing()
# Save the image
png_data = d2d.GetDrawingText() # Store drawing data in a variable
with open('diglycine_2d_w_bader-charge.png', 'wb') as f:
f.write(png_data)
# Display in the Jupyter notebook
display(Image(png_data))
(c) Quantitative Bar Plot Analysis
Plot the charges as a bar chart to quantitatively compare the values across all atoms.
The charge variable defined in the preparation step is used for this plot.
import matplotlib.pyplot as plt
import numpy as np
# 1. Prepare the data
labels = [f"{atom.symbol}{i}" for i, atom in enumerate(diglycine)]
indices = np.arange(len(charge))
# 2. Create the plot
fig, ax = plt.subplots(figsize=(8, 4))
ax.bar(indices, charge, color='#555555', width=0.7, zorder=3)
# 3. Axis and grid settings
ax.axhline(0, color='black', linewidth=1.0, zorder=4)
ax.grid(axis='y', linestyle='-', linewidth=0.5, color='#e0e0e0', zorder=1)
ax.set_ylabel('Bader Charge ($e$)', fontsize=10)
ax.set_xticks(indices)
ax.set_xticklabels(labels, rotation=45, fontsize=9)
# Refine appearance
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.tick_params(axis='y', direction='out', length=4)
ax.tick_params(axis='x', length=0)
# Display charge values on the bars
for i, c in enumerate(charge):
v_pos = c + 0.02 if c >= 0 else c - 0.02
ax.text(i, v_pos, f'{c:+.2f}',
ha='center',
va='bottom' if c >= 0 else 'top',
fontsize=8,
zorder=5, # Ensure text is on top
# Add white background to text for readability against grid lines
bbox=dict(facecolor='white', edgecolor='none', pad=0.5, alpha=0.8))
# Set Y-axis limits
y_limit = max(abs(charge)) * 1.2
ax.set_ylim(-y_limit, y_limit)
plt.tight_layout()
plt.savefig('bader_charge_bar-plot.png', dpi=300, bbox_inches='tight')
plt.show()