Tutorial 11: Basic MD Potentials in LAMMPS

Objective: Understand how adding complexity to MD potentials enables prediction of different material properties through hands-on LAMMPS simulations.

Using Pre-Run Simulations¶

For this tutorial, we have pre-run all simulations in a shared class directory. You can either:

Use pre-run results (recommended for in-class work): Copy outputs and run analysis
Run simulations yourself (for practice later): Use the LAMMPS input files and geometry builders provided

Pre-run simulation location:

/anvil/projects/x-chm250117/class_examples/tutorial11_completed/

Directory structure:

/anvil/projects/x-chm250117/class_examples/tutorial11_completed/
├── part1_argon/
│   ├── argon_lj.in
│   ├── argon_lj.log
│   ├── argon_lj.lammpstrj
│   ├── argon_rdf.dat
│   └── argon_final.data
├── part2_nitrogen/
│   ├── lj_only/
│   │   ├── n2_lj_only.in
│   │   ├── n2_lj_only.lammpstrj
│   │   └── log.lammps
│   ├── n2_system.data
│   ├── nitrogen_bonded.in
│   ├── n2.lammpstrj
│   └── n2_rdf.dat
└── part3_ethane/
    ├── ethane.data
    ├── ethane_dihedral.in
    ├── ethane.log
    └── ethane.lammpstrj

File Organization¶

For your own work, create a well-organized directory structure in your scratch space:

cd $SCRATCH
mkdir tutorial11
cd tutorial11

# Create subdirectories for each part
mkdir part1_argon
mkdir part2_nitrogen
mkdir part3_ethane

File Paths: Running Jupyter from Your Home Directory

The Python scripts in this tutorial assume you are running Jupyter from the same directory as your simulation files. However, most students launch Jupyter from their home directory ($HOME).

If your notebook is not in the same folder as your simulation, you must add explicit paths to all file operations.

Option 1: Add a path variable at the top of your notebook

# Set this to the pre-run simulation directory
sim_dir = "/anvil/projects/x-chm250117/tutorials/tutorial11/part1_argon"

# Then use it for all file operations:
with open(f'{sim_dir}/argon_rdf.dat', 'r') as f:
    # ... read file ...

Option 2: Change directory at the start of your notebook

import os
os.chdir("/anvil/projects/x-chm250117/tutorials/tutorial11/part1_argon")

# Now all file operations use the simulation directory
# The rest of the script works as written

Quick Check: Run this cell to see where your notebook thinks it is:

import os
print("Current directory:", os.getcwd())
print("Files here:", os.listdir('.'))

Using ASE (Atomic Simulation Environment) for Structure Building¶

Throughout this tutorial, we use ASE to build molecular structures. ASE is the industry-standard Python library for atomistic simulations.

Professional Tool:

Used in research labs and industry worldwide
Works with multiple simulation codes (LAMMPS, VASP, Quantum Espresso, CP2K)
Skills transfer to other computational chemistry/physics projects

Easier and More Reliable:

Built-in molecule library with correct geometries
Automatic bond detection and topology generation
No manual tetrahedral geometry calculations
Fewer bugs than custom scripts

Quick ASE Tutorial:

from ase import Atoms
from ase.build import molecule
from ase.io import write
import numpy as np

# Build a single water molecule
water = molecule('H2O')
print(f"Water has {len(water)} atoms")
print(f"Positions:\n{water.positions}")

# Available molecules
# Common ones: 'H2O', 'N2', 'CH4', 'C2H6', 'CO2', 'NH3', etc.

Creating a simulation box:

# Build multiple molecules
n2_mol = molecule('N2')

# Create a box with replicated molecules
molecules = []
for i in range(10):
    mol = molecule('N2')
    mol.positions += np.random.rand(3) * 20  # 20 Å box
    molecules.append(mol)

# Combine into one system
system = molecules[0]
for mol in molecules[1:]:
    system += mol

# Set periodic boundary conditions
system.set_cell([20, 20, 20])
system.set_pbc(True)

Documentation: https://wiki.fysik.dtu.dk/ase/

ASE is already installed in your class conda environment.

Part 1: Pure Lennard-Jones Fluid¶

1.1 The System: Liquid Argon¶

We’ll simulate 500 argon atoms in a periodic box at 300 K and 1 atm.

What LJ Can Predict:

✅ Liquid-vapor phase diagram
✅ Radial distribution function (local structure)
✅ Diffusion coefficient
✅ Pressure and density

What LJ Cannot Predict:

❌ Anything requiring bonds (molecules don’t exist)
❌ Directional interactions (all atoms identical)

1.2 LAMMPS Input File¶

The LAMMPS input file argon_lj.in:

# Liquid Argon Simulation - Pure Lennard-Jones
# Tutorial 11 - Part 1

# ========================================
# 1. Initialization
# ========================================
units           real                    # Angstrom, kcal/mol, fs
atom_style      atomic                  # No bonds, just atoms
boundary        p p p                   # Periodic in all directions

