Part 4: Visualization and analysis in VMD

Part 4 of the Getting Started with Membrane MD series. Last updated: April 2026. Tested with VMD 1.9.4a57 on Apple Silicon and Python 3.9 with matplotlib 3.9.

After Part 3, you have a real MD trajectory: a starting structure, six equilibration segments, and however many production chunks you’ve run. Now we’ll make use of it: loading the trajectory in VMD, rendering high-quality images, computing the two most important bilayer metrics (area per lipid and bilayer thickness), and validating against published experimental values.

This is also the most rewarding part — there’s something genuinely satisfying about watching your simulated bilayer self-correct to its physically-correct equilibrium structure, and seeing the data match decades of careful experimental measurements.

Step 1: Open VMD and load the system

Launch VMD:

vmd

Three windows open: VMD Main (the control panel), OpenGL Display (the 3D rendering window), and VMD Console (text output).

To load the system, you’ll bring in the topology (.psf) and the coordinates/trajectory (.pdb or .dcd) — in that order.

Load the topology file:

In VMD Main, click File → New Molecule…. A dialog opens. Click Browse… and navigate to your simulations folder, then select step5_input.psf. The “Determine file type” dropdown should auto-detect as CHARMM, XPLOR, NAMD PSF. Click Load.

The OpenGL display will be empty at this point — you’ve loaded the connectivity but not yet any atomic positions.

Load the equilibration trajectory:

Without closing the Molecule File Browser, click Browse… again and select step6.1_equilibration.dcd. The “Load files for:” dropdown at the top should already say “0: step5_input.psf” — make sure you’re loading onto the existing molecule, not a new one. Click Load.

The OpenGL display lights up with your bilayer system. By default VMD plays through all frames, so you’ll see the trajectory animating.

Append the rest of the equilibration trajectories:

Repeat the Browse → Load procedure for step6.2_equilibration.dcd through step6.6_equilibration.dcd, in order. Each new DCD appends to the previous trajectory rather than replacing it.

Then append the production trajectories:

Repeat for each of step7_production_1.dcd, step7_production_2.dcd, etc. — however many chunks you’ve run.

When done, close the Molecule File Browser. The frame counter at the bottom of VMD Main shows the total frame count. For a complete equilibration + 10 ns production:

Equilibration: 75 frames (steps 6.1-6.3 at 25 frames each) + 150 frames (steps 6.4-6.6 at 50 frames each) = 225 frames
Production: 100 frames (10 chunks × 10 frames each = 100 frames at 50,000-step DCD frequency)
Plus 1 initial frame from the PSF load
Total: 326 frames

You can now scrub through the entire ~12 ns trajectory using the playback controls at the bottom of VMD Main.

Step 2: Set up clean visual representations

The default representation (Lines, all atoms) is fast but ugly. Let’s build up something clearer.

Open the representations dialog: Graphics → Representations…

You’ll see one default rep listed (Lines / Name / all). We’ll modify this and add more.

Modify the default rep to show only lipids:

With the default rep selected, change “Selected Atoms” from all to:

resname DOPC

Press Enter or click Apply. Now only the lipids are shown.

Switch to Licorice for nicer rendering:

In the Drawing Method dropdown, change “Lines” to Licorice. The bilayer becomes much more visually substantial — cylinders for bonds and spheres at junctions.

Add a second rep highlighting phosphate atoms:

Click Create Rep to clone the current rep. With the new rep selected:

Change Selected Atoms to name P, press Enter
Change Drawing Method to VDW

Big tan/gold spheres now appear at the upper and lower membrane interfaces, marking the headgroup positions of each lipid. The bilayer’s two leaflets become unmistakable.

Set background to white:

In VMD Main, click Graphics → Colors…. In the dialog: Categories → Display, Names → Background, Colors → 8 white. The OpenGL window inverts.

You should now have a clean image: white background, two leaflets of DOPC drawn in licorice, and big spheres marking the phosphate headgroups. This is essentially the canonical “lipid bilayer” image you see in textbooks.

