Skip to content

Periodic Systems: Crystals and Surfaces

What You'll Learn

  • How to set up periodic boundary conditions (PBC) with AIMNet2 and ASE
  • How the calculator handles long-range electrostatics in periodic systems
  • When to use DSF, Ewald, or PME for Coulomb interactions
  • How to optimize unit cell parameters using the stress tensor

Prerequisites

  • AIMNet2 installed with ASE support (pip install "aimnet[ase]")
  • Familiarity with ASE Atoms objects (see Geometry Optimization)
  • Basic understanding of periodic boundary conditions in atomistic simulation

Step 1: Setting Up a Periodic System

Periodic systems require a unit cell and PBC flags. With ASE, you define these on the Atoms object:

from ase import Atoms
from aimnet.calculators import AIMNet2ASE

# Formaldehyde molecular crystal (P2_1/c, simplified)
# 4 molecules in the unit cell
atoms = Atoms(
    symbols="C4O4H8",
    positions=[
        # Molecule 1
        [0.00, 0.00, 0.00],   # C
        [1.20, 0.00, 0.00],   # O
        [-0.55, 0.94, 0.00],  # H
        [-0.55, -0.94, 0.00], # H
        # Molecule 2
        [3.80, 2.50, 1.80],   # C
        [5.00, 2.50, 1.80],   # O
        [3.25, 3.44, 1.80],   # H
        [3.25, 1.56, 1.80],   # H
        # Molecule 3
        [1.90, 5.00, 3.60],   # C
        [3.10, 5.00, 3.60],   # O
        [1.35, 5.94, 3.60],   # H
        [1.35, 4.06, 3.60],   # H
        # Molecule 4
        [5.70, 2.50, 5.40],   # C
        [6.90, 2.50, 5.40],   # O
        [5.15, 3.44, 5.40],   # H
        [5.15, 1.56, 5.40],   # H
    ],
    cell=[7.60, 5.00, 7.20],
    pbc=True,
)

calc = AIMNet2ASE("aimnet2")
atoms.calc = calc

Note

When pbc=True, the calculator always uses sparse mode with neighbor lists, regardless of system size or nb_threshold. This ensures periodic images are handled correctly through cell shift vectors.

Step 2: Understanding the Coulomb Auto-Switch

When you run a periodic calculation for the first time, the calculator detects that your system has a unit cell and automatically switches the Coulomb method:

energy = atoms.get_potential_energy()
# UserWarning: Switching to DSF Coulomb for PBC

Why does this happen? The default Coulomb method is "simple" (all-pairs 1/r sum), which does not account for periodic images. For periodic systems, this would give incorrect electrostatics. The calculator detects PBC and automatically switches to the Damped Shifted Force (DSF) method, which properly handles the infinite lattice sum with a smooth cutoff.

Warning

The auto-switch warning appears once per calculator when PBC is first detected. To avoid the warning entirely, set the Coulomb method explicitly before your first periodic calculation (see Step 3).

Step 3: Choosing a Coulomb Method for PBC

For periodic systems, you have three options for long-range electrostatics. Set the method on the underlying AIMNet2Calculator via the base_calc attribute:

DSF (Damped Shifted Force) is the recommended method for routine periodic calculations. It uses a neighbor-list-based cutoff with smooth damping, giving O(N) scaling:

calc.base_calc.set_lrcoulomb_method("dsf", cutoff=15.0, dsf_alpha=0.2)

energy = atoms.get_potential_energy()
forces = atoms.get_forces()

The default parameters (cutoff=15.0, dsf_alpha=0.2) work well for most molecular crystals. For dense systems, you can reduce the cutoff to 12 Angstrom. For surfaces or dilute systems, consider increasing to 18-20 Angstrom.

Ewald (High-Accuracy Benchmarks)

Ewald summation splits the Coulomb sum into real-space and reciprocal-space components. The splitting parameter and real/reciprocal cutoffs are estimated per call from ewald_accuracy (default 1e-6, matching the nvalchemiops default):

calc.base_calc.set_lrcoulomb_method("ewald")

energy = atoms.get_potential_energy()
forces = atoms.get_forces()
stress = atoms.get_stress()

DSF supports inference forces and stress, but force/stress losses (train=True with forces=True or stress=True) and Hessians are not supported. Ewald and PME support force/stress losses in training, but not calculator Hessian requests. See Long-Range Methods → Derivative Support.

When to use Ewald:

  • Validating DSF results for a new system type
  • Computing precise lattice energies for benchmarking
  • Systems where electrostatics dominate (ionic molecular crystals)

PME (Large Cells)

Particle Mesh Ewald uses the same accuracy machinery as Ewald but evaluates the reciprocal sum on a B-spline mesh via FFT, which scales better on large cells:

calc.base_calc.set_lrcoulomb_method("pme")

energy = atoms.get_potential_energy()
stress = atoms.get_stress()

