GrowPAL

Growth Pattern Algorithm with diversity-preserving selection for global optimization of atomic clusters.

GrowPAL implements a progressive-growth strategy for exploring potential energy surfaces (PESs) of atomic clusters. Its key innovation is incorporating unsupervised structural classification into the selection stage, preserving morphological diversity across growth steps rather than selecting exclusively by energy.

Table of Contents

• Background
• How It Works
• Installation
• Dependencies
• Module Overview
• Usage
• Examples
• Key Results
• Citation
• Authors
• License

Background

Global optimization algorithms based on progressive growth exploit structural correlations between neighboring cluster sizes and incorporate information obtained during earlier stages of the search. However, selection schemes based exclusively on energy may reduce structural diversity when multiple competing funnels are present.

The organization of potential energy surfaces into competing funnels — regions of configuration space that converge toward distinct structural motifs, separated by significant energy barriers — is a fundamental feature of cluster energy landscapes. For example, in the LJ75–LJ78 size range the global minimum corresponds to a decahedral Marks structure, while a competing icosahedral funnel contains many low-energy local minima of comparable energy. Algorithms that rank structures solely by energy tend to oversample the dominant funnel, restricting access to alternative regions of the PES.

GrowPAL integrates data-driven structural classification with the growth-based framework. Structural descriptors and hierarchical clustering are used to preserve representatives from distinct morphological families throughout the growth process, enabling simultaneous exploration of multiple regions of configuration space.

This package is a continuation of the original GrowPAL algorithm (López-Castro, Ortiz-Chi, Merino, J. Chem. Theory Comput. 2024, 20, 4939–4948), which established the growth strategy and demonstrated its efficiency for LJ clusters up to 74 atoms. The present version introduces the diversity-preserving selection scheme that extends the method and finds the correct global minima for LJ clusters up to 130 atoms.

How It Works

Each growth iteration comprises four stages:

1. Structure generation. From each seed cluster of size N, candidate structures of size N+1 are generated by systematic interstitial insertion. Internal triangles are identified, a new atom is placed at each inequivalent triangle centroid (symmetry-equivalent sites are filtered), and the cluster is radially expanded around the insertion point. Each seed generates as many candidates as there are inequivalent internal triangles.

2. Structural discrimination. Duplicate and near-duplicate structures are removed via a two-stage USR-based filter (deduplicate_by_usr in aegon.discrimination.cluster_usr): a Manhattan-distance prefilter followed by a k-d tree neighborhood search. Two structures are considered redundant when their L1 descriptor distance falls below a threshold (default 0.99).

3. MBTR descriptor computation. Each surviving structure is represented by a combined Many-Body Tensor Representation vector encoding two-body (k=2, pairwise distances) and three-body (k=3, angular correlations) information, computed over 100 grid points with Gaussian broadenings of 0.01 Å (distances) and 0.01 rad (angles). For monoatomic systems a fast Numba-compiled path is used; for multi-species systems the DScribe library is used instead.

4. Clustering-based selection. Descriptor vectors are partitioned by agglomerative Ward clustering with a k-NN connectivity constraint. The selection_count lowest-energy structures from each cluster are retained as seeds for the next iteration. This selection_count parameter controls the trade-off between computational cost and search thoroughness: larger values increase population diversity and improve recovery of challenging global minima at the expense of more local optimizations.

The algorithm is fully deterministic: repeated runs with identical inputs produce byte-for-byte identical results.

Installation

pip install growpal

This automatically installs all required packages: aegon (which brings in ase, numpy, scipy, py3Dmol, numba, joblib), molsympy (symmetry-inequivalent-site detection used by the growth engine itself), scikit-learn (Ward clustering), and dscribe (MBTR descriptors for multi-species systems).

Requires Python >= 3.10.

Google Colab

!pip install growpal

Dependencies

Module Overview

Usage

from growpal.libgrowpal import growpal_parallel
from growpal.libselect_clustering import select_by_clustering

# poscarlist: your own list of N-atom ase.Atoms seed structures (e.g. from
# aegon.io.xyz.readxyzs) -- illustrative API call, not a runnable script on
# its own. See "Examples" below for a complete, runnable growth loop.

# Grow a list of N-atom clusters to N+1
candidates = growpal_parallel(poscarlist, specie='Mo', dtol=1.2, n_cores=2)

# Select preserving morphological diversity
selected = select_by_clustering(candidates, selection_count=160, n_clusters=8)