Step 3: Render a high-quality image

The OpenGL display is fast but visually basic — no shadows, no proper lighting, no anti-aliased edges. For publication-quality images, VMD includes a built-in ray tracer called Tachyon.

Set the projection to orthographic (so molecules don’t appear distorted by perspective):

Display → Orthographic (toggle on)

Position the view how you want it:

Click and drag to rotate, right-click-drag (or option-drag) to pan, scroll to zoom. Press = (equals key) at any point to auto-fit the system in the window. A clean side-view of the bilayer (so you see edge-on) is typically what you want for a bilayer image.

Resize the OpenGL Display window as large as you can — Tachyon’s “in-memory” mode renders at the same pixel size as the OpenGL window. On a 34-inch monitor, you can easily get to 1500×1500 pixels or larger.

Render with Tachyon:

File → Render… opens the File Render Controls dialog.

In the dropdown labeled “Render the current scene using:”, select Tachyon (internal, in-memory rendering). There are two Tachyon options — pick the internal, in-memory one (the other one needs a standalone Tachyon binary that may not be installed).

In the Filename field, give the output a path:

~/Dropbox/path/to/dopc_test/dopc_bilayer.tga

If the Render Command field has been modified from the default, click Restore default. The command should be just opening the file with the default OS app after rendering.

Click Start Rendering. You’ll see “Rendering complete” within 30-90 seconds depending on system size.

The output is a .tga file (Targa format — old but lossless). Most modern image viewers can open it. To convert to PNG or JPEG:

sips -s format png ~/Dropbox/path/to/dopc_test/dopc_bilayer.tga \
     --out ~/Dropbox/path/to/dopc_test/dopc_bilayer.png

sips is built into macOS.

Step 4: Save your VMD session

If you’ve put effort into setting up nice representations, you’ll want to save the session so you don’t have to rebuild it next time.

File → Save Visualization State…

Save it next to your simulation files, e.g., ~/Dropbox/path/to/dopc_test/dopc_inspection.vmd.

To reload later:

vmd -e dopc_inspection.vmd

This opens VMD with all your representations and views restored — though you’ll still need to reload trajectory files.

Step 5: Compute area per lipid using the Tk Console

VMD’s Tk Console is a Tcl scripting interface where you can run analysis commands programmatically. It’s the most powerful interface VMD offers.

Open the Tk Console:

Extensions → Tk Console

A new window with a > prompt appears.

Compute APL frame by frame and write to a data file:

set nframes [molinfo top get numframes]
set fp [open "apl_vs_time.dat" w]
puts $fp "# frame  A(A)  B(A)  APL(A^2)"
for {set i 0} {$i < $nframes} {incr i} {
    molinfo top set frame $i
    set a [molinfo top get a]
    set b [molinfo top get b]
    set apl [expr ($a * $b) / 100.0]
    puts $fp "$i  $a  $b  $apl"
}
close $fp

What this does:

molinfo top get numframes — gets the total number of trajectory frames
The for loop iterates over each frame
molinfo top set frame $i — sets VMD’s “current frame” to the i-th frame
molinfo top get a and ... b — read the X and Y dimensions of the simulation box for that frame (NAMD writes them into the DCD)
expr ($a * $b) / 100.0 — computes box area divided by 100 (lipids per leaflet) = APL in Å²
The result is written to apl_vs_time.dat in the current directory

You can verify the output:

exec head apl_vs_time.dat

You should see a header line and then numerical data for each frame.

A practical note: paste this kind of multi-line script carefully into the Tk Console. Some terminal-style behaviors can swallow characters or stop after newlines. If a script doesn’t seem to do anything, try splitting it into smaller chunks and running each separately.

Step 6: Compute bilayer thickness frame by frame

Bilayer thickness (D_PP) is the distance between the average z-position of upper-leaflet phosphates and the average z-position of lower-leaflet phosphates.