PME and Ewald share the same ewald_accuracy, neighbor-list construction, and calculator derivative behavior; PME tends to be cheaper as the cell grows.

Comparing Methods

You can verify convergence by comparing DSF, Ewald, and PME on the same system:

calc.base_calc.set_lrcoulomb_method("dsf", cutoff=15.0)
energy_dsf = atoms.get_potential_energy()

calc.base_calc.set_lrcoulomb_method("ewald", ewald_accuracy=1e-6)
energy_ewald = atoms.get_potential_energy()

calc.base_calc.set_lrcoulomb_method("pme", ewald_accuracy=1e-6)
energy_pme = atoms.get_potential_energy()

print(f"DSF energy:   {energy_dsf:.6f} eV")
print(f"Ewald energy: {energy_ewald:.6f} eV")
print(f"PME energy:   {energy_pme:.6f} eV")

Step 4: Computing the Stress Tensor

The stress tensor is essential for optimizing unit cell parameters. Request it via stress=True in the ASE calculator:

from ase.constraints import ExpCellFilter
from ase.optimize import BFGS

# Set up the calculator with DSF for PBC
calc = AIMNet2ASE("aimnet2")
calc.base_calc.set_lrcoulomb_method("dsf", cutoff=15.0)
atoms.calc = calc

# Verify stress is available
stress = atoms.get_stress()  # Voigt notation (6,) in ASE
print(f"Stress (Voigt): {stress}")

Tip

ASE returns stress in Voigt notation as a 6-element array (xx, yy, zz, yz, xz, xy). The AIMNet2 calculator internally computes the full 3x3 stress tensor, and ASE handles the conversion.

Step 5: Cell Optimization

To optimize both atomic positions and cell parameters simultaneously, wrap the Atoms object in ExpCellFilter (or FrechetCellFilter):

from ase.constraints import ExpCellFilter
from ase.optimize import BFGS

# ExpCellFilter allows simultaneous optimization of positions and cell
ecf = ExpCellFilter(atoms)

opt = BFGS(ecf, logfile="cell_opt.log")
opt.run(fmax=0.05)  # eV/Angstrom force convergence

print(f"Optimized cell: {atoms.cell.lengths()}")
print(f"Optimized angles: {atoms.cell.angles()}")
print(f"Final energy: {atoms.get_potential_energy():.6f} eV")

Why ExpCellFilter? It maps cell degrees of freedom onto an expanded coordinate space so that a standard optimizer (BFGS, FIRE) can optimize both positions and cell shape in a single run. The fmax criterion applies to both atomic forces and stress-derived cell forces.

Step 6: Worked Example -- Molecular Crystal Optimization

Here is a complete workflow for optimizing a molecular crystal structure:

from ase import Atoms
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS
from aimnet.calculators import AIMNet2ASE

# --- 1. Build the crystal ---
atoms = Atoms(
    symbols="C4O4H8",
    positions=[
        [0.00, 0.00, 0.00], [1.20, 0.00, 0.00],
        [-0.55, 0.94, 0.00], [-0.55, -0.94, 0.00],
        [3.80, 2.50, 1.80], [5.00, 2.50, 1.80],
        [3.25, 3.44, 1.80], [3.25, 1.56, 1.80],
        [1.90, 5.00, 3.60], [3.10, 5.00, 3.60],
        [1.35, 5.94, 3.60], [1.35, 4.06, 3.60],
        [5.70, 2.50, 5.40], [6.90, 2.50, 5.40],
        [5.15, 3.44, 5.40], [5.15, 1.56, 5.40],
    ],
    cell=[7.60, 5.00, 7.20],
    pbc=True,
)

# --- 2. Attach calculator with DSF Coulomb ---
calc = AIMNet2ASE("aimnet2")
calc.base_calc.set_lrcoulomb_method("dsf", cutoff=15.0)
atoms.calc = calc

# --- 3. Optimize cell + positions ---
ecf = ExpCellFilter(atoms)
opt = BFGS(ecf, logfile="crystal_opt.log")
opt.run(fmax=0.05)

# --- 4. Report results ---
print(f"Converged in {opt.get_number_of_steps()} steps")
print(f"Final energy: {atoms.get_potential_energy():.4f} eV")
print(f"Cell lengths: {atoms.cell.lengths()}")
print(f"Cell angles:  {atoms.cell.angles()}")

Tip

For reading crystal structures from CIF files, use ase.io.read("structure.cif"). ASE will parse the cell parameters and symmetry automatically.

Summary of Coulomb Methods for PBC

Method Scaling PBC Support Accuracy control Best For
simple O(N^2) No Exact (non-PBC) Auto-switches to DSF for PBC
dsf O(N) Yes cutoff, dsf_alpha Routine PBC, MD, optimization
ewald ~O(N^{3/2}) Yes ewald_accuracy (default 1e-6) High-accuracy reference
pme ~O(N log N) Yes ewald_accuracy (default 1e-6) Large cells

What's Next