The selection_count parameter (analogous to N·f in the manuscript) controls population diversity. Higher values improve recovery of difficult global minima at the cost of more local optimizations per iteration.

Examples

1. Load and visualize a Sutton-Chen cluster (Google Colab)

The SC database in aegon covers 10 transition metals (Al, Ag, Au, Cu, Ir, Ni, Pb, Pd, Pt, Rh) for cluster sizes up to 90 atoms, using the original Sutton-Chen parameterization.

viewmol_ASE (used below) renders with py3Dmol and requires an active Jupyter/Colab notebook — it raises ImportError if run from a plain script or terminal. The data-access lines above it (get_sc_cluster, SUTTON_CHEN_PARAMS) work anywhere; only drop the viewmol_ASE(mol) call outside a notebook.

from aegon.data.sc import get_sc_cluster
from aegon.optimization.potentials.sc import SUTTON_CHEN_PARAMS
from aegon.io.gcolab import viewmol_ASE

mol = get_sc_cluster(58, symbol='Pb')  # ASE Atoms object
print(mol.info['i'])                   # cluster ID
print(mol.info['e'])                   # energy in eV

params = SUTTON_CHEN_PARAMS['Pb']
epsilon = params['epsilon']
print(mol.info['e'] / epsilon)         # energy in epsilon units

viewmol_ASE(mol)

2. Load and visualize a Lennard-Jones cluster (Google Colab)

The LJ database in aegon covers cluster sizes from 3 to 150 atoms (ε = 1 eV, equilibrium distance 2^(1/6) σ = 3 Å). All stored minima were validated against the Wales reference database.

As in example 1, viewmol_ASE requires an active Jupyter/Colab notebook and raises ImportError outside one; get_lj_cluster itself works anywhere.

from aegon.data.lj import get_lj_cluster
from aegon.io.gcolab import viewmol_ASE

atoms = get_lj_cluster(30)
viewmol_ASE(atoms)
print("%s  Energy=%f" % (atoms.info['i'], atoms.info['e']))

3. Run a GrowPAL optimization (Lennard-Jones potential, also runs in Google Colab)

Prepare an initial LJ010.xyz file with seed structures for the starting size — the first two pre-optimized 10-atom clusters from growpal's own examples/LJ/LJ010.xyz, giving the growth loop two independent starting points instead of one:

10
 -28.42253189 C3v
Mo      0.000000000      0.000000000      0.920238717
Mo      0.856362866     -1.483263994     -1.570679322
Mo      0.856362866      1.483263994     -1.570679322
Mo     -1.712725733     -0.000000000     -1.570679322
Mo     -2.933911374      0.000000000      1.128123544
Mo      2.826583404      0.000000000      0.135809539
Mo      1.466955687      2.540841783      1.128123544
Mo      1.466955687     -2.540841783      1.128123544
Mo     -1.413291702      2.447893033      0.135809539
Mo     -1.413291702     -2.447893033      0.135809539
10
 -27.55586304 Cs
Mo      0.011945677     -1.210020432      1.477438460
Mo     -2.261012587      0.085349055      0.000000000
Mo     -1.174656083      1.296071202      2.526807935
Mo      1.824519667      0.896630777      2.522222101
Mo     -1.963985406     -2.878475716      0.000000000
Mo     -1.174656083      1.296071202     -2.526807935
Mo      0.011945677     -1.210020432     -1.477438460
Mo      2.532352869     -0.597306362      0.000000000
Mo      0.369026603      1.425069930      0.000000000
Mo      1.824519667      0.896630777     -2.522222101

Then run the growth loop:

import os
from aegon.io.xyz                     import readxyzs, writexyzs
from aegon.population                 import rename, cutter_energy
from aegon.optimization.potentials.lj import opt_LJ_parallel
from aegon.discrimination.cluster_usr import deduplicate_by_usr
from growpal.libgrowpal               import growpal_parallel, display_info
from growpal.libselect_clustering     import select_by_clustering

nproc    = 2
growatom = 'Mo'
ecut     = 50.0

