Building Systems#

The build module purpose is to automate the construction of complex atomistic systems, from simple crystal bulks to intricate solid-liquid interfaces.

Workflow Overview #

Building a simulation-ready system generally follows these three steps:

Blueprint: Define composition and target density with a builder (SolutionBuilder, GlassBuilder, …). A blueprint describes what to build, not how.
Construction: Call the blueprint’s build() method (or an AtomicSystem method such as add_liquid_layer()) to generate the system, usually via Packmol.
Refinement: Assign force-field parameters, apply topology rules, or repair broken bonds as needed.

Note

Every builder that calls Packmol requires the packmol binary to be installed and available on your $PATH (see Installation).

The water template’s own O/H charges come along with it (here from its built-in .lt file). Na/Cl come from a generated single-atom template with no assigned charge (0.0) – assign real force-field charges afterward with set_ff_from_database() or set_charges() (see the force-field guide).

Pure water, or a hemispherical droplet instead of a box, are also available:

# Pure water
water = SolutionBuilder.from_water(density=1.0)
box_of_water = water.build(box=[20, 20, 20])

# A standalone hemispherical droplet
droplet = water.build_hemisphere(radius=15.0, axis="z")

Building a Glass Structure #

GlassBuilder builds an amorphous structure from a composition given as elements, oxides, or a mix of both:

from cemd.build import GlassBuilder

# By oxides
blueprint = GlassBuilder(
    density=2.3,
    composition={"SiO2": 3, "Al2O3": 2, "Na2O": 2},
)
glass = blueprint.build(box=[25, 25, 25])

print(glass.density)  # close to 2.3 g/cm³

# Equivalently, by elements
blueprint = GlassBuilder.from_elements(
    density=2.3,
    elements={"Si": 3, "Al": 4, "Na": 4, "O": 14},
)

If the composition is not charge-neutral, GlassBuilder raises a UserWarning when the blueprint is created (it still builds the structure).

Building C-S-H Models #

For cementitious materials, CSHBuilder generates C-S-H structures from a target Calcium/Silicate (Ca/Si) and Water/Silicate (H2O/Si) ratio, starting from a tobermorite crystal:

from cemd.build import CSHBuilder

builder = CSHBuilder(cs_ratio=1.5, ws_ratio=1.2)
csh_system = builder.build(
    supercell=[3, 5, 1],           # replication of the base cell
    model="tob11a_hamid.cif",      # or "tob11a_merlino.cif"
    min_mcl=3.0,                   # minimum mean chain length
)

# The builder keeps an analysis of the built structure
analysis = builder.analyze()
print(f"Real Ca/Si: {analysis['Ca/Si']:.3f}")
print(f"MCL: {analysis['MCL']:.2f}")
print(f"Qⁿ distribution: {analysis['Qn_distribution']}")

CSHBuilder can also wrap an existing C-S-H AtomicSystem (e.g. one loaded from a file) to analyze or convert it without rebuilding it:

system = AtomicSystem.from_file("csh_structure.data")
builder = CSHBuilder.from_system(system)
analysis = builder.analyze()

C-A-S-H Models (Aluminum Substitution)#

Convert a C-S-H structure to C-A-S-H by substituting a fraction of the bridging silicon atoms with aluminum:

cash_system = builder.to_cash(as_ratio=0.10)  # target Al/Si ratio

Note

Substitution targets bridging silicon atoms specifically (matching real Al-for-Si substitution chemistry), so it requires enough of them to be present – a very short or non-periodic silicate chain may have none, in which case the system is returned unchanged.

Building AFt/AFm Structures #

AFBuilder builds AFt (ettringite) or AFm structures by hydrating a sulfate-containing framework:

from cemd.build import AFBuilder

# AFt (ettringite)
builder = AFBuilder(ws_ratio=10.0)
aft_system = builder.build_aft(supercell=[2, 2, 1])

# AFm
builder = AFBuilder(ws_ratio=8.0)
afm_system = builder.build_afm(supercell=[3, 1, 1])

Creating Surfaces #

SurfaceBuilder cuts crystalline surfaces from a bulk structure along a chosen crystallographic plane and adds a vacuum layer.

from cemd import AtomicSystem
from cemd.build import SurfaceBuilder

# Load a bulk structure
bulk = AtomicSystem.from_file("tobermorite.cif")

# Generate surfaces cut along the (1, 0, 4) plane
builder = SurfaceBuilder(bulk)
surfaces, shifts, dipoles, broken_bonds = builder.build(
    miller_indices=[1, 0, 4],
    min_slab_size=25.0,   # Å
    min_vacuum_size=15.0,  # Å
)

# Select the first generated surface
my_surface = surfaces[0]

print(f"Generated {len(surfaces)} surfaces")
print(f"Broken bonds: {broken_bonds}")

Understanding the Parameters #