# ========================================
# 2. Create Atoms
# ========================================
lattice         fcc 5.26                # FCC lattice, density ~ liquid Ar
region          box block 0 4 0 4 0 4   # 4x4x4 unit cells
create_box      1 box                   # 1 atom type
create_atoms    1 box                   # Fill box with type 1

mass            1 39.948                # Argon atomic mass (g/mol)

# ========================================
# 3. Lennard-Jones Potential
# ========================================
pair_style      lj/cut 10.0             # LJ with 10 Angstrom cutoff
pair_coeff      1 1 0.238 3.405         # epsilon=0.238 kcal/mol, sigma=3.405 A

# ========================================
# 4. Settings
# ========================================
neighbor        2.0 bin                 # Neighbor list skin distance
neigh_modify    every 1 delay 0 check yes

# ========================================
# 5. Equilibration (NVT)
# ========================================
velocity        all create 300.0 12345  # Initialize velocities at 300K
fix             1 all nvt temp 300.0 300.0 100.0  # Nose-Hoover thermostat

timestep        2.0                     # 2 fs timestep
thermo          1000                    # Print every 1000 steps
thermo_style    custom step temp pe ke etotal press density

run             50000                   # 100 ps equilibration

# ========================================
# 6. Production (NVT)
# ========================================
reset_timestep  0
dump            1 all custom 1000 argon_lj.lammpstrj id type x y z
dump_modify     1 sort id

# Compute radial distribution function
compute         rdf all rdf 100 1 1     # 100 bins, type 1 with type 1
fix             2 all ave/time 100 10 1000 c_rdf[*] file argon_rdf.dat mode vector

run             100000                  # 200 ps production

write_data      argon_final.data        # Save final configuration

1.3 Running the Simulation¶

#!/bin/bash
#SBATCH -A chm250117
#SBATCH -p shared
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -t 00:10:00
#SBATCH -J argon_lj

module load gcc/11.2.0 openmpi/4.0.6
module load lammps/20210310

lmp < argon_lj.in > argon_lj.log

Using pre-run results:

# Copy pre-run outputs to your working directory
cp /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part1_argon/argon_rdf.dat .
cp /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part1_argon/argon_lj.log .

Or analyze directly from the shared directory by setting the path in your notebook.

1.4 Analysis: Radial Distribution Function¶

After obtaining the simulation outputs, analyze the results in a Jupyter notebook.

In Anvil OnDemand:

Go to Interactive Apps → Jupyter Notebook
Create a new notebook: analysis.ipynb

Notebook Cell 1: Load and Plot RDF

import numpy as np
import matplotlib.pyplot as plt

# Set path to simulation outputs
sim_dir = "/anvil/projects/x-chm250117/class_examples/tutorial11_completed/part1_argon"

# Load RDF data - handle LAMMPS fix ave/time format
# File has multiple frames, we want the last one (final averaged RDF)
with open(f'{sim_dir}/argon_rdf.dat', 'r') as f:
    lines = f.readlines()

# Find the last data block
data_lines = []
i = len(lines) - 1

# Read backwards to find the last complete data block
while i >= 0:
    line = lines[i].strip()
    if line and not line.startswith('#'):
        parts = line.split()
        if len(parts) == 4:  # Data line: Row r g(r) coord
            data_lines.insert(0, line)
        elif len(parts) == 2:  # Timestep line
            break
    i -= 1

# Parse the data
data = np.array([list(map(float, line.split())) for line in data_lines])
r = data[:, 1]      # Distance (column 2)
g_r = data[:, 2]    # g(r) (column 3)

plt.figure(figsize=(10, 6))
plt.plot(r, g_r, 'b-', linewidth=2)
plt.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='Ideal gas')
plt.xlabel('Distance r (Å)', fontsize=12)
plt.ylabel('g(r)', fontsize=12)
plt.title('Radial Distribution Function: Liquid Argon at 300K', fontsize=13, fontweight='bold')
plt.xlim(0, 10)
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('argon_rdf.png', dpi=150)
plt.show()

# Find first peak position
first_peak_idx = np.argmax(g_r[:50])  # Search first 50 points
print(f"\nFirst peak at r = {r[first_peak_idx]:.2f} Å")
print(f"Peak height g(r) = {g_r[first_peak_idx]:.2f}")
print(f"\nExperimental first peak: ~3.8 Å")

Expected Result:

First peak at ~3.7-3.8 Å (nearest neighbor distance)
Peak height ~2.5-3.0 (liquid structure)
g(r) → 1 at large r (bulk average)

Part 2: Diatomic Molecules (N₂) — The Need for Bonds¶

2.1 Experiment A: N₂ with LJ Only (This Will Fail!)¶

Before adding bonds, let’s see what happens if we try to simulate N₂ with only Lennard-Jones interactions.

The LAMMPS input n2_lj_only.in:

# N2 with ONLY LJ - Will molecules stay together?
# Tutorial 11 - Part 2a (Demonstration of Failure)

units           real
atom_style      atomic                      # No bonds!
boundary        p p p

# Create two N atoms close together (as if bonded)
region          box block 0 30 0 30 0 30
create_box      1 box

