18 - 3D Visualization#
This script demonstrates visualization of volumetric data from VASP:
Charge density (CHGCAR)
Electrostatic potential (LOCPOT)
Electron localization function (ELFCAR) python run.py
import os
import numpy as np
from ase.build import bulk
from vasp import Vasp
# Try to import matplotlib for plotting
try:
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
HAS_MATPLOTLIB = True
except ImportError:
HAS_MATPLOTLIB = False
print("Note: matplotlib not found. Plots will be skipped.")
# Try to import plotly for interactive 3D
try:
import plotly.graph_objects as go
HAS_PLOTLY = True
except ImportError:
HAS_PLOTLY = False
print("Note: plotly not found. Interactive 3D plots will be skipped.")
def read_vasp_volumetric(filename):
"""Read VASP volumetric data (CHGCAR, LOCPOT, ELFCAR format)."""
with open(filename) as f:
lines = f.readlines()
# Skip header (system name, scale, lattice, atoms)
# Find the grid dimensions line
idx = 0
for i, line in enumerate(lines):
parts = line.split()
if len(parts) == 3:
try:
ngx, ngy, ngz = int(parts[0]), int(parts[1]), int(parts[2])
idx = i + 1
break
except ValueError:
continue
# Read data
data = []
for line in lines[idx:]:
# Stop at augmentation charges or spin section
if 'augmentation' in line.lower():
break
values = line.split()
for v in values:
try:
data.append(float(v))
except ValueError:
break
if len(data) >= ngx * ngy * ngz:
break
# Reshape to 3D grid
data = np.array(data[:ngx * ngy * ngz])
data = data.reshape((ngx, ngy, ngz), order='F')
return data, (ngx, ngy, ngz)
def plot_2d_slice(data, axis=2, slice_idx=None, title='', filename=None):
"""Plot a 2D slice through volumetric data."""
if slice_idx is None:
slice_idx = data.shape[axis] // 2
if axis == 0:
slice_data = data[slice_idx, :, :]
xlabel, ylabel = 'y', 'z'
elif axis == 1:
slice_data = data[:, slice_idx, :]
xlabel, ylabel = 'x', 'z'
else:
slice_data = data[:, :, slice_idx]
xlabel, ylabel = 'x', 'y'
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.contourf(slice_data.T, levels=50, cmap='viridis')
plt.colorbar(im, ax=ax, label='results/Value')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
ax.set_aspect('equal')
if filename:
plt.savefig(filename, dpi=150, bbox_inches='tight')
print(f"Saved: {filename}")
plt.close()
def plot_isosurface_matplotlib(data, isovalue, title='', filename=None):
"""Plot isosurface using matplotlib (simple version)."""
from skimage import measure
# Use marching cubes to get isosurface
try:
verts, faces, _, _ = measure.marching_cubes(data, isovalue)
except Exception:
print(f"Could not generate isosurface at value {isovalue}")
return
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2],
cmap='viridis', alpha=0.7)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title(title)
if filename:
plt.savefig(filename, dpi=150, bbox_inches='tight')
print(f"Saved: {filename}")
plt.close()
def plot_isosurface_plotly(data, isovalue, title='', filename=None):
"""Plot isosurface using plotly (interactive)."""
X, Y, Z = np.mgrid[0:data.shape[0], 0:data.shape[1], 0:data.shape[2]]
fig = go.Figure(data=go.Isosurface(
x=X.flatten(),
y=Y.flatten(),
z=Z.flatten(),
value=data.flatten(),
isomin=isovalue * 0.9,
isomax=isovalue * 1.1,
surface_count=2,
colorscale='Viridis',
caps=dict(x_show=False, y_show=False, z_show=False),
opacity=0.6,
))
fig.update_layout(
title=title,
scene=dict(
xaxis_title='X',
yaxis_title='Y',
zaxis_title='Z',
aspectmode='data'
)
)
if filename:
fig.write_html(filename)
print(f"Saved: {filename}")
return fig
Main Calculation#
print("=" * 60)
print("3D Volumetric Data Visualization")
print("=" * 60)
print()
============================================================
3D Volumetric Data Visualization
============================================================
Step 1: Run VASP calculation with volumetric output#
print("Step 1: Running VASP calculation")
print()
# Use silicon as example
atoms = bulk('Si', 'diamond', a=5.43)
calc = Vasp(
label='results/viz/si',
atoms=atoms,
xc='PBE',
encut=300, # Lower for speed
kpts=(4, 4, 4),
ismear=0,
sigma=0.1,
lcharg=True, # Write CHGCAR
lvtot=True, # Write LOCPOT
lelf=True, # Write ELFCAR
prec='Normal',
)
energy = calc.potential_energy
print(f" Total energy: {energy:.6f} eV")
print()
Step 1: Running VASP calculation
Total energy: -10.835611 eV
Step 2: Read volumetric data#
print("Step 2: Reading volumetric data")
print()
calc_dir = calc.directory
# Check which files exist
files_to_read = ['CHGCAR', 'LOCPOT', 'ELFCAR']
available_files = {}
for fname in files_to_read:
fpath = os.path.join(calc_dir, fname)
if os.path.exists(fpath):
try:
data, grid = read_vasp_volumetric(fpath)
available_files[fname] = data
print(f" {fname}: grid = {grid}, range = [{data.min():.3f}, {data.max():.3f}]")
except Exception as e:
print(f" {fname}: could not read ({e})")
else:
print(f" {fname}: not found (file will be generated by VASP)")
print()
Step 2: Reading volumetric data
CHGCAR: grid = (32, 32, 32), range = [-7.503, 23.235]
LOCPOT: grid = (32, 32, 32), range = [-130.898, 4.806]
ELFCAR: grid = (16, 16, 16), range = [0.014, 0.954]
Step 3: Visualize charge density#
if HAS_MATPLOTLIB:
print("Step 3: Visualizing volumetric data")
print()
if 'CHGCAR' in available_files:
chg = available_files['CHGCAR']
# Normalize by volume
vol = atoms.cell.volume
chg_normalized = chg / vol
# 2D slices
plot_2d_slice(
chg_normalized,
axis=2,
title='Charge Density - XY plane (z=0.5)',
filename='charge_density_xy.png'
)
plot_2d_slice(
chg_normalized,
axis=0,
title='Charge Density - YZ plane (x=0.5)',
filename='charge_density_yz.png'
)
# Find appropriate isosurface value
mean_chg = chg_normalized.mean()
isovalue = mean_chg * 2 # Above average density
print(f" Charge density: mean = {mean_chg:.4f} e/ų")
print(f" Using isosurface value: {isovalue:.4f} e/ų")
# Try isosurface with scikit-image
try:
from skimage import measure
plot_isosurface_matplotlib(
chg_normalized,
isovalue,
title=f'Si Charge Density Isosurface (ρ = {isovalue:.3f} e/ų)',
filename='charge_density_3d.png'
)
except ImportError:
print(" Note: scikit-image not found, skipping 3D isosurface")
if 'LOCPOT' in available_files:
pot = available_files['LOCPOT']
plot_2d_slice(
pot,
axis=2,
title='Local Potential - XY plane (z=0.5)',
filename='potential_xy.png'
)
print(f" Local potential: range = [{pot.min():.2f}, {pot.max():.2f}] eV")
if 'ELFCAR' in available_files:
elf = available_files['ELFCAR']
plot_2d_slice(
elf,
axis=2,
title='ELF - XY plane (z=0.5)',
filename='elf_xy.png'
)
# ELF isosurface at high localization
try:
from skimage import measure
plot_isosurface_matplotlib(
elf,
isovalue=0.8,
title='Si ELF Isosurface (ELF = 0.8)',
filename='elf_3d.png'
)
except ImportError:
pass
print(f" ELF: range = [{elf.min():.3f}, {elf.max():.3f}]")
print()
Step 3: Visualizing volumetric data
Saved: charge_density_xy.png
Saved: charge_density_yz.png
Charge density: mean = 0.1999 e/ų
Using isosurface value: 0.3997 e/ų
Saved: charge_density_3d.png
Saved: potential_xy.png
Local potential: range = [-130.90, 4.81] eV
Saved: elf_xy.png
Saved: elf_3d.png
ELF: range = [0.014, 0.954]
Step 4: Interactive 3D with Plotly (if available)#
if HAS_PLOTLY and 'CHGCAR' in available_files:
print("Step 4: Creating interactive 3D visualization")
print()
chg = available_files['CHGCAR']
vol = atoms.cell.volume
chg_normalized = chg / vol
isovalue = chg_normalized.mean() * 2
fig = plot_isosurface_plotly(
chg_normalized,
isovalue,
title='Si Charge Density (interactive)',
filename='charge_density_interactive.html'
)
print()
Step 4: Creating interactive 3D visualization
Saved: charge_density_interactive.html
Summary#
print("=" * 60)
print("3D Visualization complete!")
print("=" * 60)
print()
print("Generated files:")
for ext in ['png', 'html']:
import glob
files = glob.glob(f'*.{ext}')
for f in files:
print(f" - {f}")
print()
print("Key points:")
print(" - CHGCAR: electron density, analyze bonding")
print(" - LOCPOT: electrostatic potential, work function")
print(" - ELFCAR: electron localization, lone pairs/bonds")
print(" - Use VESTA or ParaView for publication-quality rendering")
print()
print("Recommended tools for better visualization:")
print(" - VESTA: https://jp-minerals.org/vesta/")
print(" - ParaView: https://www.paraview.org/")
print(" - py3Dmol: pip install py3Dmol")
print()
print("Congratulations! You've completed the tutorial series.")
============================================================
3D Visualization complete!
============================================================
Generated files:
- charge_density_xy.png
- convergence_plots.png
- silicon_phonons.png
- charge_density_3d.png
- potential_xy.png
- equation_of_state.png
- silicon_bands.png
- elf_xy.png
- silicon_dos.png
- elf_3d.png
- charge_density_yz.png
- charge_density_interactive.html
Key points:
- CHGCAR: electron density, analyze bonding
- LOCPOT: electrostatic potential, work function
- ELFCAR: electron localization, lone pairs/bonds
- Use VESTA or ParaView for publication-quality rendering
Recommended tools for better visualization:
- VESTA: https://jp-minerals.org/vesta/
- ParaView: https://www.paraview.org/
- py3Dmol: pip install py3Dmol
Congratulations! You've completed the tutorial series.