molecu0 = False
if __name__ == "__main__":
    for iii in range(10, 20 + 1):
        base_name = 'LJ' + str(iii).zfill(3)
        if os.path.isfile(base_name + '.xyz'):
            print("%s exists" % base_name)
            anterior = base_name + '.xyz'
            continue
        if not molecu0:
            molecu0 = readxyzs(anterior)
        molecu0 = rename(molecu0, base_name, 5)
        molecu1 = growpal_parallel(molecu0, growatom, dtol=1.2, n_cores=nproc)
        # Optimize and filter duplicates
        molecu0 = opt_LJ_parallel(molecu1, n_jobs=nproc)
        molecu1 = deduplicate_by_usr(molecu0, tols=0.99, tole=0.1, mono=True)
        molecu1 = cutter_energy(molecu1, ecut)
        molecu0 = select_by_clustering(molecu1, selection_count=round(iii * 1.8), n_clusters=8)
        display_info(molecu0[0:5], base_name)
        writexyzs(molecu0, base_name + '.xyz')

4. Verify LJ results against the reference database

from aegon.optimization.potentials.lj import lj_energy
from aegon.data.lj import get_lj_cluster
from aegon.io.xyz import readxyzs
import os

ii = 10
file = 'LJ' + str(ii).zfill(3) + '.xyz'
while os.path.isfile(file):
    mol    = readxyzs(file)
    nm     = len(mol)
    efili  = lj_energy(mol[0].get_positions())
    etrue  = get_lj_cluster(ii).info['e']
    deltae = efili - etrue
    print('%s (%s) %11.6f %5.2f' % (file, f"{nm:02d}", efili, deltae))
    ii  += 1
    file = 'LJ' + str(ii).zfill(3) + '.xyz'

5. Full production search (LJ010 → LJ100)

This is a much heavier run than examples 3–4: 90 sizes instead of 10, five re-optimization rounds per size instead of zero, and no upper bound on pool growth until the energy cutoff kicks in. Budget several hours and run it on a machine with several CPU cores and enough RAM for hundreds of concurrent local optimizations — a laptop-class 2–4 core machine will be very slow. nproc below auto-detects available cores; on a shared or memory-constrained machine, lower it explicitly instead.

This is the actual script behind the LJ75–LJ130 results in Key Results below (growpal/examples/LJ/x_run_growpal.py), starting from the same LJ010.xyz seed as example 3. Beyond just growing the pool at each size, it tracks a separate "special" lineage — the running best structure at each size is carried forward and re-injected into growth even if it wouldn't otherwise survive clustering-based selection, and a safety check after every size confirms the true global minimum (get_lj_cluster) was actually recovered before continuing, stopping early if not:

import os
import gc
from joblib.externals.loky        import get_reusable_executor
from aegon.io.xyz                 import readxyzs, writexyzs
from aegon.population             import rename, cutter_energy, sort_by_energy
from growpal.libgrowpal           import growpal_parallel
from growpal.libselect_clustering import select_by_clustering
from aegon.discrimination.cluster_usr import deduplicate_by_usr, filter_against_reference_usr
from aegon.optimization.potentials.lj import opt_LJ_parallel, lj_energy
from aegon.data.lj                import get_lj_cluster

nproc    = os.cpu_count() - 2
growatom = 'Mo'