# Place two N atoms at "bond distance" (1.1 Å apart)
create_atoms    1 single 15.0 15.0 15.0
create_atoms    1 single 16.1 15.0 15.0    # 1.1 Å separation

mass            1 14.007

# Only LJ potential (no bonds!)
pair_style      lj/cut 10.0
pair_coeff      1 1 0.0690 3.31             # N atom LJ

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

# Start at low temperature
velocity        all create 100.0 347923
fix             1 all nvt temp 100.0 100.0 100.0

timestep        1.0
thermo          100

# Track distance between the two atoms
compute         1 all property/atom id x y z
dump            1 all custom 10 n2_lj_only.lammpstrj id type x y z

thermo_style    custom step temp pe ke

run             5000                        # 5 ps

write_data      n2_lj_only_final.data

Using pre-run results:

# Path to pre-run LJ-only simulation
ls /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen/lj_only/

Analysis notebook: check_separation.ipynb

import numpy as np
import matplotlib.pyplot as plt

# Set path to pre-run simulation
sim_dir = "/anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen/lj_only"

# Parse trajectory to get atom positions
positions_atom1 = []
positions_atom2 = []
timesteps = []

with open(f'{sim_dir}/n2_lj_only.lammpstrj', 'r') as f:
    lines = f.readlines()

i = 0
while i < len(lines):
    if 'ITEM: TIMESTEP' in lines[i]:
        timestep = int(lines[i+1].strip())
        # Skip to atoms section
        while i < len(lines) and 'ITEM: ATOMS' not in lines[i]:
            i += 1
        i += 1  # Skip the ATOMS header
        
        # Read two atoms
        atom1_data = lines[i].split()
        atom2_data = lines[i+1].split()
        
        pos1 = np.array([float(atom1_data[2]), float(atom1_data[3]), float(atom1_data[4])])
        pos2 = np.array([float(atom2_data[2]), float(atom2_data[3]), float(atom2_data[4])])
        
        positions_atom1.append(pos1)
        positions_atom2.append(pos2)
        timesteps.append(timestep)
    i += 1

# Calculate distances
distances = [np.linalg.norm(p1 - p2) for p1, p2 in zip(positions_atom1, positions_atom2)]

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Distance vs time
ax1.plot(timesteps, distances, 'b-', linewidth=2)
ax1.axhline(y=1.1, color='r', linestyle='--', linewidth=2, label='Initial separation (1.1 Å)')
ax1.axhline(y=3.31, color='g', linestyle='--', linewidth=2, label='LJ σ (3.31 Å)')
ax1.set_xlabel('Timestep', fontsize=12)
ax1.set_ylabel('N-N Distance (Å)', fontsize=12)
ax1.set_title('N₂ with LJ Only: Atoms Fly Apart!', fontsize=13, fontweight='bold', color='red')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)

# LJ potential curve
r_plot = np.linspace(1.0, 10, 200)
epsilon = 0.0690  # kcal/mol
sigma = 3.31      # Angstrom
U_LJ = 4 * epsilon * ((sigma/r_plot)**12 - (sigma/r_plot)**6)

ax2.plot(r_plot, U_LJ, 'k-', linewidth=2.5, label='LJ Potential')
ax2.axvline(x=1.1, color='r', linestyle='--', linewidth=2, label='Initial distance (1.1 Å)')
ax2.axhline(y=0, color='gray', linestyle=':', alpha=0.5)
ax2.axvline(x=3.71, color='g', linestyle='--', alpha=0.7, label='LJ minimum (3.71 Å)')
ax2.set_xlabel('Distance (Å)', fontsize=12)
ax2.set_ylabel('Energy (kcal/mol)', fontsize=12)
ax2.set_title('Why It Failed: No Attractive Well at 1.1 Å', fontsize=13, fontweight='bold')
ax2.set_xlim(1, 8)
ax2.set_ylim(-0.1, 2)
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('n2_lj_failure.png', dpi=150)
plt.show()

print("\n" + "="*60)
print("RESULTS: LJ-ONLY FOR N₂")
print("="*60)
print(f"Initial distance: 1.1 Å (N≡N bond length)")
print(f"Final distance: {distances[-1]:.2f} Å")
print(f"\nLJ minimum at: {2**(1/6) * sigma:.2f} Å = {2**(1/6) * 3.31:.2f} Å")
print(f"Energy at 1.1 Å: {4 * epsilon * ((sigma/1.1)**12 - (sigma/1.1)**6):.2f} kcal/mol (REPULSIVE!)")
print("\n→ At 1.1 Å, atoms are INSIDE the LJ repulsive wall")
print("→ They push apart and separate")
print("→ LJ minimum (~3.7 Å) is too far for a covalent bond")
print("\n✗ LJ ALONE CANNOT HOLD N₂ TOGETHER")
print("="*60)

What You’ll See:

Atoms start at 1.1 Å (bond length)
They rapidly push apart to ~3-4 Å (LJ equilibrium)
The “molecule” dissociates!
Lesson: LJ well is at ~3.7 Å, not 1.1 Å. We need a bond potential.

2.2 Experiment B: N₂ with Bonds (This Works!)¶