A subtlety: with a flexible NPT box, the bilayer can drift slightly along z, so we want to compute the bilayer center each frame rather than assuming z = 0.

set nframes [molinfo top get numframes]
set fp [open "thickness_vs_time.dat" w]
puts $fp "# frame  thickness(A)"
for {set i 0} {$i < $nframes} {incr i} {
    molinfo top set frame $i
    set selDOPC [atomselect top "resname DOPC"]
    set zmem [lindex [measure center $selDOPC] 2]
    set selP_top [atomselect top "name P and z>$zmem"]
    set selP_bot [atomselect top "name P and z<$zmem"]
    set zP_top [lindex [measure center $selP_top] 2]
    set zP_bot [lindex [measure center $selP_bot] 2]
    set thickness [expr $zP_top - $zP_bot]
    puts $fp "$i  $thickness"
    $selDOPC delete
    $selP_top delete
    $selP_bot delete
}
close $fp

What this does:

For each frame, computes the geometric center of all DOPC atoms; the z-coordinate of that center is the bilayer midplane
Selects phosphates above the midplane (upper leaflet) and below the midplane (lower leaflet)
Computes the z-coordinate of the centroid of each selection
Subtracts to get the thickness

The $sel delete calls at the end of each iteration are important: VMD selections accumulate in memory unless explicitly deleted, and over hundreds of frames this can add up to noticeable memory bloat.

Verify:

exec head thickness_vs_time.dat

Step 7: Plot APL and thickness vs time

The data files are easy to plot in any plotting tool. We’ll use Python with matplotlib for flexibility.

In a regular Terminal (not the Tk Console), navigate to your simulation folder and create a plotting script:

cd ~/simulations/dopc_test
cat > plot_metrics.py << 'EOF'
#!/usr/bin/env python3
"""Plot APL and bilayer thickness over the simulation trajectory."""

import numpy as np
import matplotlib.pyplot as plt

apl_data = np.loadtxt("apl_vs_time.dat", comments="#")
thick_data = np.loadtxt("thickness_vs_time.dat", comments="#")

frames = apl_data[:, 0]
apl = apl_data[:, 3]
thick = thick_data[:, 1]

fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=True)

axes[0].plot(frames, apl, color="C0", lw=1.5)
axes[0].axhline(67.5, color="gray", linestyle=":", alpha=0.7,
                label="Experimental DOPC at 30°C (~67.5 Å²)")
axes[0].set_ylabel("Area per lipid (Å²)")
axes[0].legend(loc="best", frameon=False)
axes[0].grid(alpha=0.3)

axes[1].plot(frames, thick, color="C1", lw=1.5)
axes[1].axhspan(37, 39, color="gray", alpha=0.15,
                label="Experimental D_PP range")
axes[1].set_xlabel("Frame index")
axes[1].set_ylabel("Bilayer thickness D_PP (Å)")
axes[1].legend(loc="best", frameon=False)
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig("metrics_vs_time.png", dpi=200)
plt.show()
EOF

python3 plot_metrics.py

A matplotlib window will pop up showing two stacked plots — APL and thickness over time — with horizontal reference lines marking the experimental values for DOPC.

If matplotlib is not installed:

pip3 install matplotlib --user

then re-run the script.

What you should see, broadly:

APL plot: flat at 69.7 Å² for the first ~50 frames (steps 6.1 and 6.2 of equilibration, before the Langevin piston engages), then dropping during step 6.3 and settling around 65-67 Å² for the rest. After production starts, APL should fluctuate around its equilibrated value with small thermal noise.
Thickness plot: more excursive during the early restrained equilibration steps, then settling around 39-40 Å for the rest.

Both plots should show that the system has equilibrated — that is, after the early transient, the values fluctuate around stable means rather than drifting steadily. Drift suggests insufficient equilibration; stable fluctuations suggest you’re in a real equilibrium state.

Step 8: Validate against experiment

For DOPC at 30°C, published experimental values are:

