Optimize¶

The optimize module contains geometry optimization algorithms.

Base¶

Base classes for geometry optimization.

OptimizationResult `dataclass` ¶

OptimizationResult(structure, energy, fmax, nsteps, converged, reason)

Bases: Generic[S]

Result of a geometry optimization.

Attributes:

Name	Type	Description
`structure`	`S`	Final optimized structure (same object if inplace=True)
`energy`	`float`	Final potential energy
`fmax`	`float`	Final maximum force component
`nsteps`	`int`	Number of optimization steps taken
`converged`	`bool`	Whether convergence criteria were met
`reason`	`str`	Human-readable termination reason

Optimizer ¶

Optimizer(potential, *, entity_type=Entity)

Bases: ABC, Generic[S]

Base class for structure optimizers.

Works with any StructLike (Struct, Atomistic, CoarseGrain, etc.) that: - Has .entities (TypeBucket) containing entities with "xyz" field - Has .to_frame() method to convert to Frame format

The optimizer calls potential.calc_energy(frame) and potential.calc_forces(frame) directly - each potential is responsible for extracting what it needs from Frame.

Parameters:

Name	Type	Description	Default
`potential`	`PotentialLike`	Potential with calc_energy/calc_forces methods	required
`entity_type`	`type[Entity]`	Type of entity to optimize (default: Entity for all)	`Entity`

Example

from molpy.optimize import LBFGS from molpy.potential.bond import Harmonic

potential = Harmonic(k=100.0, r0=1.5) opt = LBFGS(potential, maxstep=0.04, memory=20) result = opt.run(struct, fmax=0.01, steps=500)

Source code in src/molpy/optimize/base.py

def __init__(
    self,
    potential: PotentialLike,
    *,
    entity_type: type[Entity] = Entity,
) -> None:
    self.potential = potential
    self.entity_type = entity_type
    self._callbacks: list[tuple[Callable, int, dict]] = []

attach ¶

attach(func, interval=1, **kwargs)

Attach a callback function.

The callback will be called every interval steps with the optimizer instance and current structure as arguments.

Parameters:

Name	Type	Description	Default
`func`	`Callable`	Callback function(optimizer, structure, **kwargs)	required
`interval`	`int`	Call every N steps	`1`
`**kwargs`	`Any`	Additional arguments for callback	`{}`

Source code in src/molpy/optimize/base.py

def attach(self, func: Callable, interval: int = 1, **kwargs: Any) -> None:
    """Attach a callback function.

    The callback will be called every `interval` steps with the optimizer
    instance and current structure as arguments.

    Args:
        func: Callback function(optimizer, structure, **kwargs)
        interval: Call every N steps
        **kwargs: Additional arguments for callback
    """
    self._callbacks.append((func, interval, kwargs))

get_energy ¶

get_energy(structure)

Compute energy for structure.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to evaluate	required

Returns:

Type	Description
`float`	Potential energy

Source code in src/molpy/optimize/base.py

def get_energy(self, structure: S) -> float:
    """Compute energy for structure.

    Args:
        structure: Structure to evaluate

    Returns:
        Potential energy
    """
    energy, _ = self.get_energy_and_forces(structure)
    return energy

get_energy_and_forces ¶

get_energy_and_forces(structure)

Compute energy and forces via Frame interface.

Converts structure to Frame and calls potential methods directly. Each potential is responsible for extracting what it needs from Frame.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to evaluate	required

Returns:

Type	Description
`tuple[float, ndarray]`	(energy, forces) where forces is (N, 3) array

Source code in src/molpy/optimize/base.py

def get_energy_and_forces(self, structure: S) -> tuple[float, np.ndarray]:
    """Compute energy and forces via Frame interface.

    Converts structure to Frame and calls potential methods directly.
    Each potential is responsible for extracting what it needs from Frame.

    Args:
        structure: Structure to evaluate

    Returns:
        (energy, forces) where forces is (N, 3) array
    """
    # Convert structure to Frame
    frame = structure.to_frame()

    # Call potential methods directly - they handle Frame extraction
    energy = self.potential.calc_energy(frame)
    forces = self.potential.calc_forces(frame)

    return energy, forces

get_forces ¶

get_forces(structure)

Compute forces for structure as (N, 3) array.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to evaluate	required

Returns:

Type	Description
`ndarray`	(N, 3) array of forces