Now let’s add a harmonic bond potential to keep the N atoms together.

2.3 What Bonds Add¶

With bonds, we can now predict:

✅ Bond vibration frequency (IR spectroscopy)
✅ Heat capacity (vibrational contribution)
✅ Molecular identity (N₂ stays intact)

Still Cannot Predict:

❌ Molecular shape (linear vs. bent) — needs angles
❌ Internal rotation — needs dihedrals

2.4 LAMMPS Input File¶

The LAMMPS input file nitrogen_bonded.in:

# Nitrogen Gas (N2) - LJ + Harmonic Bonds
# Tutorial 11 - Part 2

units           real
atom_style      full                    # Includes charges (even if zero)
boundary        p p p

# ========================================
# Create Atoms Manually
# ========================================
read_data       n2_system.data          # Read pre-built structure

# ========================================
# Force Field
# ========================================
# Non-bonded: LJ (even for bonded atoms, applies to inter-molecular)
pair_style      lj/cut 10.0
pair_coeff      1 1 0.0690 3.31         # N atom LJ parameters

# Bonded: Harmonic bond
bond_style      harmonic
bond_coeff      1 200.0 1.10            # k=200 kcal/mol/A^2, r0=1.10 A (N≡N)

# CRITICAL: Exclude 1-2 interactions from LJ
special_bonds   lj 0.0 0.0 0.5          # Scale 1-2, 1-3, 1-4 by 0, 0, 0.5

# ========================================
# Settings
# ========================================
neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

# ========================================
# Equilibration
# ========================================
velocity        all create 300.0 54321
fix             1 all nvt temp 300.0 300.0 100.0

timestep        1.0                     # 1 fs (smaller for bond vibrations)
thermo          1000
thermo_style    custom step temp pe ke etotal press

run             50000                   # 50 ps equilibration

# ========================================
# Production
# ========================================
reset_timestep  0

# Output trajectory for analysis in OVITO
dump            1 all custom 100 n2.lammpstrj id mol type x y z
dump_modify     1 sort id

run             100000                  # 100 ps production

2.5 Creating the Initial Structure¶

The data file n2_system.data is already provided in the shared directory:

ls /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen/n2_system.data

from ase.build import molecule
import numpy as np

print("="*60)
print("BUILDING N2 SYSTEM: ASE + LAMMPS File Format")
print("="*60)

# System parameters
n_molecules = 100
box_size = 35.0
min_dist = 3.5  # Minimum distance between molecule centers (Angstrom)

np.random.seed(42)  # For reproducibility

# Build N2 molecules with collision detection
molecules = []
centers = []

print(f"\nPlacing {n_molecules} N2 molecules with overlap checking...")

for i in range(n_molecules):
    placed = False
    attempts = 0
    max_attempts = 1000
    
    while not placed and attempts < max_attempts:
        # Random center position
        center = np.random.rand(3) * box_size
        
        # Check distance to all existing molecules
        too_close = False
        for existing_center in centers:
            # Minimum image distance (periodic)
            delta = center - existing_center
            delta = delta - box_size * np.round(delta / box_size)
            dist = np.linalg.norm(delta)
            if dist < min_dist:
                too_close = True
                break
        
        if not too_close:
            # Create and place molecule
            n2 = molecule('N2')
            
            # Random rotation
            n2.rotate(np.random.rand() * 360, 'x')
            n2.rotate(np.random.rand() * 360, 'y')
            n2.rotate(np.random.rand() * 360, 'z')
            
            # Center molecule at origin, then translate
            n2.positions -= n2.get_center_of_mass()
            n2.positions += center
            
            # Wrap positions into box
            n2.positions = n2.positions % box_size
            
            molecules.append(n2)
            centers.append(center)
            placed = True
        
        attempts += 1
    
    if not placed:
        print(f"Warning: Could not place molecule {i+1} after {max_attempts} attempts")

print(f"✓ Successfully placed {len(molecules)} molecules")

# Combine into one system
system = molecules[0].copy()
for mol in molecules[1:]:
    system += mol

system.set_cell([box_size, box_size, box_size])
system.set_pbc(True)

# Extract positions
positions = system.get_positions()

