Integrating Satisfiability Modulo Theories (SMT) solvers with large language models to enforce strict physical and chemical stoichiometry constraints reduces hallucinated compositions in generative inorganic crystal discovery.

Constraining an LLM's generative sampling process for inorganic crystal compositions with an SMT solver (encoding charge-neutrality, oxidation-state feasibility, stoichiometric ratio validity, and ionic radius/Pauling-rule constraints) reduces the rate of chemically invalid ("hallucinated") compositions from a measured baseline (unconstrained LLM) by a statistically significant margin, without reducing the rate of novel, synthesizable, property-relevant compositions among the valid outputs. Specifically: for a fixed LLM (e.g., a fine-tuned GPT-family or Llama-family model prompted/fine-tuned for crystal composition generation), the SMT-constrained pipeline will achieve ≥90% chemical validity (charge-balance + oxidation-state consistency, verified against ICSD/Materials Project ground truth and DFT-checked subsample) versus a baseline unconstrained validity rate of ≤60%, while retaining ≥70% of the unconstrained model's diversity (measured via unique composition count and Shannon entropy over element combinations) and ≥50% overlap with DFT-confirmed formation-energy-negative (stable/metastable, Ehull ≤ 50 meV/atom) structures relative to unconstrained sampling.

Minimum viable test (MVT):

1. Fine-tune or prompt-engineer one open LLM (Llama-3-8B or similar) on a composition-generation task using Materials Project (MP) training split (~80%).
2. Generate n=1,000 candidate compositions unconstrained (baseline arm) and n=1,000 with SMT-solver-in-the-loop constrained decoding (Z3 or CVC5 encoding charge balance + oxidation state tables) (treatment arm), using identical seeds/prompts.
3. Score validity: (a) rule-based — charge neutrality + oxidation state lookup against Materials Project ICSD; (b) DFT-based — random subsample of n=100 per arm relaxed via VASP/Quantum ESPRESSO or fast surrogate (M3GNet/CHGNet) to compute Ehull.
4. Compare validity %, diversity metrics, and DFT-confirmed stability rate between arms with pre-registered statistical tests.
5. Replicate on second LLM backbone and second held-out dataset (OQMD) to test generalization.

# Pseudocode: SMT-Constrained Crystal Composition Generation

def build_smt_constraints(elements, oxidation_states_db):
    solver = z3.Solver()
    counts = {e: z3.Int(f"count_{e}") for e in elements}
    ox_states = {e: z3.Int(f"ox_{e}") for e in elements}
    for e in elements:
        solver.add(counts[e] >= 0)
        solver.add(z3.Or([ox_states[e] == v for v in oxidation_states_db[e]]))
    # charge neutrality
    solver.add(z3.Sum([counts[e] * ox_states[e] for e in elements]) == 0)
    # stoichiometry bounds, radius-ratio rules appended similarly
    return solver, counts, ox_states

def llm_generate_candidate(prompt, model, temperature=0.8):
    return model.generate(prompt, temperature=temperature)  # returns composition string

def constrained_generation_loop(model, prompt, max_retries=5):
    for attempt in range(max_retries):
        candidate = llm_generate_candidate(prompt, model)
        elements, ratios = parse_composition(candidate)
        solver, counts, ox_states = build_smt_constraints(elements, OXIDATION_DB)
        # fix candidate ratios as additional constraints, check satisfiability
        solver.push()
        for e in elements:
            solver.add(counts[e] == ratios[e])
        if solver.check() == z3.sat:
            return candidate, "valid"
        solver.pop()
        # optional: ask solver for nearest satisfiable assignment (MaxSMT/soft constraints)
        prompt = augment_prompt_with_feedback(prompt, candidate, reason="charge_imbalance")
    return None, "rejected_after_max_retries"

def run_pipeline(n_candidates, model, dataset_split):
    results = []
    for i in range(n_candidates):
        candidate, status = constrained_generation_loop(model, base_prompt(dataset_split))
        results.append((candidate, status))
    return results

def evaluate(results, ground_truth_db, dft_surrogate):
    validity_rate = fraction_valid(results, ground_truth_db)
    subsample = random.sample([r for r in results if r[1]=="valid"], 100)
    ehull_scores = [dft_surrogate.relax_and_score(c) for c, _ in subsample]
    diversity = shannon_entropy(unique_element_sets(results))
    return validity_rate, ehull_scores, diversity