Source code in src/molpy/optimize/base.py

def get_forces(self, structure: S) -> np.ndarray:
    """Compute forces for structure as (N, 3) array.

    Args:
        structure: Structure to evaluate

    Returns:
        (N, 3) array of forces
    """
    _, forces = self.get_energy_and_forces(structure)
    return forces

get_positions ¶

get_positions(structure)

Extract positions as (N, 3) array from entities.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to extract positions from	required

Returns:

Type	Description
`ndarray`	(N, 3) numpy array of positions

Source code in src/molpy/optimize/base.py

def get_positions(self, structure: S) -> np.ndarray:
    """Extract positions as (N, 3) array from entities.

    Args:
        structure: Structure to extract positions from

    Returns:
        (N, 3) numpy array of positions
    """
    entities = structure.entities.all()
    if not entities:
        return np.empty((0, 3), dtype=float)

    # Use x, y, z fields (never use xyz)
    x_list = entities["x"]
    y_list = entities["y"]
    z_list = entities["z"]
    positions = np.column_stack([x_list, y_list, z_list])

    if positions.ndim == 1:
        positions = positions.reshape(-1, 3)
    return positions

run ¶

run(structure, fmax=0.01, steps=1000, *, inplace=True)

Run optimization until convergence or max steps.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to optimize	required
`fmax`	`float`	Convergence threshold (max force component)	`0.01`
`steps`	`int`	Maximum number of steps	`1000`
`inplace`	`bool`	If True, modify structure in-place; if False, work on copy	`True`

Returns:

Type	Description
`OptimizationResult[S]`	OptimizationResult with final state

Source code in src/molpy/optimize/base.py

def run(
    self,
    structure: S,
    fmax: float = 0.01,
    steps: int = 1000,
    *,
    inplace: bool = True,
) -> OptimizationResult[S]:
    """Run optimization until convergence or max steps.

    Args:
        structure: Structure to optimize
        fmax: Convergence threshold (max force component)
        steps: Maximum number of steps
        inplace: If True, modify structure in-place; if False, work on copy

    Returns:
        OptimizationResult with final state
    """
    # Make copy if needed
    if inplace:
        working_structure = structure
    else:
        from copy import deepcopy

        working_structure = deepcopy(structure)

    energy = 0.0
    current_fmax = float("inf")
    nsteps = 0

    for i in range(steps):
        energy, current_fmax = self.step(working_structure)
        nsteps += 1

        # Call callbacks
        for callback, interval, kwargs in self._callbacks:
            if nsteps % interval == 0:
                callback(self, working_structure, **kwargs)

        # Check convergence
        if current_fmax < fmax:
            return OptimizationResult(
                structure=working_structure,
                energy=energy,
                fmax=current_fmax,
                nsteps=nsteps,
                converged=True,
                reason=f"Converged: fmax={current_fmax:.6f} < {fmax}",
            )

    # Max steps reached
    return OptimizationResult(
        structure=working_structure,
        energy=energy,
        fmax=current_fmax,
        nsteps=nsteps,
        converged=False,
        reason=f"Max steps reached: {steps}",
    )

set_positions ¶

set_positions(structure, positions)

Write positions back into structure entities.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to update	required
`positions`	`ndarray`	(N, 3) array of new positions	required

Source code in src/molpy/optimize/base.py

def set_positions(self, structure: S, positions: np.ndarray) -> None:
    """Write positions back into structure entities.

    Args:
        structure: Structure to update
        positions: (N, 3) array of new positions
    """
    entities = structure.entities.all()
    positions = positions.reshape(-1, 3)
    for i, entity in enumerate(entities):
        pos = positions[i]
        # Update x, y, z fields (never use xyz)
        entity["x"] = float(pos[0])
        entity["y"] = float(pos[1])
        entity["z"] = float(pos[2])

step `abstractmethod` ¶

step(structure)

Perform one optimization step.

Parameters:

Name	Type	Description	Default
`structure`	`S`	Structure to optimize (modified in-place)	required

Returns:

Type	Description
`tuple[float, float]`	(energy, fmax) tuple where: energy: potential energy after step fmax: maximum force component after step

Source code in src/molpy/optimize/base.py

@abstractmethod
def step(self, structure: S) -> tuple[float, float]:
    """Perform one optimization step.

    Args:
        structure: Structure to optimize (modified in-place)

    Returns:
        (energy, fmax) tuple where:
            energy: potential energy after step
            fmax: maximum force component after step
    """
    pass