# Write LAMMPS data file
with open('n2_system.data', 'w') as f:
    f.write('N2 System Built with ASE\n\n')
    f.write(f'{len(positions)} atoms\n')
    f.write(f'{len(molecules)} bonds\n\n')
    
    f.write('1 atom types\n')
    f.write('1 bond types\n\n')
    
    f.write(f'0.0 {box_size} xlo xhi\n')
    f.write(f'0.0 {box_size} ylo yhi\n')
    f.write(f'0.0 {box_size} zlo zhi\n\n')
    
    f.write('Masses\n\n')
    f.write('1 14.007\n\n')
    
    f.write('Atoms  # full\n\n')
    for i in range(len(positions)):
        mol_id = (i // 2) + 1
        atom_id = i + 1
        pos = positions[i]
        f.write(f'{atom_id} {mol_id} 1 0.0 {pos[0]:.6f} {pos[1]:.6f} {pos[2]:.6f}\n')
    
    f.write('\nBonds\n\n')
    for i in range(len(molecules)):
        bond_id = i + 1
        atom1 = i * 2 + 1
        atom2 = i * 2 + 2
        f.write(f'{bond_id} 1 {atom1} {atom2}\n')

print(f"✓ Wrote n2_system.data")
print(f"✓ {len(positions)} atoms, {len(molecules)} bonds")
print(f"✓ Box: {box_size} × {box_size} × {box_size} Å³")

2.6 Analysis: Bond Coordination¶

Using pre-run results:

# Path to pre-run bonded N2 simulation
ls /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen/n2.lammpstrj

Notebook Cell 1: Instructions for OVITO Analysis

print("="*60)
print("N₂ BOND ANALYSIS USING OVITO")
print("="*60)
print("\nThe trajectory file 'n2.lammpstrj' is ready for analysis.")
print(f"\nFile location:")
print("  /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen/n2.lammpstrj")
print("\nTo analyze N-N bond distances in OVITO:")
print("\n1. Open OVITO (available on Anvil via Open OnDemand)")
print("   - Go to Interactive Apps → OVITO")
print("\n2. Load trajectory:")
print("   - File → Load File → n2.lammpstrj")
print("\n3. Add Coordination Analysis modifier:")
print("   - Add modification → Structure identification → Coordination analysis")
print("   - Set cutoff: 1.5 Å (to capture only bonded N-N pairs)")
print("\n4. Add Histogram modifier:")
print("   - Add modification → Visualization → Histogram")
print("   - Property: Coordination")
print("   - You should see a sharp peak at coordination = 1")
print("   - (Each N atom bonded to exactly 1 other N)")
print("\n5. Visual inspection:")
print("   - Particles → Display → Radius: 0.5 Å")
print("   - Particles → Bonds → Enabled")
print("   - Play animation to see molecules vibrating")
print("\n" + "="*60)
print("EXPECTED OBSERVATIONS:")
print("="*60)
print("✓ Each N atom has coordination number = 1 (bonded to 1 neighbor)")
print("✓ Bond lengths fluctuate around 1.10 Å due to thermal motion")
print("✓ Molecules rotate and translate but bonds don't break")
print("✓ At 300K, thermal energy (~0.6 kcal/mol) << bond strength (~200 kcal/mol)")
print("="*60)

Notebook Cell 2: Python Analysis from Trajectory

import numpy as np
import matplotlib.pyplot as plt

# Set path to pre-run simulation
sim_dir = "/anvil/projects/x-chm250117/class_examples/tutorial11_completed/part2_nitrogen"
n_molecules = 100

# Simple bond length calculator from trajectory
def calculate_bond_lengths(trajfile, n_molecules=100):
    """Calculate N-N bond lengths for each molecule at each timestep"""
    
    bond_lengths = []
    
    with open(trajfile, 'r') as f:
        lines = f.readlines()
    
    i = 0
    while i < len(lines):
        if 'ITEM: TIMESTEP' in lines[i]:
            i += 1
            timestep = int(lines[i].strip())
            
            # Skip to atoms section
            while i < len(lines) and 'ITEM: ATOMS' not in lines[i]:
                i += 1
            i += 1  # Skip ATOMS header
            
            # Read all atoms for this timestep
            positions = {}  # mol_id -> [pos1, pos2]
            
            for _ in range(n_molecules * 2):
                parts = lines[i].split()
                atom_id = int(parts[0])
                mol_id = int(parts[1])
                x, y, z = float(parts[3]), float(parts[4]), float(parts[5])
                
                if mol_id not in positions:
                    positions[mol_id] = []
                positions[mol_id].append(np.array([x, y, z]))
                i += 1
            
            # Calculate bond length for each molecule
            for mol_id in positions:
                if len(positions[mol_id]) == 2:
                    dist = np.linalg.norm(positions[mol_id][0] - positions[mol_id][1])
                    bond_lengths.append(dist)
        else:
            i += 1
    
    return np.array(bond_lengths)

# Calculate bond lengths
print("Calculating N-N bond lengths from trajectory...")
bond_lengths = calculate_bond_lengths(f'{sim_dir}/n2.lammpstrj')

# Statistics
mean_length = np.mean(bond_lengths)
std_length = np.std(bond_lengths)

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Time series (first 1000 bonds)
ax1.plot(bond_lengths[:1000], 'b-', linewidth=0.5, alpha=0.7)
ax1.axhline(y=1.10, color='r', linestyle='--', linewidth=2, label='Equilibrium (1.10 Å)')
ax1.axhline(y=mean_length, color='g', linestyle='-', linewidth=2, 
            label=f'Mean ({mean_length:.4f} Å)')
ax1.set_xlabel('Sample', fontsize=12)
ax1.set_ylabel('N-N Bond Length (Å)', fontsize=12)
ax1.set_title('Bond Length Fluctuations', fontsize=13, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# Distribution
ax2.hist(bond_lengths, bins=50, density=True, alpha=0.7, edgecolor='black')
ax2.axvline(x=1.10, color='r', linestyle='--', linewidth=2, label='r₀ = 1.10 Å')
ax2.axvline(x=mean_length, color='g', linestyle='-', linewidth=2,
            label=f'Mean = {mean_length:.4f} Å')
ax2.set_xlabel('Bond Length (Å)', fontsize=12)
ax2.set_ylabel('Probability Density', fontsize=12)
ax2.set_title('Bond Length Distribution', fontsize=13, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('n2_bond_analysis.png', dpi=150)
plt.show()

print(f"\nBond Statistics:")
print(f"  Mean length: {mean_length:.4f} Å")
print(f"  Std dev: {std_length:.4f} Å")
print(f"  Equilibrium: 1.10 Å")
print(f"\n→ Bonds vibrate around equilibrium")
print(f"→ Standard deviation related to temperature")

Key Observation: The bond vibrates around 1.10 Å but never breaks. At 300 K, the thermal energy (~0.6 kcal/mol) is much smaller than the bond strength (~200 kcal/mol from k=200 kcal/mol/Å²).

Part 3: Ethane - Adding Dihedrals¶

3.1 Why Dihedrals Matter¶

Ethane (H₃C-CH₃) can rotate around the C-C bond. The “staggered” conformation (H atoms far apart) is more stable than “eclipsed” (H atoms aligned).

With dihedrals, we can predict:

✅ Conformational preferences
✅ Rotational barriers
✅ Free energy landscapes
✅ Polymer chain statistics

3.2 LAMMPS Input File¶

The LAMMPS input file ethane_dihedral.in:

# Ethane (C2H6) - Full Molecular Force Field
# Tutorial 11 - Part 3

units           real
atom_style      full
boundary        p p p

read_data       ethane.data

# ========================================
# Force Field (OPLS-AA style)
# ========================================
pair_style      lj/cut/coul/cut 10.0
pair_coeff      1 1 0.066 3.50          # C (sp3)
pair_coeff      2 2 0.030 2.50          # H

bond_style      harmonic
bond_coeff      1 268.0 1.529           # C-C
bond_coeff      2 340.0 1.090           # C-H

angle_style     harmonic
angle_coeff     1 37.5 110.7            # C-C-H
angle_coeff     2 33.0 107.8            # H-C-H

dihedral_style  opls
dihedral_coeff  1 0.0 0.0 0.3 0.0       # H-C-C-H (V3 = 0.3 kcal/mol)

special_bonds   lj/coul 0.0 0.0 0.5

# ========================================
# Settings
# ========================================
neighbor        2.0 bin
neigh_modify    delay 0 every 1

# ========================================
# Equilibration
# ========================================
velocity        all create 300.0 44444
fix             1 all nvt temp 300.0 300.0 100.0

timestep        0.5                     # 0.5 fs (stiff bonds/angles)
thermo          1000
thermo_style    custom step temp pe ke etotal

run             50000                   # 25 ps equilibration

# ========================================
# Production - Track Dihedral Angle
# ========================================
reset_timestep  0

# Compute dihedral angles
compute         dihedrals all property/local dtype datom1 datom2 datom3 datom4 dphi
fix             2 all ave/time 1 1 10 c_dihedrals[*] file dihedrals.dat mode vector

dump            1 all custom 500 ethane.lammpstrj id mol type x y z
run             200000                  # 100 ps production

3.3 Building Ethane¶

The data file ethane.data is already provided in the shared directory:

ls /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part3_ethane/ethane.data

from ase.build import molecule
import numpy as np

print("="*60)
print("BUILDING ETHANE MOLECULE: ASE + LAMMPS File Format")
print("="*60)
print("\nStep 1: Use ASE for molecular geometry")
print("-" * 60)

# Single ethane molecule in a box
box_size = 15.0

# Get ethane from ASE database
ethane = molecule('C2H6')

# Center molecule in box
center = np.array([box_size/2, box_size/2, box_size/2])
ethane.positions += center - ethane.get_center_of_mass()

print(f"✓ ASE created ethane with {len(ethane)} atoms")

# Check geometry
c_atoms = [i for i, atom in enumerate(ethane) if atom.symbol == 'C']
c_c_dist = np.linalg.norm(ethane.positions[c_atoms[0]] - ethane.positions[c_atoms[1]])
print(f"✓ C-C bond: {c_c_dist:.3f} Å (correct from ASE)")
print(f"✓ Staggered conformation (minimum energy)")

print("\nStep 2: Write LAMMPS data file")
print("-" * 60)

# Extract data from ASE
positions = ethane.get_positions()
symbols = ethane.get_chemical_symbols()

# Write LAMMPS data file
with open('ethane.data', 'w') as f:
    # Header
    f.write('Ethane Molecule Built with ASE\n\n')
    f.write('8 atoms\n')
    f.write('7 bonds\n')
    f.write('12 angles\n')
    f.write('9 dihedrals\n\n')
    
    f.write('2 atom types  # 1=C, 2=H\n')
    f.write('2 bond types  # 1=C-C, 2=C-H\n')
    f.write('2 angle types # 1=C-C-H, 2=H-C-H\n')
    f.write('1 dihedral types\n\n')
    
    # Box
    f.write(f'0.0 {box_size} xlo xhi\n')
    f.write(f'0.0 {box_size} ylo yhi\n')
    f.write(f'0.0 {box_size} zlo zhi\n\n')
    
    # Masses
    f.write('Masses\n\n')
    f.write('1 12.011  # C\n')
    f.write('2 1.008   # H\n\n')
    
    # Atoms (full style)
    f.write('Atoms  # full\n\n')
    for i in range(len(positions)):
        atom_id = i + 1
        
        # Determine type and charge
        if symbols[i] == 'C':
            atom_type = 1
            charge = -0.18
        else:  # H
            atom_type = 2
            charge = 0.06
        
        pos = positions[i]
        f.write(f'{atom_id} 1 {atom_type} {charge} {pos[0]:.6f} {pos[1]:.6f} {pos[2]:.6f}\n')
    
    # Bonds - ASE orders atoms as: C, C, H, H, H, H, H, H
    f.write('\nBonds\n\n')
    f.write('1 1 1 2\n')        # C-C
    f.write('2 2 1 3\n')        # C-H
    f.write('3 2 1 4\n')        # C-H
    f.write('4 2 1 5\n')        # C-H
    f.write('5 2 2 6\n')        # C-H
    f.write('6 2 2 7\n')        # C-H
    f.write('7 2 2 8\n')        # C-H
    
    # Angles
    f.write('\nAngles\n\n')
    angles = [
        (1, 3, 1, 2), (1, 4, 1, 2), (1, 5, 1, 2),  # H-C-C
        (1, 6, 2, 1), (1, 7, 2, 1), (1, 8, 2, 1),  # H-C-C
        (2, 3, 1, 4), (2, 3, 1, 5), (2, 4, 1, 5),  # H-C-H
        (2, 6, 2, 7), (2, 6, 2, 8), (2, 7, 2, 8),  # H-C-H
    ]
    for aid, (atype, a1, a2, a3) in enumerate(angles, 1):
        f.write(f'{aid} {atype} {a1} {a2} {a3}\n')
    
    # Dihedrals (H-C-C-H)
    f.write('\nDihedrals\n\n')
    dihedrals = [
        (1, 3, 1, 2, 6), (1, 3, 1, 2, 7), (1, 3, 1, 2, 8),
        (1, 4, 1, 2, 6), (1, 4, 1, 2, 7), (1, 4, 1, 2, 8),
        (1, 5, 1, 2, 6), (1, 5, 1, 2, 7), (1, 5, 1, 2, 8),
    ]
    for did, (dtype, a1, a2, a3, a4) in enumerate(dihedrals, 1):
        f.write(f'{did} {dtype} {a1} {a2} {a3} {a4}\n')

print(f"✓ Wrote ethane.data with full topology:")
print(f"   - 7 bonds (1 C-C + 6 C-H)")
print(f"   - 12 angles (6 C-C-H + 6 H-C-H)")
print(f"   - 9 dihedrals (all H-C-C-H)")
print(f"✓ File ready for LAMMPS!\n")

print("="*60)
print("KEY POINTS:")
print("="*60)
print("• ASE provides correct ethane geometry (staggered)")
print("• All bonds, angles, dihedrals written explicitly")
print("• Dihedral potentials will show rotation barrier")
print("="*60)

3.4 Analysis: Dihedral Rotation¶

Using pre-run results:

ls /anvil/projects/x-chm250117/class_examples/tutorial11_completed/part3_ethane/ethane.lammpstrj

Notebook Cell 1: Dihedral Angle Analysis from Trajectory

import numpy as np
import matplotlib.pyplot as plt

# Set path to pre-run simulation
sim_dir = "/anvil/projects/x-chm250117/class_examples/tutorial11_completed/part3_ethane"

def calculate_dihedral(p1, p2, p3, p4):
    """Calculate dihedral angle between 4 points in degrees"""
    b1 = p2 - p1
    b2 = p3 - p2
    b3 = p4 - p3
    
    n1 = np.cross(b1, b2)
    n2 = np.cross(b2, b3)
    
    n1 /= np.linalg.norm(n1)
    n2 /= np.linalg.norm(n2)
    
    m1 = np.cross(n1, b2 / np.linalg.norm(b2))
    
    x = np.dot(n1, n2)
    y = np.dot(m1, n2)
    
    return np.degrees(np.arctan2(y, x))

# Read trajectory
print("Reading trajectory and calculating dihedrals...")
dihedrals = []

with open(f'{sim_dir}/ethane.lammpstrj', 'r') as f:
    lines = f.readlines()

i = 0
while i < len(lines):
    if 'ITEM: TIMESTEP' in lines[i]:
        # Skip to atoms section
        while i < len(lines) and 'ITEM: ATOMS' not in lines[i]:
            i += 1
        i += 1
        
        # Read 8 atoms (ethane molecule)
        positions = []
        for j in range(8):
            parts = lines[i].split()
            x, y, z = float(parts[3]), float(parts[4]), float(parts[5])
            positions.append(np.array([x, y, z]))
            i += 1
        
        # Calculate one H-C-C-H dihedral (atoms 2, 0, 1, 5)
        # Atom ordering from ASE: C(0), C(1), H(2), H(3), H(4), H(5), H(6), H(7)
        phi = calculate_dihedral(positions[2], positions[0], positions[1], positions[5])
        dihedrals.append(phi)
    else:
        i += 1

dihedrals = np.array(dihedrals)

plt.figure(figsize=(12, 5))

# Plot 1: Time series
plt.subplot(1, 2, 1)
plt.plot(dihedrals, 'b-', linewidth=0.5, alpha=0.7)
plt.axhline(y=60, color='g', linestyle='--', label='Staggered (60°)')
plt.axhline(y=-60, color='g', linestyle='--')
plt.axhline(y=180, color='r', linestyle='--', label='Staggered (180°)')
plt.axhline(y=0, color='orange', linestyle='--', label='Eclipsed (0°)')
plt.xlabel('Frame', fontsize=12)
plt.ylabel('Dihedral Angle (degrees)', fontsize=12)
plt.title('H-C-C-H Dihedral Angle vs Time', fontsize=13, fontweight='bold')
plt.ylim(-180, 180)
plt.legend(fontsize=9)
plt.grid(alpha=0.3)

# Plot 2: Histogram
plt.subplot(1, 2, 2)
plt.hist(dihedrals, bins=72, density=True, alpha=0.7, edgecolor='black', color='skyblue')
plt.axvline(x=60, color='g', linestyle='--', linewidth=2, label='Staggered')
plt.axvline(x=-60, color='g', linestyle='--', linewidth=2)
plt.axvline(x=180, color='g', linestyle='--', linewidth=2)
plt.axvline(x=0, color='orange', linestyle='--', linewidth=2, label='Eclipsed')
plt.xlabel('Dihedral Angle (degrees)', fontsize=12)
plt.ylabel('Probability Density', fontsize=12)
plt.title('Dihedral Angle Distribution', fontsize=13, fontweight='bold')
plt.xlim(-180, 180)
plt.legend(fontsize=10)
plt.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('ethane_dihedral.png', dpi=150)
plt.show()

# Calculate populations
staggered_count = np.sum((dihedrals > 45) & (dihedrals < 75)) + \
                  np.sum((dihedrals < -45) & (dihedrals > -75)) + \
                  np.sum(np.abs(dihedrals) > 165)
eclipsed_count = np.sum(np.abs(dihedrals) < 15)

total = len(dihedrals)
print(f"\nConformational Analysis:")
print(f"  Staggered conformations: {100*staggered_count/total:.1f}%")
print(f"  Eclipsed conformations: {100*eclipsed_count/total:.1f}%")
print(f"\n→ Staggered is strongly preferred (barrier ~3 kcal/mol)")
print(f"→ At 300K, eclipsed is rarely sampled")

Expected Result:

Strong peaks at ±60° and 180° (staggered)
Almost no population at 0°, ±120° (eclipsed)
Occasional barrier crossings (jumps between staggered states)

Summary and Key Lessons¶

What We Learned¶

Potential	System	New Physics	Predicted Properties
LJ only	Argon liquid	van der Waals forces	Structure (RDF), density, diffusion
LJ only	N₂ (fails!)	None — molecules dissociate	❌ Cannot hold covalent bonds
+ Bonds	N₂ gas	Permanent bonds	✅ Molecular identity, vibrations
+ Dihedrals	Ethane	Internal rotation	Conformational equilibria

Critical Observations¶

LJ cannot create covalent bonds — Two N atoms at 1.1 Å are INSIDE the LJ repulsive wall (σ = 3.31 Å). They fly apart!
Bonds create molecules — Harmonic bond at r₀ = 1.1 Å keeps N₂ together. But molecules are still floppy without angles.
Angles create shapes — Water’s 104.5° angle is why ice floats
Dihedrals create conformational complexity — Essential for proteins, polymers

Computational Cost Scaling¶

From the log files, compare:

Argon (LJ): ~X timesteps/second
N₂ (bonds): ~0.8X timesteps/second (slightly slower)
Water (angles + Coulomb): ~0.3X timesteps/second (PME dominates)
Ethane (full FF): ~0.6X timesteps/second

Lesson: Electrostatics (long-range Coulomb) is the main bottleneck, not bonded terms.

Exercise Questions¶

Argon: What happens to the RDF first peak if you increase temperature to 500 K? Run and compare.
N₂: Calculate the vibrational frequency from bond length fluctuations. Does it match experimental IR spectroscopy?
Water: What is the coordination number (average neighbors within first shell)? Should be ~4-5 for tetrahedral H-bonding.
Ethane: Heat to 600 K. Does the barrier crossing rate increase? Why?

Going Further¶

In Tutorial 12, we will explore:

Complete force field usage (AMBER/CHARMM parameter files)
EAM for metal nanoparticles
ReaxFF for methane combustion at high temperature

These potentials enable even more complex physics—but the fundamental lesson remains: use the simplest potential that captures your phenomenon.