molecular-visualization-vmd

name	molecular-visualization-vmd
description	Geometry-first VMD pipeline for publication-quality molecular images: canonicalize the structure with PCA, generate candidate camera views on an azimuth/elevation/zoom grid, score them with objective visibility metrics, and render the best. Use when camera placement for a molecular scene must be deterministic rather than guessed.

Molecular Visualization with VMD

Overview

Geometry-first pipeline for generating publication-quality molecular images. Avoids LLM-based camera guessing by using deterministic geometry analysis.

Architecture

structure → canonical orientation → candidate cameras → score → render best

Stage A: Canonicalize Structure

• Align to principal axes (PCA)
• Center on object of interest
• Define domain axis (surface normal, bond axis, etc.)

Stage B: Generate Candidate Views

Sample a grid of views:

• Azimuth: 0°, 45°, 90°, ... 315°
• Elevation: -30°, 0°, 15°, 30°, 45°, 60°
• Zoom levels: 2.0, 2.5, 4.0 (closeup)

Stage C: Score with Objective Metrics

• Projected area of target motif
• Fraction of target atoms visible
• Depth separation between groups
• Reactive site visibility
• Molecular plane visibility

Stage D: Render Best Views

Use Tachyon ray tracer for high-quality output.

Environment Setup

# Create vmd-python environment
conda create -n vmd-env python=3.12
conda activate vmd-env
conda install -c conda-forge vmd-python pillow scikit-learn

Key Code Patterns

Load and Analyze Structure

import numpy as np
from sklearn.decomposition import PCA

def read_pdb_atoms(pdb_path):
    atoms = {'pt': [], 'c': [], 'n': [], 'h': []}
    atom_names = {'pt': [], 'c': [], 'n': [], 'h': []}

    with open(pdb_path) as f:
        for line in f:
            if line.startswith('ATOM'):
                name = line[12:16].strip()
                pos = np.array([float(line[30:38]),
                               float(line[38:46]),
                               float(line[46:54])])
                # classify by element prefix
                ...
    return atoms, atom_names

Compute Optimal Camera

def compute_optimal_camera(atoms):
    """Camera that views molecule from surface side."""
    mol_center = atoms['c'].mean(axis=0)
    pt_center = atoms['pt'].mean(axis=0)

    # View direction: from Pt center toward molecule
    view_dir = mol_center - pt_center
    view_dir = view_dir / np.linalg.norm(view_dir)

    # Convert to spherical (VMD uses z-forward)
    azimuth = np.degrees(np.arctan2(view_dir[0], view_dir[2]))
    elevation = np.degrees(np.arcsin(view_dir[1]))

    return azimuth, elevation

Find Reactive Hydrogens

def find_reactive_h(atoms, atom_names):
    """Find H atoms near nitrogen (reactive sites)."""
    n_pos = atoms['n'][0]

    # C atoms bonded to N (within 1.6 Å)
    c_near_n = [i for i, c in enumerate(atoms['c'])
                if np.linalg.norm(c - n_pos) < 1.6]

    # H atoms on those carbons (within 1.2 Å)
    reactive_h = []
    for i, h in enumerate(atoms['h']):
        for c_idx in c_near_n:
            if np.linalg.norm(h - atoms['c'][c_idx]) < 1.2:
                reactive_h.append(atom_names['h'][i])
                break

    return reactive_h[:4]  # typically 4 reactive H

VMD Setup and Rendering

import vmd
from vmd import molecule, molrep, evaltcl

def setup_molecule(pdb_path, atoms, atom_names, reactive_h):
    mol_id = molecule.load("pdb", str(pdb_path))

    # Clear default rep
    while molrep.num(mol_id) > 0:
        molrep.delrep(mol_id, 0)

    pt_sel = " ".join(sorted(set(atom_names['pt'])))

    # Pt - gray VDW
    molrep.addrep(mol_id, style="VDW 0.8 20", color="ColorID 2",
                  selection=f"name {pt_sel}", material="AOShiny")

    # C - dark gray CPK
    molrep.addrep(mol_id, style="CPK 1.0 0.3 20 20", color="ColorID 16",
                  selection=f"name {c_sel}", material="AOShiny")

    # N - blue VDW
    molrep.addrep(mol_id, style="VDW 0.6 20", color="ColorID 0",
                  selection=f"name {n_sel}", material="AOShiny")

    # Reactive H - RED (highlighted)
    molrep.addrep(mol_id, style="VDW 0.55 20", color="ColorID 1",
                  selection=f"name {reactive_h_sel}", material="AOShiny")

    return mol_id

def set_camera(azimuth, elevation, zoom):
    evaltcl("display resetview")
    evaltcl("display projection Orthographic")
    evaltcl("display depthcue off")
    evaltcl("color Display Background white")
    evaltcl("axes location Off")
    evaltcl(f"rotate y by {azimuth}")
    evaltcl(f"rotate x by {elevation}")
    evaltcl(f"scale by {zoom}")

def render_tachyon(output_path, width=1200, height=900):
    scene_path = output_path.with_suffix('.dat')
    evaltcl(f'render Tachyon "{scene_path}"')

    tga_path = output_path.with_suffix('.tga')
    subprocess.run(["tachyon", str(scene_path),
                   "-res", str(width), str(height),
                   "-format", "TGA", "-o", str(tga_path)])

    # Convert to PNG
    from PIL import Image
    img = Image.open(tga_path)
    img.save(output_path)

VMD Color IDs

ID	Color
0	Blue
1	Red
2	Gray
8	White
16	Black

Recommended View Types

Purpose	Azimuth	Elevation	Zoom
Orientation comparison	(computed)	15°	2.5
Plane view (tilt angle)	(computed)	0°	2.5
Closeup (reactive H)	(computed)	10°	4.0
Top-down	(computed)	60°	2.0

Troubleshooting

Molecule Hidden Behind Surface

• Compute optimal camera based on molecule position relative to Pt center
• Try multiple azimuth angles (0°, 90°, 180°, 270°)
• The molecule may be on the "back" of the nanoparticle

Edge-On Molecules Hard to See

• Use plane-view camera (perpendicular to molecular plane)
• Try different configs with similar orientation angles
• 80-85° orientations often render better than 89°

vmd-python Not Available

conda activate vmd-env
# or create: conda create -n vmd-env python=3.12 && conda install -c conda-forge vmd-python

Files

• Pipeline script: vmd_configs/render_smart.py
• Edge-on fix: vmd_configs/render_edge_on_fix.py
• Output: figures/renders/

Usage

cd /path/to/project/figures/vmd_configs
conda run -n vmd-env python render_smart.py

Outputs for your molecular system

Image	Description
orient_flat.png	Flat molecule (6°) on surface
orient_tilted.png	Tilted molecule (36°)
edge_edge_on_82_side.png	Edge-on molecule (82°)
plane_*.png	Side views showing tilt angle
closeup_*.png	Reactive H highlighted (4x zoom)
top_*.png	Top-down view

Key Principles

1. Geometry first, VLM second - Don't ask LLM to hallucinate camera angles
2. Sample candidate views - Grid search over azimuth/elevation
3. Score objectively - Projected area, visibility, depth separation
4. Reproducible - Save camera params in metadata.json
5. Highlight reactive sites - Red spheres for reactive H atoms