#!/usr/bin/env python3 """Periodic carbon test sites: wires, sheets, bulk, defects and confinement. Prepare one case (the default is geometry generation only): python 40_extended_systems.py --system cnt --guest water --output /path/to/run python 40_extended_systems.py --system graphene --vacancy --guest water \ --evaluate --method mace --output /path/to/run These are independently generated starting structures, not the optimized geometries of a literature benchmark. See the extended-systems acceptance protocol in docs/development/extended_systems_acceptance.md. A periodic axial scan is a potential-energy probe, not a diffusion coefficient or a demonstration of trapping. Vacancy structures require reconstruction and electronic/spin-state checks before scientific interpretation. """ from __future__ import annotations import argparse import json from pathlib import Path import numpy as np from ase.build import bulk, graphene, molecule, nanotube BOHR = 0.529177210903 SYSTEMS = ('cnt', 'graphene', 'diamond', 'slit') def to_periodic_system(atoms): """Map ASE row vectors to columns, putting periodic vectors first. Cartesian atom coordinates stay fixed. The legacy ASE SCF bridge only accepts 3D, whereas the direct MACE and semiempirical routes support lower-dimensional PeriodicSystem inputs. """ from vibeqc import Atom, PeriodicSystem periodic = np.flatnonzero(atoms.pbc) if not len(periodic): raise ValueError('A periodic system needs at least one periodic axis') order = np.r_[periodic, np.flatnonzero(~atoms.pbc)] lattice = atoms.cell.array[order].T/BOHR return PeriodicSystem(len(periodic), lattice, [Atom(int(z), p/BOHR) for z, p in zip(atoms.numbers, atoms.positions)], charge=0, multiplicity=1) def periodic_gradient(method, system, result): """Reuse the converged analytic result where its API exposes the kernel.""" if not getattr(result, 'converged', True): raise RuntimeError('An unconverged result cannot supply accepted forces') if method == 'mace': return np.asarray(result.gradient()) from vibeqc._vibeqc_core import semiempirical as se if method == 'gfn2_xtb': from vibeqc.semiempirical.methods.gfn2_params import load_gfn2_params return np.asarray(se.compute_periodic_gfn2_gradient( system, result.native_result, load_gfn2_params())) if method == 'dftb0': from vibeqc.semiempirical.parameters import default_parameters return np.asarray(se.compute_periodic_dftb0_gradient( system, result, default_parameters())) # SCC-DFTB currently needs 6N displaced SCFs. This is reached only with # explicit --forces, never during an energy-only size probe. from vibeqc.semiempirical.periodic import evaluate_periodic_energy_gradient energy, gradient = evaluate_periodic_energy_gradient(method, system) if abs(energy-float(result.energy)) > 1e-8: raise RuntimeError('Derivative evaluation changed the central energy; audit the SCF states') return np.asarray(gradient) def build_site(system='cnt', *, repeat=1, chirality=(10, 0), vacancy=False, guest='none', phase=0.0, radial_offset=0.0): """Return host, guest and combined ASE geometries with frozen fragments. CNT repetition grows the length; sheets grow in both in-plane directions; diamond grows in all three directions. Exactly one guest and at most one vacancy are inserted after repetition, so coverage changes with size. """ if system not in SYSTEMS or type(repeat) is not int or repeat < 1: raise ValueError('Choose a known system and a positive integer repeat') if guest not in ('none', 'water', 'methane'): raise ValueError('guest must be none, water or methane') if not np.isfinite([phase, radial_offset]).all(): raise ValueError('phase and radial_offset must be finite') if system == 'diamond' and (guest != 'none' or phase or radial_offset): raise ValueError('Bulk diamond probe has no guest or transport coordinate') if system != 'cnt' and (phase or radial_offset): raise ValueError('phase and radial_offset apply only to the nanotube') if system == 'cnt': n, m = chirality if not all(type(x) is int for x in chirality) or not n >= m >= 0 or n == 0: raise ValueError('chirality must be integers n >= m >= 0, with n > 0') host = nanotube(n, m, length=2*repeat, bond=1.42) host.center(vacuum=10.0, axis=(0, 1)) origin = np.array([host.cell[0, 0]/2, host.cell[1, 1]/2, host.cell[2, 2]*(0.5+phase)]) origin[0] += radial_offset pbc = (False, False, True) elif system in ('graphene', 'slit'): host = graphene(a=2.46, size=(3*repeat, 3*repeat, 1), vacuum=10.0) pbc = (True, True, False) center = (host.cell[0]+host.cell[1])/2 origin = np.array([center[0], center[1], host.positions[0, 2]+3.37]) if system == 'slit': upper = host.copy(); upper.positions[:, 2] += 6.8 host += upper host.cell[2, 2] += 6.8 origin[2] = host.positions[:, 2].min()+3.4 else: host = bulk('C', 'diamond', a=3.567, cubic=True)*((repeat,)*3) pbc = (True, True, True) origin = np.zeros(3) host.pbc = pbc removed = None if vacancy: # For a slit choose the lower sheet. Keep its site coordinate for the # paired adsorbed case; do not pretend the unreconstructed hole is stable. candidates = np.arange(len(host)) if system == 'slit': candidates = candidates[host.positions[:, 2] < origin[2]] target = origin.copy() if system == 'cnt': target[0] = host.positions[:, 0].max() target[2] = 0.5*host.cell[2, 2] i = min(candidates, key=lambda j: np.linalg.norm(host.positions[j]-target)) removed = host.positions[i].copy().tolist() if system == 'graphene': origin[:2] = host.positions[i, :2] del host[i] if guest == 'none': from ase import Atoms adsorbate = Atoms(cell=host.cell) else: adsorbate = molecule('H2O' if guest == 'water' else 'CH4') anchor = 0 # ASE H2O/CH4 places O/C first. adsorbate.translate(origin-adsorbate.positions[anchor]) adsorbate.cell = host.cell combined = host+adsorbate combined.pbc = pbc return {'combined': combined, 'host': host, 'guest': adsorbate}, { 'system': system, 'repeat': repeat, 'chirality': list(chirality) if system == 'cnt' else None, 'vacancy': vacancy, 'removed_position_angstrom': removed, 'guest': guest, 'guest_indices': list(range(len(host), len(combined))), 'phase_cell_lengths': phase, 'radial_offset_angstrom': radial_offset, 'pbc': list(pbc), 'dimension': sum(pbc), 'charge': 0, 'multiplicity': 1, 'relaxed': False, 'electronic_ground_state_established': False, 'literature_geometry': False, } def main(): p = argparse.ArgumentParser(description=__doc__) p.add_argument('--system', choices=SYSTEMS, default='cnt') p.add_argument('--repeat', type=int, default=1) p.add_argument('--chirality', nargs=2, type=int, default=(10, 0)) p.add_argument('--vacancy', action='store_true') p.add_argument('--guest', choices=['none', 'water', 'methane'], default='none') p.add_argument('--phase', type=float, default=0.0) p.add_argument('--radial-offset', type=float, default=0.0) p.add_argument('--component', choices=['combined', 'host', 'guest'], default='combined') p.add_argument('--method', choices=['mace', 'gfn2_xtb', 'dftb0', 'scc_dftb', 'seccm'], default='mace') p.add_argument('--weight-sha256') p.add_argument('--evaluate', action='store_true') p.add_argument('--forces', action='store_true', help='Also write forces; periodic SCC-DFTB needs 6N displaced SCFs') p.add_argument('--output', type=Path, required=True) args = p.parse_args() if args.weight_sha256 and args.method != 'mace': p.error('--weight-sha256 requires mace') if args.forces and not args.evaluate: p.error('--forces requires --evaluate') parts, record = build_site(args.system, repeat=args.repeat, chirality=tuple(args.chirality), vacancy=args.vacancy, guest=args.guest, phase=args.phase, radial_offset=args.radial_offset) from vibeqc import Atom, Molecule, run_job, run_periodic_job from vibeqc.output import OutputChannel, write from vibeqc.output.formats.xyz import write_xyz def mol(atoms): return Molecule([Atom(int(z), xyz/BOHR) for z, xyz in zip(atoms.numbers, atoms.positions)], 0, 1) args.output.mkdir(parents=True, exist_ok=True) for name, atoms in parts.items(): if not len(atoms): continue lattice = ' '.join(f'{v:.15g}' for v in atoms.cell.array.ravel()) flags = ' '.join('T' if b else 'F' for b in atoms.pbc) write_xyz(args.output/name, mol(atoms), comment=f'Lattice="{lattice}" Properties=species:S:1:pos:R:3 pbc="{flags}" unrelaxed=true') write(json.dumps(record, sort_keys=True)+'\n') with OutputChannel.to_file(args.output/'site.json'): write(json.dumps(record, indent=2)+'\n') if not args.evaluate: return atoms = parts[args.component] if not len(atoms): p.error('The requested guest component is empty') stem = args.output/(args.component+'-result') if args.method == 'seccm': from vibeqc.semiempirical.methods.msindo_ccm import CCMOptions options = CCMOptions(translations=atoms.cell.array[atoms.pbc].tolist()) if not atoms.pbc.any(): p.error('Use the molecular MSINDO example for an isolated guest') result = run_job(mol(atoms), method='seccm', ccm_options=options, output=stem, record_hostname=False) else: options = {} if args.method == 'mace': from vibeqc.mlip import MLIPOptions options['mlip_options'] = MLIPOptions(model='medium-mpa-0', device='cpu', dtype='float64', weight_sha256=args.weight_sha256) if atoms.pbc.any(): system = to_periodic_system(atoms) result = run_periodic_job(system, None, method=args.method, output=stem, record_hostname=False, **options) else: result = run_job(mol(atoms), method=args.method, output=stem, record_hostname=False, **options) summary = {'energy_ha': float(result.energy), 'converged': bool(getattr(result, 'converged', True)), 'accepted_equilibrium_geometry': False} if args.forces: gradient = (periodic_gradient(args.method, system, result) if atoms.pbc.any() and args.method != 'seccm' else np.asarray(result.gradient())) summary.update(forces_ha_per_bohr=(-gradient).tolist(), max_force_ha_per_bohr=float(np.linalg.norm(gradient, axis=1).max()), gradient_potential=('mermin_free_energy' if args.method == 'gfn2_xtb' and atoms.pbc.any() else 'energy')) if hasattr(result, 'free_energy'): summary['free_energy_ha'] = float(result.free_energy) with OutputChannel.to_file(args.output/(args.component+'-summary.json')): write(json.dumps(summary, indent=2, allow_nan=False)+'\n') write(json.dumps(summary, allow_nan=False)+'\n') if __name__ == '__main__': main()