Area per lipid: ~67.5 Å² (Kučerka et al. 2011, Biochim. Biophys. Acta 1808:2761)
Bilayer thickness D_PP: 37–39 Å (same source)

Compute averages over your equilibrated frames (typically the last half or last quarter of your trajectory):

set nframes [molinfo top get numframes]
set start_frame [expr $nframes / 2]

# APL
set apl_sum 0.0
set count 0
for {set i $start_frame} {$i < $nframes} {incr i} {
    molinfo top set frame $i
    set apl_sum [expr $apl_sum + [molinfo top get a] * [molinfo top get b] / 100.0]
    incr count
}
puts "Average APL: [expr $apl_sum / $count] A^2"

# Thickness
set thick_sum 0.0
set count 0
for {set i $start_frame} {$i < $nframes} {incr i} {
    molinfo top set frame $i
    set selDOPC [atomselect top "resname DOPC"]
    set zmem [lindex [measure center $selDOPC] 2]
    set selP_top [atomselect top "name P and z>$zmem"]
    set selP_bot [atomselect top "name P and z<$zmem"]
    set thickness [expr [lindex [measure center $selP_top] 2] - [lindex [measure center $selP_bot] 2]]
    set thick_sum [expr $thick_sum + $thickness]
    incr count
    $selDOPC delete
    $selP_top delete
    $selP_bot delete
}
puts "Average thickness: [expr $thick_sum / $count] A"

For our test bilayer (200 lipids, ~12 ns total trajectory), expect roughly:

APL: 65–67 Å² (slightly below the experimental value but well within typical CHARMM36 range)
Thickness: 39–40 Å (slightly above the experimental upper bound, again within typical CHARMM36 range)

The fact that both metrics fall within the experimental range — without you tuning anything — is a powerful validation. Decades of careful neutron-scattering and X-ray-scattering experiments have constrained DOPC’s structural parameters, and your simulation independently reproduces them.

Step 9 (advanced): The APL-thickness coupling

If both APL and thickness are properly equilibrated, you should see a negative correlation between them: when APL fluctuates upward (lipids spread out laterally), thickness tends to decrease (the bilayer thins), and vice versa. This is the bilayer equivalent of volume conservation — at constant lipid count, lipids must thicken when the area shrinks.

To test this, plot APL vs. thickness using only the well-equilibrated frames (skipping the early transient):

cat > plot_apl_vs_thickness.py << 'EOF'
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt

apl_data = np.loadtxt("apl_vs_time.dat", comments="#")
thick_data = np.loadtxt("thickness_vs_time.dat", comments="#")

# Use only frames after equilibration is complete
frame_start = 100
apl = apl_data[frame_start:, 3]
thick = thick_data[frame_start:, 1]

slope, intercept = np.polyfit(apl, thick, 1)
r = np.corrcoef(apl, thick)[0, 1]

fig, ax = plt.subplots(figsize=(6.5, 5.5))
sc = ax.scatter(apl, thick, c=np.arange(len(apl)), cmap="viridis", s=20,
                alpha=0.7, edgecolor="none")
fit_x = np.array([apl.min(), apl.max()])
ax.plot(fit_x, slope * fit_x + intercept, color="red", lw=1.5,
        label=f"slope = {slope:.2f} Å/Å²\nr = {r:.3f}")

plt.colorbar(sc, ax=ax, label="Frame index")
ax.set_xlabel("Area per lipid (Å²)")
ax.set_ylabel("Bilayer thickness D_PP (Å)")
ax.legend(loc="best", frameon=False)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("apl_vs_thickness.png", dpi=200)
plt.show()
EOF

python3 plot_apl_vs_thickness.py

For a well-equilibrated DOPC bilayer, you should see a slope around -0.25 to -0.35 Å/Å² with a Pearson correlation coefficient r ≈ -0.5 to -0.8.