Miller Indices: Define the cutting plane orientation ([1,0,0], [1,1,0], [1,1,1], etc.)
min_slab_size: Minimum thickness of the material (too thin = bulk properties not representative)
min_vacuum_size: “Gap” in PBC to prevent surface-surface interactions
Broken Bonds: Number of bonds broken during the cut (diagnostic metric); build() automatically tries increasing tolerances up to max_broken_bonds until it finds a valid cut

Creating Solid/Liquid Interfaces #

Once you have a solid surface (and optionally a SolutionBuilder blueprint), add liquid directly onto it with AtomicSystem methods:

from cemd.build import SolutionBuilder

blueprint = SolutionBuilder(density=1.0, counts={"Na": 10, "Cl": 10})

# Add a flat liquid layer on top of the surface (in-place, returns self)
my_surface.add_liquid_layer(
    blueprint,
    thickness=20.0,  # Å
    distance=2.0,    # gap between solid and liquid
    vacuum=10.0,     # vacuum above the liquid
)

# Or a hemispherical droplet instead of a flat layer
droplet_blueprint = SolutionBuilder(density=1.0, counts={"Ca": 5, "Cl": 10})
my_surface.add_droplet(
    droplet_blueprint,
    radius=15.0,
    distance=2.0,
)

Adding Custom Structures #

To place any existing AtomicSystem on top of a surface (aligned by center of mass, at a fixed distance), use add_structure():

from cemd import AtomicSystem

custom = AtomicSystem.from_file("custom_molecule.pdb")

my_surface.add_structure(
    custom,
    distance=3.0,
    axis="z",
    vacuum=15.0,
)

Splitting Systems #

Splitter cuts a system open along an axis, creating a gap, and can optionally fill that gap with a solution in the same call (fluent API):

from cemd.build import Splitter, SolutionBuilder

# Simple split: cut at z=15 and open a 20 A gap there
split_system = Splitter(system, coordinate=15.0, axis="z", gap_size=20.0).split()

# Split and fill the gap with a solution
blueprint = SolutionBuilder.from_water(density=1.0)
filled_system = (
    Splitter(system, coordinate=15.0, axis="z", gap_size=30.0)
    .add_solution(blueprint, padding=2.0)
    .split()
)

Any bond stretched across the newly opened gap (and any angle/dihedral/ improper built on it) is automatically removed.

Choosing where to cut #

A cut can also sever contacts that carry no explicit bond – the Si-O framework under ClayFF/CSHFF is treated as non-bonded, so nothing in the topology marks it, yet slicing through it leaves broken tetrahedra. count_broken_bonds() and scan_broken_bonds() find those contacts geometrically, from a cutoff table, so you can pick a cut plane that goes through an interlayer rather than through a silicate sheet:

splitter = Splitter(csh_system, coordinate=0.0, axis="z", gap_size=20.0)

scan = splitter.scan_broken_bonds(step=0.5)
print(scan.head())

   coordinate  n_broken
  0.003344         0
  0.503344        15
  1.003344        19
  1.503344        24
  2.003344        37

# Cut where the structure suffers least
best = scan.loc[scan["n_broken"].idxmin(), "coordinate"]
result = splitter.with_coordinate(best).split()

The cutoff table (bonds_dict) accepts atom types first, elements as a fallback – the same layered convention as SurfaceBuilder. Types matter because elements alone cannot separate structure-bearing oxygens from ones that merely sit nearby:

# Default table keys on elements, so water O-H counts as breakable too
splitter.count_broken_bonds()

# Framework only: water is free to sit anywhere across the cut
splitter.count_broken_bonds(bonds_dict={("Si", "O"): 1.8})

Repairing what the cut breaks #

When a cut through the framework is unavoidable, repair=True caps the atoms it leaves under-coordinated: an exposed oxygen is protonated into a hydroxyl, and an exposed cation gets a hydroxyl group back in the direction its partner used to occupy.

splitter = Splitter(csh_system, coordinate=7.5, axis="z", gap_size=20.0)
result = splitter.split(repair=True)

print(splitter.repair_report)

{'broken': 66, 'capped': 192, 'skipped': 3}

skipped counts caps that would have landed on top of an existing atom — a bridging oxygen that loses both its cations would otherwise have each of them restore a copy of it at the same site; those atoms are left under-coordinated rather than duplicated.

Note

The added atoms are placed collinearly with the broken contact and carry a zero charge: this is a starting geometry to be relaxed. Re-run set_topology() and set_ff_from_database() afterwards so the new atoms get proper types and charges.

Building Systems#

Workflow Overview #

Construction Examples #

Building a Liquid Solution #

Building a Glass Structure #

Building C-S-H Models #

C-A-S-H Models (Aluminum Substitution)#

Building AFt/AFm Structures #

Creating Surfaces #

Understanding the Parameters #

Creating Solid/Liquid Interfaces #

Adding Custom Structures #

Splitting Systems #

Choosing where to cut #

Repairing what the cut breaks #