PotentialLike ¶

Bases: Protocol

Protocol for potential functions compatible with calc_energy_from_frame.

LBFGS¶

L-BFGS geometry optimizer.

LBFGS ¶

LBFGS(potential, *, maxstep=0.04, memory=20, damping=1.0, entity_type=None)

Bases: Optimizer[S]

Limited-memory BFGS geometry optimizer.

Implements the L-BFGS algorithm for quasi-Newton optimization with limited memory storage. Uses two-loop recursion to compute search directions efficiently.

Parameters:

Name	Type	Description	Default
`potential`		Potential with calc_energy/calc_forces methods	required
`maxstep`	`float`	Maximum step size (as displacement norm)	`0.04`
`memory`	`int`	Number of previous steps to store for Hessian approximation	`20`
`damping`	`float`	Damping factor for step size	`1.0`
`entity_type`		Type of entity to optimize	`None`

Attributes:

Name	Type	Description
`maxstep`		Maximum allowed step size
`memory`		LBFGS memory size
`damping`		Step damping factor
`s_history`	`list[ndarray]`	Position difference history
`y_history`	`list[ndarray]`	Gradient difference history
`rho_history`	`list[float]`	Curvature history (1 / y·s)

Example

from molpy.core.atomistic import Atomistic from molpy.potential.bond import Harmonic from molpy.optimize import LBFGS

struct = Atomistic()

... add atoms and bonds ...¶

potential = Harmonic(k=100.0, r0=1.5) opt = LBFGS(potential, maxstep=0.04, memory=20) result = opt.run(struct, fmax=0.01, steps=500)

Source code in src/molpy/optimize/lbfgs.py

def __init__(
    self,
    potential,
    *,
    maxstep: float = 0.04,
    memory: int = 20,
    damping: float = 1.0,
    entity_type=None,
) -> None:
    # Handle entity_type default
    if entity_type is None:
        from molpy.core.entity import Entity

        entity_type = Entity

    super().__init__(potential, entity_type=entity_type)
    self.maxstep = maxstep
    self.memory = memory
    self.damping = damping

    # LBFGS state (reset for each new structure/run)
    self._current_structure_id = None
    self.s_history: list[np.ndarray] = []  # position differences
    self.y_history: list[np.ndarray] = []  # gradient differences
    self.rho_history: list[float] = []  # 1 / (y · s)

    self._prev_positions: np.ndarray | None = None
    self._prev_gradient: np.ndarray | None = None

step ¶

step(structure)

Perform one L-BFGS optimization step.

Args: structure: Structure to optimize (modified in-place)

Returns:

Type	Description
`tuple[float, float]`	(energy, fmax) tuple where: energy: potential energy after step fmax: maximum force component after step

Source code in src/molpy/optimize/lbfgs.py

def step(self, structure: S) -> tuple[float, float]:
    """Perform one L-BFGS optimization step.

     Args:
         structure: Structure to optimize (modified in-place)

    Returns:
         (energy, fmax) tuple where:
             energy: potential energy after step
             fmax: maximum force component after step
    """
    # Reset state if this is a new structure
    self._reset_state(id(structure))

    # Get current state
    positions = self.get_positions(structure)
    forces = self.get_forces(structure)
    gradient = -forces  # Forces are -∇E, we need ∇E for minimization
    energy = self.get_energy(structure)

    # Flatten for linear algebra
    x = positions.reshape(-1)
    g = gradient.reshape(-1)

    # Update LBFGS history if we have a previous step
    if self._prev_positions is not None:
        s = x - self._prev_positions
        y = g - self._prev_gradient

        sy = np.dot(s, y)
        if sy > 1e-10:  # Only update if curvature condition holds
            self.s_history.append(s)
            self.y_history.append(y)
            self.rho_history.append(1.0 / sy)

            # Keep only last 'memory' updates
            if len(self.s_history) > self.memory:
                self.s_history.pop(0)
                self.y_history.pop(0)
                self.rho_history.pop(0)

    # Compute search direction using two-loop recursion
    # LBFGS returns H^{-1} * g, which is the direction for minimization
    search_dir = self._lbfgs_direction(g)

    # Apply step-size control
    step_size = self._compute_step_size(search_dir)

    # Update positions: move in the direction of -∇E (minimization)
    new_x = x - step_size * search_dir
    new_positions = new_x.reshape(positions.shape)
    self.set_positions(structure, new_positions)

    # Store for next iteration (before updating positions)
    self._prev_positions = x.copy()
    self._prev_gradient = g.copy()

    # Recompute energy and forces at new positions
    new_energy = self.get_energy(structure)
    new_forces = self.get_forces(structure)
    fmax = float(np.max(np.abs(new_forces)))

    return new_energy, fmax