What does the magnitude tell you? If lipid volume were perfectly conserved (a single lipid is a rigid block of fixed volume), the slope would be d(D_PP)/d(APL) = -D_PP/APL ≈ -0.6 Å/Å². The fact that the measured slope is shallower (around -0.25 to -0.35) means the bilayer is partially compressible in its thickness direction — as it stretches laterally, it doesn’t thin proportionally. Lipid tails have some flexibility in how they pack, and water can penetrate slightly into the headgroup region. Both effects buffer the thickness response. This compressibility is real and well-known; published values from CHARMM36 simulations of fluid-phase DOPC are in the same range.

What you have now

After completing this part, you should have:

A loaded VMD session showing your full trajectory, with publication-grade representations
High-quality rendered images of the bilayer
Time-series data for area per lipid and bilayer thickness
Equilibrated averages that validate your simulation against published experimental values
An understanding of the area-thickness coupling that emerges from physical lipid behavior

This completes the basic membrane MD workflow. You now have a complete loop: build a bilayer in CHARMM-GUI, simulate it in NAMD, visualize and analyze in VMD, and validate against experiment.

Where to go from here

This series covered a single-component DOPC bilayer — the simplest possible membrane MD project. From here, you might explore:

Mixed-lipid bilayers to study lipid-lipid interactions, raft-like behavior, etc.
Membrane proteins embedded in bilayers (CHARMM-GUI’s Bilayer Builder supports this)
Asymmetric bilayers with different leaflets, more biologically realistic
Different force fields (Slipids, AMBER Lipid, MARTINI for coarse-grained work)
Different MD engines (GROMACS or OpenMM both run faster than NAMD on Apple Silicon thanks to Metal GPU support, and CHARMM-GUI generates inputs for both)
Analyses beyond APL and thickness: lipid order parameters, headgroup orientation distributions, lipid lateral diffusion, water penetration profiles, etc.

There’s a vast literature on membrane MD methods, but the most practical next step is usually to read recent papers on the system you’re interested in and copy their methodology — most authors are quite explicit about force field choices, equilibration protocols, and analysis approaches.

Troubleshooting

molinfo top get a returns 0

The DCD didn’t carry box dimensions for the loaded frames, or VMD didn’t parse them. This sometimes happens with older PDB files. Workaround: manually set the box dimensions before computing APL:

molinfo top set a 83.4865259
molinfo top set b 83.4865259
molinfo top set c 75.0

(Adjust values to match your system, e.g., what’s in step5_assembly.str.)

This makes APL an estimate only — for a properly variable-cell analysis, the box dimensions need to be in the DCD. NAMD writes them with dcdUnitCell yes (which CHARMM-GUI sets by default).

Tk Console scripts stop midway through

Some terminals or VMD versions can swallow newlines on copy-paste. If a script seems incomplete, try splitting into smaller chunks and running each separately. Also, ending each command with a semicolon can help on some systems.

Tachyon rendering produces a tiny low-resolution image

The “internal, in-memory” Tachyon renders at whatever size your OpenGL Display window currently is. To get high resolution, resize the OpenGL window as large as you can before clicking Start Rendering. The standalone Tachyon (the other dropdown option) has a -res flag for explicit resolution but requires the standalone binary to be installed.

APL doesn’t equilibrate even after long simulation

If you see APL drifting steadily across hundreds of nanoseconds without leveling off, something deeper may be wrong — possibly the system was built with the wrong number of lipids per leaflet for the box size, or the force field isn’t behaving as expected. Cross-check the system: 100 DOPC at 67-70 Å²/lipid should give a box around 82-84 Å square. If your initial box was very different, the system has too much “room” or too little.

Negative correlation isn’t visible in APL vs. thickness scatter

Make sure you’re using only the equilibrated portion of the trajectory. The early frames are dominated by the transient drop from initial APL ~70 to equilibrated APL ~66, which produces a positive correlation (both decreasing together) that masks the equilibrium relationship. Skip the first ~50% of frames at minimum.

End of series. Comments, corrections, and feedback welcome — please reach out to the lab.

Previous: Part 3 — Running your first simulation in NAMD