molecu0      = False
special_mols = []
prev_comp    = 'LJ010.xyz'
prev_spec    = None
if __name__ == "__main__":
    for iii in range(11, 100 + 1):
        base_name = 'LJ' + str(iii).zfill(3)
        spec_file = base_name + '.xyz'
        comp_file = 'LJ_' + str(iii).zfill(3) + '_complement.xyz'
        if os.path.isfile(spec_file) and os.path.isfile(comp_file):
            print('%s exists' % base_name)
            prev_spec = spec_file
            prev_comp = comp_file
            molecu0 = False
            continue
        if not molecu0:
            molecu0      = readxyzs(prev_comp)
            special_mols = readxyzs(prev_spec) if prev_spec else []
        print('-------------------------------------------------------------------')
        print('%s START' % base_name)

        # Filter specials to the current size and tag them before growth
        special_mols = [m for m in special_mols if len(m) == iii - 1]
        sp_tags = []
        for i, sp in enumerate(special_mols):
            sp.info['i'] = f'_sp{i:03d}'
            sp_tags.append(f'_sp{i:03d}')
        molecu1 = rename(molecu0, base_name, 5)

        # Grow the regular pool and the special lineage together
        molecu0 = growpal_parallel(molecu1 + special_mols, growatom, dtol=1.2, n_cores=nproc)
        del molecu1
        molecu1 = opt_LJ_parallel(molecu0, nproc)
        del molecu0

        special_candidates = []
        for tag in sp_tags:
            children = [m for m in molecu1 if m.info['i'].startswith(tag + '_')]
            if children:
                best = min(children, key=lambda a: a.info['e'])
                special_candidates.append(best)
        special_mols = deduplicate_by_usr(special_candidates, tols=0.99, tole=0.1, mono=True)

        molecu0 = deduplicate_by_usr(molecu1, tols=0.99, tole=0.1, mono=True)
        del molecu1
        for counter in range(5):
            molecu1 = opt_LJ_parallel(molecu0, nproc)
            del molecu0
            molecu0 = deduplicate_by_usr(molecu1, tols=0.99, tole=0.1, mono=True)
            del molecu1
            gc.collect()

        # The running global minimum is guaranteed to end up in special_mols[0]
        global_min = molecu0[0]
        ans_type = filter_against_reference_usr([global_min], special_mols, tols=0.99, tole=0.1, mono=True, flag=1)
        if ans_type != []:
            special_mols.append(global_min)
            special_mols = sort_by_energy(special_mols, opt=1)
        print('Special cluster  : %d members (n=%d)' % (len(special_mols), iii))
        writexyzs(special_mols, spec_file)

        # Safety check: confirm the true global minimum was actually found
        e_best  = lj_energy(special_mols[0].get_positions())
        e_true  = get_lj_cluster(iii).info['e']
        delta_e = e_best - e_true
        print('Global min check: E_best=%.6f  E_true=%.6f  dE=%.4f' % (e_best, e_true, delta_e))
        if delta_e > 0.01:
            print('STOP: global minimum not found for %s (dE=%.4f > 0.01 eV)' % (base_name, delta_e))
            print('-------------------------------------------------------------------')
            break

        # Selection on the main pool for the next iteration's seed
        main_pool = filter_against_reference_usr(molecu0, special_mols, tols=0.99, tole=0.1, mono=True, flag=1)
        del molecu0
        ecut_n    = 5.0 * (iii ** (2.0 / 3.0))
        molecu1   = cutter_energy(main_pool, ecut_n)
        del main_pool
        nclusters = max(1, len(special_mols))
        sel_count = round(iii * 1.8)
        molecu0   = select_by_clustering(molecu1, selection_count=sel_count, n_clusters=nclusters,
                                         n_jobs=nproc, use_connectivity=True, mono=True)
        del molecu1
        gc.collect()
        writexyzs(molecu0, comp_file)
    get_reusable_executor().shutdown(wait=True)

Key Results

Lennard-Jones clusters (10–130 atoms)

• All reported global minima recovered, including the challenging LJ75–LJ78, LJ98 (tetrahedral, Td symmetry), and LJ102–LJ104 cases.
• Recovery of LJ98 required 1288 seed structures; genealogical analysis shows that its growth pathway passes through an intermediate at LJ55 that lies 8.4 eV above the putative global minimum — a structure that would be discarded by any purely energy-based selection.
• The selection_count parameter controls search thoroughness: larger values explore harder cases at the cost of more local optimizations.
• Efficiency comparison for LJ41: GrowPAL requires ~12,000 local optimizations vs. ~34,000 for SCG and ~118,000 for SFAEA (LJ38).

Sutton-Chen clusters (5–90 atoms, 10 metals)

• Nine parameter sets benchmarked (Al, Ag, Au, Cu, Ir, Ni, Pb, Pd, Pt, Rh); all previously reported global minima reproduced.
• New putative global minima identified for:
  • Palladium (n–m = 12-7): N = 26 (D3h, –0.092 eV vs. prior), 29 (Cs, –0.181 eV vs. prior), 37 (C2v), and 62 (Cs)
  • Aluminum (n–m = 7-6): N = 12 and 20
• First complete datasets of putative global minima provided for the unexplored n–m = 10-7 and 14-6 parameterizations (5–90 atoms).

Citation

If you use GrowPAL in your research, please cite:

López-Castro C., Ortiz-Chi F., Merino G. An Efficient Growth Pattern Algorithm (GrowPAL) for Cluster Structure Prediction. J. Chem. Theory Comput. 2024, 20, 4939–4948.

Authors

• Isaac Gutiérrez-Campos — Departamento de Física Aplicada, Cinvestav-IPN, Mérida, México
• Gabriel Merino — Departamento de Física Aplicada, Cinvestav-IPN, Mérida, México
• Filiberto Ortiz-Chi — Secihti-Departamento de Física Aplicada, Cinvestav-IPN, Mérida, México

License

MIT