Potential Wrappers¶

Potential wrappers for optimizer compatibility.

Since potentials in molpy work with raw arrays (positions, indices, types), not Frame objects, we need simple wrappers to extract data from Frames.

Each wrapper implements calc_energy(frame) and calc_forces(frame) by extracting the appropriate data and calling the underlying potential.

AnglePotentialWrapper ¶

AnglePotentialWrapper(potential)

Wrapper for angle potentials to work with Frame interface.

Extracts (r, angle_idx, angle_types) from Frame and calls underlying potential.

Source code in src/molpy/optimize/potential_wrappers.py

def __init__(self, potential):
    self.potential = potential
    self.type = "angle"
    self.name = getattr(potential, "name", "angle")

calc_energy ¶

calc_energy(frame)

Extract angle data from Frame and compute energy.

Source code in src/molpy/optimize/potential_wrappers.py

def calc_energy(self, frame):
    """Extract angle data from Frame and compute energy."""
    # Get coordinates from x, y, z fields (never use xyz)
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    angles = frame["angles"]
    angle_idx = angles[["atomi", "atomj", "atomk"]]
    angle_types = angles["type"]
    return self.potential.calc_energy(r, angle_idx, angle_types)

calc_forces ¶

calc_forces(frame)

Extract angle data from Frame and compute forces.

Source code in src/molpy/optimize/potential_wrappers.py

def calc_forces(self, frame):
    """Extract angle data from Frame and compute forces."""
    # Get coordinates from x, y, z fields (never use xyz)
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    angles = frame["angles"]
    angle_idx = angles[["atomi", "atomj", "atomk"]]
    angle_types = angles["type"]
    return self.potential.calc_forces(r, angle_idx, angle_types)

BondPotentialWrapper ¶

BondPotentialWrapper(potential)

Wrapper for bond potentials to work with Frame interface.

Extracts (r, bond_idx, bond_types) from Frame and calls underlying potential.

Source code in src/molpy/optimize/potential_wrappers.py

def __init__(self, potential):
    self.potential = potential
    self.type = "bond"
    self.name = getattr(potential, "name", "bond")

calc_energy ¶

calc_energy(frame)

Extract bond data from Frame and compute energy.

Source code in src/molpy/optimize/potential_wrappers.py

def calc_energy(self, frame):
    """Extract bond data from Frame and compute energy."""
    # Get coordinates from x, y, z fields (never use xyz)
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    bonds = frame["bonds"]
    bond_idx = bonds[["atomi", "atomj"]]
    bond_types = bonds["type"]
    return self.potential.calc_energy(r, bond_idx, bond_types)

calc_forces ¶

calc_forces(frame)

Extract bond data from Frame and compute forces.

Source code in src/molpy/optimize/potential_wrappers.py

def calc_forces(self, frame):
    """Extract bond data from Frame and compute forces."""
    # Get coordinates from x, y, z fields (never use xyz)
    x = frame["atoms"]["x"]
    y = frame["atoms"]["y"]
    z = frame["atoms"]["z"]
    r = np.column_stack([x, y, z])
    bonds = frame["bonds"]
    bond_idx = bonds[["atomi", "atomj"]]
    bond_types = bonds["type"]
    return self.potential.calc_forces(r, bond_idx, bond_types)

Optimize¶

Base¶

OptimizationResult dataclass ¶

Optimizer ¶

attach ¶

get_energy ¶

get_energy_and_forces ¶

get_forces ¶

get_positions ¶

run ¶

set_positions ¶

step abstractmethod ¶

PotentialLike ¶

LBFGS¶

LBFGS ¶

... add atoms and bonds ...¶

step ¶

Potential Wrappers¶

AnglePotentialWrapper ¶

calc_energy ¶

calc_forces ¶

BondPotentialWrapper ¶

calc_energy ¶

calc_forces ¶

OptimizationResult `dataclass` ¶

step `abstractmethod` ¶