diff --git a/pyproject.toml b/pyproject.toml index eacd84d..711c5bd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,7 @@ mlip = ["mace-torch"] plot = ["matplotlib>=3.5"] pacmof2 = ["pacmof2"] graspa = ["pyyaml>=6.0"] +pygraspa = ["pyyaml>=6.0"] all = ["rdkit", "mace-torch", "matplotlib>=3.5", "pacmof2", "pyyaml>=6.0"] dev = ["pytest>=7.0", "ruff>=0.4"] @@ -51,4 +52,7 @@ include-package-data = true "matkit.graspa.files.template_mixture" = ["*"] "matkit.graspa.files.template_co2_n2" = ["*"] "matkit.graspa.files.template_mixture_isotherm" = ["*"] +"matkit.pygraspa.files.template" = ["*"] +"matkit.pygraspa.files.template_mixture" = ["*"] +"matkit.pygraspa.files.template_mixture_isotherm" = ["*"] "matkit.zeopp.files" = ["*"] diff --git a/scripts/pygraspa/isotherm_config.yaml b/scripts/pygraspa/isotherm_config.yaml new file mode 100644 index 0000000..eb6b54f --- /dev/null +++ b/scripts/pygraspa/isotherm_config.yaml @@ -0,0 +1,42 @@ +# pygRASPA isotherm simulation config (ML-potential GCMC) +# Usage: python setup_isotherms.py isotherm_config.yaml + +cif_dir: ./cifs +outdir: ./batch_results + +# mode: single — each adsorbate runs independently (separate sims) +# mode: mixture — all adsorbates run together in one sim (requires MolFraction) +mode: single + +adsorbates: [CO2, N2] + +# ML configuration (required for pygRASPA). +# ecomps must contain one entry per adsorbate name above. +# Compute these once with `matkit pygraspa compute-ecomp` and paste here. +ml: + model_path: /path/to/model.pt + model_type: FAIRChem-esen # FAIRChem-esen | FAIRChem-uma | FAIRChem-AllScAIP | mace_polar | Allegro + task: null # required for FAIRChem-uma / FAIRChem-AllScAIP + run_mode: run-auto # run | run-auto + save_poscar: false + ecomps: + CO2: -22.97915 + N2: -17.5 + # H2O: -14.2 + # CH4: -11.3 + +temperatures: [298] + +pressures: + - 0.001 + - 0.01 + - 0.05 + - 0.1 + - 0.5 + - 1.0 + +pressure_unit: bar # bar, kPa, atm, or Pa + +cutoff: 12.8 +cycles: 50000 +max_workers: 4 # parallel threads for CIF processing (null = auto) diff --git a/scripts/pygraspa/setup_isotherms.py b/scripts/pygraspa/setup_isotherms.py new file mode 100644 index 0000000..a4d2ab4 --- /dev/null +++ b/scripts/pygraspa/setup_isotherms.py @@ -0,0 +1,140 @@ +"""Set up pygRASPA isotherm simulations (ML-potential GCMC) from YAML. + +Usage: + python setup_isotherms.py + python setup_isotherms.py # defaults to isotherm_config.yaml + +Supports two modes via the 'mode' field in the YAML config: + - single: each adsorbate runs as independent single-component sims + - mixture: all adsorbates run together in one sim (requires MolFraction) +""" + +import sys +from pathlib import Path + +import yaml + +from matkit.pygraspa import setup_batch + +PRESSURE_TO_PA = { + "pa": 1.0, + "bar": 1e5, + "kpa": 1e3, + "atm": 101325.0, +} + + +def main(): + if len(sys.argv) > 1: + config_path = Path(sys.argv[1]) + else: + config_path = Path(__file__).parent / "isotherm_config.yaml" + + with open(config_path) as f: + cfg = yaml.safe_load(f) + + pressure_unit = cfg.get("pressure_unit", "Pa").lower() + if pressure_unit not in PRESSURE_TO_PA: + raise ValueError( + f"Unknown pressure_unit '{cfg['pressure_unit']}'. " + f"Supported: {list(PRESSURE_TO_PA.keys())}" + ) + factor = PRESSURE_TO_PA[pressure_unit] + pressures_pa = [p * factor for p in cfg["pressures"]] + + mode = cfg.get("mode", "single") + ml = cfg["ml"] + model_path = ml["model_path"] + model_type = ml.get("model_type", "FAIRChem-esen") + task = ml.get("task") + ecomp_lookup = ml["ecomps"] # {adsorbate_name: float} + run_mode = ml.get("run_mode", "run-auto") + save_poscar = ml.get("save_poscar", False) + + print(f"Config: {config_path}") + print(f"Mode: {mode}") + print(f"ML model: {model_path} ({model_type})") + print(f"CIF dir: {cfg['cif_dir']}") + print(f"Temperatures: {cfg['temperatures']}") + print( + f"Pressures: {len(pressures_pa)} points " + f"({cfg['pressures'][0]}-{cfg['pressures'][-1]} " + f"{cfg.get('pressure_unit', 'Pa')})" + ) + + common_kwargs = dict( + cif_dir=cfg["cif_dir"], + temperatures=cfg["temperatures"], + pressures=pressures_pa, + cutoff=cfg.get("cutoff", 12.8), + n_cycle=cfg.get("cycles", 1000), + max_workers=cfg.get("max_workers"), + model_path=model_path, + model_type=model_type, + task=task, + mode=run_mode, + save_poscar=save_poscar, + ) + + if mode == "single": + print(f"Adsorbates: {cfg['adsorbates']} (each runs independently)") + total = 0 + for ads_name in cfg["adsorbates"]: + if ads_name not in ecomp_lookup: + raise ValueError( + f"E_comp for adsorbate {ads_name!r} not provided " + "under ml.ecomps in the YAML config." + ) + ads_outdir = str(Path(cfg["outdir"]) / ads_name) + adsorbates = [{"MoleculeName": ads_name}] + manifest = setup_batch( + outpath=ads_outdir, + adsorbates=adsorbates, + E_comps=[ecomp_lookup[ads_name]], + template_dir=cfg.get("template_dir", "template"), + **common_kwargs, + ) + print( + f" {ads_name}: {len(manifest)} simulations in {ads_outdir}" + ) + total += len(manifest) + print(f"\nTotal: {total} simulations") + + elif mode == "mixture": + adsorbates = [] + for ad in cfg["adsorbates"]: + entry = {"MoleculeName": ad["name"]} + for k, v in ad.items(): + if k != "name": + entry[k] = v + adsorbates.append(entry) + names = [ad["MoleculeName"] for ad in adsorbates] + print(f"Adsorbates: {names} (mixture)") + + ecomps = [] + for name in names: + if name not in ecomp_lookup: + raise ValueError( + f"E_comp for adsorbate {name!r} not provided " + "under ml.ecomps." + ) + ecomps.append(ecomp_lookup[name]) + + template_dir = cfg.get("template_dir", "template_mixture_isotherm") + manifest = setup_batch( + outpath=cfg["outdir"], + adsorbates=adsorbates, + E_comps=ecomps, + template_dir=template_dir, + **common_kwargs, + ) + print(f"\nSet up {len(manifest)} simulations in {cfg['outdir']}") + + else: + raise ValueError(f"Unknown mode '{mode}'. Use 'single' or 'mixture'.") + + print(f"Manifest: {cfg['outdir']}/simulations.jsonl") + + +if __name__ == "__main__": + main() diff --git a/src/matkit/cli.py b/src/matkit/cli.py index 7eaa788..6449d4b 100644 --- a/src/matkit/cli.py +++ b/src/matkit/cli.py @@ -246,6 +246,381 @@ def graspa_sycl_analyze(path, unit, adsorbate, fname): click.echo(f"Error analyzing gRASPA SYCL results: {e}", err=True) +# ========================================== +# PYGRASPA COMMANDS (ML-potential GCMC) +# ========================================== +@main.group("pygraspa") +def pygraspa_cli(): + """Commands for pygRASPA simulations (ML-potential GCMC). + + Writes simulation files plus a ``run.py`` launcher. Execute on a GPU + machine with pygRASPA + the required ML backend installed. + """ + pass + + +def _parse_ecomps(ctx, param, value): + """Click callback: convert comma-separated floats into a list.""" + if value is None: + return None + parts = [p.strip() for p in value.split(",") if p.strip()] + try: + return [float(p) for p in parts] + except ValueError as e: + raise click.BadParameter( + f"--ecomps must be comma-separated floats: {e}" + ) + + +@pygraspa_cli.command("setup") +@click.option( + "--cif", + required=True, + type=click.Path(exists=True), + help="Path to input CIF file.", +) +@click.option( + "--outdir", + required=True, + type=click.Path(), + help="Directory to generate simulation files.", +) +@click.option( + "--adsorbate", + required=True, + multiple=True, + help="Adsorbate molecule name (e.g. CO2). Can be repeated.", +) +@click.option( + "--model", + "model_path", + required=True, + help="ML model checkpoint path.", +) +@click.option( + "--model-type", + default="FAIRChem-esen", + type=click.Choice( + [ + "FAIRChem-esen", + "FAIRChem-uma", + "FAIRChem-AllScAIP", + "mace_polar", + "Allegro", + "customized", + ] + ), + help="ML backend identifier.", +) +@click.option( + "--task", + default=None, + help="Task name (FAIRChem-uma / AllScAIP).", +) +@click.option( + "--ecomps", + required=True, + callback=_parse_ecomps, + help="Comma-separated E_comp values (eV), one per --adsorbate.", +) +@click.option( + "--mode", + default="run-auto", + type=click.Choice(["run", "run-auto"]), + help="pygRASPA execution mode written into run.py.", +) +@click.option("--save-poscar", is_flag=True, help="Save accepted POSCARs.") +@click.option("--temp", default=298.0, help="Temperature in Kelvin.") +@click.option("--pressure", default=1e5, help="Pressure in Pa.") +@click.option("--cutoff", default=12.8, help="Cutoff radius (Angstrom).") +@click.option("--cycles", default=1000, help="Number of MC cycles.") +@click.option( + "--template-dir", + default="template", + help=( + "Template subdir (template, template_mixture, " + "template_mixture_isotherm)." + ), +) +def pygraspa_setup( + cif, + outdir, + adsorbate, + model_path, + model_type, + task, + ecomps, + mode, + save_poscar, + temp, + pressure, + cutoff, + cycles, + template_dir, +): + """Setup a pygRASPA simulation.""" + from matkit.pygraspa import setup_simulation + + adsorbate_list = [{"MoleculeName": ads} for ads in adsorbate] + try: + setup_simulation( + cif=cif, + outpath=outdir, + adsorbates=adsorbate_list, + model_path=model_path, + model_type=model_type, + E_comps=ecomps, + task=task, + mode=mode, + save_poscar=save_poscar, + temperature=temp, + pressure=pressure, + cutoff=cutoff, + n_cycle=cycles, + template_dir=template_dir, + ) + click.echo(f"Successfully set up pygRASPA simulation in {outdir}") + except Exception as e: + click.echo(f"Error setting up pygRASPA simulation: {e}", err=True) + + +@pygraspa_cli.command("batch-setup") +@click.option( + "--cif-dir", + required=True, + type=click.Path(exists=True, file_okay=False), + help="Directory containing CIF files.", +) +@click.option( + "--outdir", + required=True, + type=click.Path(), + help="Base output directory.", +) +@click.option( + "--adsorbate", + required=True, + multiple=True, + help="Adsorbate molecule name. Can be repeated.", +) +@click.option( + "--model", + "model_path", + required=True, + help="ML model checkpoint path.", +) +@click.option( + "--model-type", + default="FAIRChem-esen", + type=click.Choice( + [ + "FAIRChem-esen", + "FAIRChem-uma", + "FAIRChem-AllScAIP", + "mace_polar", + "Allegro", + "customized", + ] + ), + help="ML backend identifier.", +) +@click.option( + "--task", + default=None, + help="Task name (FAIRChem-uma / AllScAIP).", +) +@click.option( + "--ecomps", + required=True, + callback=_parse_ecomps, + help="Comma-separated E_comp values, one per --adsorbate.", +) +@click.option( + "--mode", + default="run-auto", + type=click.Choice(["run", "run-auto"]), + help="pygRASPA execution mode written into run.py.", +) +@click.option("--save-poscar", is_flag=True, help="Save accepted POSCARs.") +@click.option( + "--temp", + required=True, + multiple=True, + type=float, + help="Temperature(s) in Kelvin. Can be repeated.", +) +@click.option( + "--pressure", + required=True, + multiple=True, + type=float, + help="Pressure(s) in Pa. Can be repeated.", +) +@click.option("--cutoff", default=12.8, help="Cutoff radius (Angstrom).") +@click.option("--cycles", default=1000, help="Number of MC cycles.") +@click.option( + "--template-dir", + default="template", + help="Template subdir.", +) +@click.option( + "--workers", + default=None, + type=int, + help="Max parallel threads for CIF processing.", +) +def pygraspa_batch_setup( + cif_dir, + outdir, + adsorbate, + model_path, + model_type, + task, + ecomps, + mode, + save_poscar, + temp, + pressure, + cutoff, + cycles, + template_dir, + workers, +): + """Set up pygRASPA simulations for all CIF x T x P.""" + from matkit.pygraspa import setup_batch + + adsorbate_list = [{"MoleculeName": ads} for ads in adsorbate] + try: + manifest = setup_batch( + cif_dir=cif_dir, + outpath=outdir, + adsorbates=adsorbate_list, + model_path=model_path, + model_type=model_type, + E_comps=ecomps, + task=task, + mode=mode, + save_poscar=save_poscar, + temperatures=list(temp), + pressures=list(pressure), + cutoff=cutoff, + n_cycle=cycles, + template_dir=template_dir, + max_workers=workers, + ) + click.echo(f"Set up {len(manifest)} simulations in {outdir}") + click.echo(f"Manifest written to {outdir}/simulations.jsonl") + except Exception as e: + click.echo(f"Error setting up batch simulations: {e}", err=True) + + +@pygraspa_cli.command("compute-ecomp") +@click.option( + "--def-file", + "def_file", + required=True, + type=click.Path(exists=True), + help="Path to adsorbate .def file.", +) +@click.option( + "--model", + "model_path", + required=True, + help="ML model checkpoint path.", +) +@click.option( + "--model-type", + default="FAIRChem-esen", + type=click.Choice( + [ + "FAIRChem-esen", + "FAIRChem-uma", + "FAIRChem-AllScAIP", + "mace_polar", + "Allegro", + "customized", + ] + ), +) +@click.option( + "--task", + default=None, + help="Task name (FAIRChem-uma / AllScAIP).", +) +@click.option("--device", default="cuda", type=click.Choice(["cuda", "cpu"])) +@click.option( + "--cache", + "cache_path", + default=None, + type=click.Path(), + help="JSON cache file to memoize results.", +) +def pygraspa_compute_ecomp( + def_file, + model_path, + model_type, + task, + device, + cache_path, +): + """Compute E_comp (isolated-adsorbate ML energy) via pygRASPA.""" + from matkit.pygraspa import compute_ecomp + + try: + value = compute_ecomp( + adsorbate_def=def_file, + model_path=model_path, + model_type=model_type, + task=task, + device=device, + cache_path=cache_path, + ) + click.echo( + json.dumps( + {"adsorbate_def": def_file, "E_comp_eV": value}, + indent=2, + ) + ) + except Exception as e: + click.echo(f"Error computing E_comp: {e}", err=True) + + +@pygraspa_cli.command("analyze") +@click.option( + "--path", + required=True, + type=click.Path(exists=True), + help="Path to simulation output directory.", +) +@click.option( + "--unit", + default="mol/kg", + type=click.Choice(["mol/kg", "mg/g", "g/L"]), + help="Unit for uptake.", +) +@click.option("--fname", default="output.log", help="Log filename.") +@click.option( + "--skip-cycles", + default=None, + type=int, + help="Override n_skip_cycles (default: from simulation.input).", +) +def pygraspa_analyze(path, unit, fname, skip_cycles): + """Analyze pygRASPA simulation results.""" + from matkit.pygraspa import get_output_data + + try: + result = get_output_data( + output_path=path, + unit=unit, + output_fname=fname, + n_skip_cycles=skip_cycles, + ) + click.echo(json.dumps(result, indent=2)) + except Exception as e: + click.echo(f"Error analyzing pygRASPA results: {e}", err=True) + + # ========================================== # RASPA2 COMMANDS # ========================================== diff --git a/src/matkit/pygraspa/__init__.py b/src/matkit/pygraspa/__init__.py new file mode 100644 index 0000000..642b28c --- /dev/null +++ b/src/matkit/pygraspa/__init__.py @@ -0,0 +1,18 @@ +"""matkit pygRASPA backend (ML-potential GCMC, setup-only). + +See :mod:`matkit.pygraspa.pygraspa` for the public API. +""" + +from matkit.pygraspa.pygraspa import ( + compute_ecomp, + get_output_data, + setup_batch, + setup_simulation, +) + +__all__ = [ + "compute_ecomp", + "get_output_data", + "setup_batch", + "setup_simulation", +] diff --git a/src/matkit/pygraspa/files/template/CO2.def b/src/matkit/pygraspa/files/template/CO2.def new file mode 100644 index 0000000..fa1ff0a --- /dev/null +++ b/src/matkit/pygraspa/files/template/CO2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +304.1282 +7377300.0 +0.22394 +#Number Of Atoms + 3 +# Number of groups +1 +# CO2-group +rigid +# number of atoms +3 +# atomic positions +0 O_co2 0.0 0.0 1.16 +1 C_co2 0.0 0.0 0.0 +2 O_co2 0.0 0.0 -1.16 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template/H2.def b/src/matkit/pygraspa/files/template/H2.def new file mode 100644 index 0000000..89598c2 --- /dev/null +++ b/src/matkit/pygraspa/files/template/H2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +33.19 +1315000.0 +-0.214 +#Number Of Atoms +3 +# Number of groups +1 +# H2-group +rigid +# number of atoms +3 +# atomic positions +0 H_h2 0.0 0.0 0.37 +1 H_com 0.0 0.0 0.0 +2 H_h2 0.0 0.0 -0.37 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template/N2.def b/src/matkit/pygraspa/files/template/N2.def new file mode 100644 index 0000000..53068c5 --- /dev/null +++ b/src/matkit/pygraspa/files/template/N2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +126.192 +3395800.0 +0.0372 +#Number Of Atoms +3 +# Number of groups +1 +# N2-group +rigid +# number of atoms +3 +# atomic positions +0 N_n2 0.0 0.0 0.55 +1 N_com 0.0 0.0 0.0 +2 N_n2 0.0 0.0 -0.55 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template/TIP4P.def b/src/matkit/pygraspa/files/template/TIP4P.def new file mode 100644 index 0000000..b6b4fd6 --- /dev/null +++ b/src/matkit/pygraspa/files/template/TIP4P.def @@ -0,0 +1,25 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +647.14 +22064000.0 +-0.217000 +# total number Of atoms +4 +# Number of groups +1 +# water-group +rigid +# number of atoms +4 +# atomic positions +0 Ow 0.0 0.0 0.0 +1 Lw 0.0 0.15 0.0 +2 Hw 0.75695 0.58588 0.0 +3 Hw -0.75695 0.58588 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +0 2 RIGID_BOND +0 3 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template/force_field.def b/src/matkit/pygraspa/files/template/force_field.def new file mode 100644 index 0000000..0e7bb78 --- /dev/null +++ b/src/matkit/pygraspa/files/template/force_field.def @@ -0,0 +1,6 @@ +# rules to overwrite +0 +# number of defined interactions +0 +# mixing rules to overwrite +0 diff --git a/src/matkit/pygraspa/files/template/force_field_mixing_rules.def b/src/matkit/pygraspa/files/template/force_field_mixing_rules.def new file mode 100644 index 0000000..a062f88 --- /dev/null +++ b/src/matkit/pygraspa/files/template/force_field_mixing_rules.def @@ -0,0 +1,121 @@ +# general rule for shifted vs truncated +truncated +# general rule tailcorrections +no +# number of defined interactions +112 +# type interaction +Ac lennard-jones 16.60 3.10 +Ag lennard-jones 18.11 2.80 +Al lennard-jones 254.09 4.01 +Am lennard-jones 7.04 3.01 +Ar lennard-jones 93.08 3.45 +As lennard-jones 155.47 3.77 +At lennard-jones 142.89 4.23 +Au lennard-jones 19.62 2.93 +B lennard-jones 90.57 3.64 +Ba lennard-jones 183.15 3.30 +Be lennard-jones 42.77 2.45 +Bi lennard-jones 260.63 3.89 +Bk lennard-jones 6.54 2.97 +Br lennard-jones 126.30 3.73 +C lennard-jones 52.83 3.43 +Ca lennard-jones 119.75 3.03 +Cd lennard-jones 114.72 2.54 +Ce lennard-jones 6.54 3.17 +Cf lennard-jones 6.54 2.95 +Cl lennard-jones 114.23 3.517 +Cm lennard-jones 6.54 2.96 +Co lennard-jones 7.04 2.56 +Cr lennard-jones 7.55 2.69 +Cs lennard-jones 22.64 4.02 +Cu lennard-jones 2.516 3.114 +Dy lennard-jones 3.52 3.05 +Er lennard-jones 3.52 3.02 +Es lennard-jones 6.04 2.94 +Eu lennard-jones 4.03 3.11 +F lennard-jones 25.16 2.997 +Fe lennard-jones 6.54 2.59 +Fm lennard-jones 6.04 2.93 +Fr lennard-jones 25.16 4.37 +Ga lennard-jones 208.81 3.90 +Gd lennard-jones 4.53 3.00 +Ge lennard-jones 190.69 3.81 +H lennard-jones 22.14 2.57 +He lennard-jones 10.9 2.64 +Hf lennard-jones 36.23 2.80 +Hg lennard-jones 193.71 2.41 +Ho lennard-jones 3.52 3.04 +I lennard-jones 170.57 4.01 +In lennard-jones 301.39 3.98 +Ir lennard-jones 36.73 2.53 +K lennard-jones 17.61 3.40 +Kr lennard-jones 166.40 3.636 +La lennard-jones 8.55 3.14 +Li lennard-jones 12.58 2.18 +Lu lennard-jones 20.63 3.24 +Md lennard-jones 5.53 2.92 +Mg lennard-jones 55.85 2.69 +Mn lennard-jones 6.54 2.64 +Mo lennard-jones 28.18 2.72 +N lennard-jones 34.72 3.26 +Na lennard-jones 15.09 2.66 +Nb lennard-jones 29.69 2.82 +Nd lennard-jones 5.03 3.18 +Ne lennard-jones 21.13 2.66 +Ni lennard-jones 7.55 2.52 +No lennard-jones 5.53 2.89 +Np lennard-jones 9.56 3.05 +O lennard-jones 30.19 3.12 +Os lennard-jones 18.62 2.78 +P lennard-jones 153.46 3.69 +Pa lennard-jones 11.07 3.05 +Pb lennard-jones 333.59 3.83 +Pd lennard-jones 24.15 2.58 +Pm lennard-jones 4.53 3.16 +Po lennard-jones 163.52 4.20 +Pr lennard-jones 5.03 3.21 +Pt lennard-jones 40.25 2.45 +Pu lennard-jones 8.05 3.05 +Ra lennard-jones 203.27 3.28 +Rb lennard-jones 20.13 3.67 +Re lennard-jones 33.21 2.63 +Rh lennard-jones 26.67 2.61 +Rn lennard-jones 124.78 4.25 +Ru lennard-jones 28.18 2.64 +S lennard-jones 137.86 3.59 +Sb lennard-jones 225.91 3.94 +Sc lennard-jones 9.56 2.94 +Se lennard-jones 146.42 3.75 +Si lennard-jones 202.27 3.83 +Sm lennard-jones 4.03 3.14 +Sn lennard-jones 285.28 3.91 +Sr lennard-jones 118.24 3.24 +Ta lennard-jones 40.75 2.82 +Tb lennard-jones 3.52 3.07 +Tc lennard-jones 24.15 2.67 +Te lennard-jones 200.25 3.98 +Th lennard-jones 13.08 3.03 +Ti lennard-jones 8.55 2.83 +Tl lennard-jones 342.14 3.87 +Tm lennard-jones 3.02 3.01 +U lennard-jones 11.07 3.02 +V lennard-jones 8.05 2.80 +W lennard-jones 33.71 2.73 +Xe lennard-jones 221.00 4.10 +Y lennard-jones 36.23 2.98 +Yb lennard-jones 114.72 2.99 +Zn lennard-jones 62.40 2.46 +Zr lennard-jones 34.72 2.783 +C_co2 lennard-jones 27.0 2.80 +O_co2 lennard-jones 79.0 3.05 +N_n2 lennard-jones 36.0 3.31 +N_com lennard-jones 0.0 0.0 +Ow lennard-jones 78.0 3.154 +Hw lennard-jones 0.0 0.0 +Lw lennard-jones 0.0 0.0 +CH4_sp3 lennard-jones 148.0 3.73 +H_h2 lennard-jones 0.0 0.0 +H_com lennard-jones 36.7 2.958 +# general mixing rule for Lennard-Jones +Lorentz-Berthelot diff --git a/src/matkit/pygraspa/files/template/methane.def b/src/matkit/pygraspa/files/template/methane.def new file mode 100644 index 0000000..5b5aec7 --- /dev/null +++ b/src/matkit/pygraspa/files/template/methane.def @@ -0,0 +1,18 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +190.564 +4599200.0 +0.01142 +# Number Of Atoms +1 +# Number Of Groups +1 +# Alkane-group +rigid +# number of atoms +1 +# atomic positions +0 CH4_sp3 0.0 0.0 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Bond/Bend Bend/Bend Bond/Torsion Bend/Torsion IntraVDW Intra ch-ch Intra ch-bd Intra bd-bd + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Number of config moves +0 \ No newline at end of file diff --git a/src/matkit/pygraspa/files/template/pseudo_atoms.def b/src/matkit/pygraspa/files/template/pseudo_atoms.def new file mode 100644 index 0000000..7d35e12 --- /dev/null +++ b/src/matkit/pygraspa/files/template/pseudo_atoms.def @@ -0,0 +1,115 @@ +#number of pseudo atoms +112 +#type print as chem oxi mass charge polarization B-factor radii connectivity anisotropic anisotropic-type tinker-type +Ac yes Ac Ac 0 227 0 0 1 1 0 0 relative 0 +Ag yes Ag Ag 0 107.868 0 0 1 1 0 0 relative 0 +Al yes Al Al 0 26.982 0 0 1 1 0 0 relative 0 +Am yes Am Am 0 243 0 0 1 1 0 0 relative 0 +Ar yes Ar Ar 0 39.948 0 0 1 1 0 0 relative 0 +As yes As As 0 74.922 0 0 1 1 0 0 relative 0 +At yes At At 0 210 0 0 1 1 0 0 relative 0 +Au yes Au Au 0 196.967 0 0 1 1 0 0 relative 0 +B yes B B 0 10.811 0 0 1 1 0 0 relative 0 +Ba yes Ba Ba 0 137.327 0 0 1 1 0 0 relative 0 +Be yes Be Be 0 9.012 0 0 1 1 0 0 relative 0 +Bi yes Bi Bi 0 208.98 0 0 1 1 0 0 relative 0 +Bk yes Bk Bk 0 247 0 0 1 1 0 0 relative 0 +Br yes Br Br 0 79.904 0 0 1 1 0 0 relative 0 +C yes C C 0 12.011 0 0 1 1 0 0 relative 0 +Ca yes Ca Ca 0 40.078 0 0 1 1 0 0 relative 0 +Cd yes Cd Cd 0 112.411 0 0 1 1 0 0 relative 0 +Ce yes Ce Ce 0 140.116 0 0 1 1 0 0 relative 0 +Cf yes Cf Cf 0 251 0 0 1 1 0 0 relative 0 +Cl yes Cl Cl 0 35.453 0 0 1 1 0 0 relative 0 +Cm yes Cm Cm 0 247 0 0 1 1 0 0 relative 0 +Co yes Co Co 0 58.933 0 0 1 1 0 0 relative 0 +Cr yes Cr Cr 0 51.996 0 0 1 1 0 0 relative 0 +Cs yes Cs Cs 0 132.905 0 0 1 1 0 0 relative 0 +Cu yes Cu Cu 0 63.546 0 0 1 1 0 0 relative 0 +Dy yes Dy Dy 0 162.5 0 0 1 1 0 0 relative 0 +Er yes Er Er 0 167.26 0 0 1 1 0 0 relative 0 +Es yes Es Es 0 252 0 0 1 1 0 0 relative 0 +Eu yes Eu Eu 0 151.964 0 0 1 1 0 0 relative 0 +F yes F F 0 18.998 0 0 1 1 0 0 relative 0 +Fe yes Fe Fe 0 55.845 0 0 1 1 0 0 relative 0 +Fm yes Fm Fm 0 257 0 0 1 1 0 0 relative 0 +Fr yes Fr Fr 0 223 0 0 1 1 0 0 relative 0 +Ga yes Ga Ga 0 69.723 0 0 1 1 0 0 relative 0 +Gd yes Gd Gd 0 157.25 0 0 1 1 0 0 relative 0 +Ge yes Ge Ge 0 72.61 0 0 1 1 0 0 relative 0 +H yes H H 0 1.008 0 0 1 1 0 0 relative 0 +He yes He He 0 4.003 0 0 1 1 0 0 relative 0 +Hf yes Hf Hf 0 178.49 0 0 1 1 0 0 relative 0 +Hg yes Hg Hg 0 200.59 0 0 1 1 0 0 relative 0 +Ho yes Ho Ho 0 164.93 0 0 1 1 0 0 relative 0 +I yes I I 0 126.904 0 0 1 1 0 0 relative 0 +In yes In In 0 114.818 0 0 1 1 0 0 relative 0 +Ir yes Ir Ir 0 192.217 0 0 1 1 0 0 relative 0 +K yes K K 0 39.098 0 0 1 1 0 0 relative 0 +Kr yes Kr Kr 0 83.8 0 0 1 1 0 0 relative 0 +La yes La La 0 138.906 0 0 1 1 0 0 relative 0 +Li yes Li Li 0 6.941 0 0 1 1 0 0 relative 0 +Lu yes Lu Lu 0 174.967 0 0 1 1 0 0 relative 0 +Md yes Md Md 0 258 0 0 1 1 0 0 relative 0 +Mg yes Mg Mg 0 24.305 0 0 1 1 0 0 relative 0 +Mn yes Mn Mn 0 54.938 0 0 1 1 0 0 relative 0 +Mo yes Mo Mo 0 95.94 0 0 1 1 0 0 relative 0 +N yes N N 0 14.007 0 0 1 1 0 0 relative 0 +Na yes Na Na 0 22.991 0 0 1 1 0 0 relative 0 +Nb yes Nb Nb 0 92.906 0 0 1 1 0 0 relative 0 +Nd yes Nd Nd 0 144.24 0 0 1 1 0 0 relative 0 +Ne yes Ne Ne 0 20.18 0 0 1 1 0 0 relative 0 +Ni yes Ni Ni 0 58.693 0 0 1 1 0 0 relative 0 +No yes No No 0 259 0 0 1 1 0 0 relative 0 +Np yes Np Np 0 237 0 0 1 1 0 0 relative 0 +O yes O O 0 15.999 0 0 1 1 0 0 relative 0 +Os yes Os Os 0 190.23 0 0 1 1 0 0 relative 0 +P yes P P 0 30.974 0 0 1 1 0 0 relative 0 +Pa yes Pa Pa 0 231.036 0 0 1 1 0 0 relative 0 +Pb yes Pb Pb 0 207.2 0 0 1 1 0 0 relative 0 +Pd yes Pd Pd 0 106.42 0 0 1 1 0 0 relative 0 +Pm yes Pm Pm 0 145 0 0 1 1 0 0 relative 0 +Po yes Po Po 0 210 0 0 1 1 0 0 relative 0 +Pr yes Pr Pr 0 140.908 0 0 1 1 0 0 relative 0 +Pt yes Pt Pt 0 195.078 0 0 1 1 0 0 relative 0 +Pu yes Pu Pu 0 244 0 0 1 1 0 0 relative 0 +Ra yes Ra Ra 0 226 0 0 1 1 0 0 relative 0 +Rb yes Rb Rb 0 85.468 0 0 1 1 0 0 relative 0 +Re yes Re Re 0 186.207 0 0 1 1 0 0 relative 0 +Rh yes Rh Rh 0 102.906 0 0 1 1 0 0 relative 0 +Rn yes Rn Rn 0 222 0 0 1 1 0 0 relative 0 +Ru yes Ru Ru 0 101.07 0 0 1 1 0 0 relative 0 +S yes S S 0 32.066 0 0 1 1 0 0 relative 0 +Sb yes Sb Sb 0 121.76 0 0 1 1 0 0 relative 0 +Sc yes Sc Sc 0 44.956 0 0 1 1 0 0 relative 0 +Se yes Se Se 0 78.96 0 0 1 1 0 0 relative 0 +Si yes Si Si 0 28.086 0 0 1 1 0 0 relative 0 +Sm yes Sm Sm 0 150.36 0 0 1 1 0 0 relative 0 +Sn yes Sn Sn 0 118.71 0 0 1 1 0 0 relative 0 +Sr yes Sr Sr 0 87.62 0 0 1 1 0 0 relative 0 +Ta yes Ta Ta 0 180.948 0 0 1 1 0 0 relative 0 +Tb yes Tb Tb 0 158.925 0 0 1 1 0 0 relative 0 +Tc yes Tc Tc 0 98 0 0 1 1 0 0 relative 0 +Te yes Te Te 0 127.6 0 0 1 1 0 0 relative 0 +Th yes Th Th 0 232.038 0 0 1 1 0 0 relative 0 +Ti yes Ti Ti 0 47.867 0 0 1 1 0 0 relative 0 +Tl yes Tl Tl 0 204.383 0 0 1 1 0 0 relative 0 +Tm yes Tm Tm 0 168.934 0 0 1 1 0 0 relative 0 +U yes U U 0 238.029 0 0 1 1 0 0 relative 0 +V yes V V 0 50.942 0 0 1 1 0 0 relative 0 +W yes W W 0 183.84 0 0 1 1 0 0 relative 0 +Xe yes Xe Xe 0 131.29 0 0 1 1 0 0 relative 0 +Y yes Y Y 0 88.906 0 0 1 1 0 0 relative 0 +Yb yes Yb Yb 0 173.04 0 0 1 1 0 0 relative 0 +Zn yes Zn Zn 0 65.39 0 0 1 1 0 0 relative 0 +Zr yes Zr Zr 0 91.224 0 0 1 1 0 0 relative 0 +C_co2 yes C C 0 12.011 0.70 0 1 1 0 0 relative 0 +O_co2 yes O O 0 15.9994 -0.35 0 1 1 0 0 relative 0 +N_n2 yes N N 0 14.0067 -0.482 0 1 1 0 0 relative 0 +N_com no N - 0 0.0 0.964 0 1 1 0 0 relative 0 +Ow yes O O 0 15.9994 0.0 0 1.0 0.5 2 0 relative 0 +Hw yes H H 0 1.008 0.52 0 1.0 0.32 1 0 relative 0 +Lw no Lw - 0 0.0 -1.04 0 1.0 1.00 1 0 relative 0 +CH4_sp3 yes C C 0 16.04246 0.0 0 1.0 1.00 0 0 relative 0 +H_h2 yes H H 0 1.00794 0.468 0 1 1 0 0 relative 0 +H_com no H H 0 0.0 -0.936 0 1 1 0 0 relative 0 diff --git a/src/matkit/pygraspa/files/template/simulation.input b/src/matkit/pygraspa/files/template/simulation.input new file mode 100644 index 0000000..7aa77a9 --- /dev/null +++ b/src/matkit/pygraspa/files/template/simulation.input @@ -0,0 +1,30 @@ +NumberOfInitializationCycles NCYCLE +NumberOfEquilibrationCycles 0 +NumberOfProductionCycles NCYCLE +UseMaxStep no +MaxStepPerCycle 1 + +RestartFile no +BMCBiasingMethod LJ_Biasing +NumberOfTrialPositions 10 +NumberOfTrialOrientations 10 +NumberOfBlocks 1 +AdsorbateAllocateSpace 30240 + +NumberOfSimulations 1 +SingleSimulation yes +DifferentFrameworks yes + +UseChargesFromCIFFile yes +InputFileType cif +FrameworkName CIFFILE +UnitCells 0 UC_X UC_Y UC_Z +ChargeMethod Ewald +Temperature TEMPERATURE +Pressure PRESSURE +OverlapCriteria 1e5 +CutOffVDW CUTOFF +CutOffCoulomb CUTOFF +EwaldPrecision 1e-6 + +__COMPONENTS__ diff --git a/src/matkit/pygraspa/files/template_mixture/CO2.def b/src/matkit/pygraspa/files/template_mixture/CO2.def new file mode 100644 index 0000000..fa1ff0a --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/CO2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +304.1282 +7377300.0 +0.22394 +#Number Of Atoms + 3 +# Number of groups +1 +# CO2-group +rigid +# number of atoms +3 +# atomic positions +0 O_co2 0.0 0.0 1.16 +1 C_co2 0.0 0.0 0.0 +2 O_co2 0.0 0.0 -1.16 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture/N2.def b/src/matkit/pygraspa/files/template_mixture/N2.def new file mode 100644 index 0000000..53068c5 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/N2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +126.192 +3395800.0 +0.0372 +#Number Of Atoms +3 +# Number of groups +1 +# N2-group +rigid +# number of atoms +3 +# atomic positions +0 N_n2 0.0 0.0 0.55 +1 N_com 0.0 0.0 0.0 +2 N_n2 0.0 0.0 -0.55 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture/TIP4P.def b/src/matkit/pygraspa/files/template_mixture/TIP4P.def new file mode 100644 index 0000000..b6b4fd6 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/TIP4P.def @@ -0,0 +1,25 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +647.14 +22064000.0 +-0.217000 +# total number Of atoms +4 +# Number of groups +1 +# water-group +rigid +# number of atoms +4 +# atomic positions +0 Ow 0.0 0.0 0.0 +1 Lw 0.0 0.15 0.0 +2 Hw 0.75695 0.58588 0.0 +3 Hw -0.75695 0.58588 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +0 2 RIGID_BOND +0 3 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture/force_field.def b/src/matkit/pygraspa/files/template_mixture/force_field.def new file mode 100644 index 0000000..0e7bb78 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/force_field.def @@ -0,0 +1,6 @@ +# rules to overwrite +0 +# number of defined interactions +0 +# mixing rules to overwrite +0 diff --git a/src/matkit/pygraspa/files/template_mixture/force_field_mixing_rules.def b/src/matkit/pygraspa/files/template_mixture/force_field_mixing_rules.def new file mode 100644 index 0000000..ce34ef4 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/force_field_mixing_rules.def @@ -0,0 +1,119 @@ +# general rule for shifted vs truncated +truncated +# general rule tailcorrections +no +# number of defined interactions +110 +# type interaction +Ac lennard-jones 16.60 3.10 +Ag lennard-jones 18.11 2.80 +Al lennard-jones 254.09 4.01 +Am lennard-jones 7.04 3.01 +Ar lennard-jones 93.08 3.45 +As lennard-jones 155.47 3.77 +At lennard-jones 142.89 4.23 +Au lennard-jones 19.62 2.93 +B lennard-jones 90.57 3.64 +Ba lennard-jones 183.15 3.30 +Be lennard-jones 42.77 2.45 +Bi lennard-jones 260.63 3.89 +Bk lennard-jones 6.54 2.97 +Br lennard-jones 126.30 3.73 +C lennard-jones 52.83 3.43 +Ca lennard-jones 119.75 3.03 +Cd lennard-jones 114.72 2.54 +Ce lennard-jones 6.54 3.17 +Cf lennard-jones 6.54 2.95 +Cl lennard-jones 114.23 3.517 +Cm lennard-jones 6.54 2.96 +Co lennard-jones 7.04 2.56 +Cr lennard-jones 7.55 2.69 +Cs lennard-jones 22.64 4.02 +Cu lennard-jones 2.516 3.114 +Dy lennard-jones 3.52 3.05 +Er lennard-jones 3.52 3.02 +Es lennard-jones 6.04 2.94 +Eu lennard-jones 4.03 3.11 +F lennard-jones 25.16 2.997 +Fe lennard-jones 6.54 2.59 +Fm lennard-jones 6.04 2.93 +Fr lennard-jones 25.16 4.37 +Ga lennard-jones 208.81 3.90 +Gd lennard-jones 4.53 3.00 +Ge lennard-jones 190.69 3.81 +H lennard-jones 22.14 2.57 +He lennard-jones 10.9 2.64 +Hf lennard-jones 36.23 2.80 +Hg lennard-jones 193.71 2.41 +Ho lennard-jones 3.52 3.04 +I lennard-jones 170.57 4.01 +In lennard-jones 301.39 3.98 +Ir lennard-jones 36.73 2.53 +K lennard-jones 17.61 3.40 +Kr lennard-jones 166.40 3.636 +La lennard-jones 8.55 3.14 +Li lennard-jones 12.58 2.18 +Lu lennard-jones 20.63 3.24 +Md lennard-jones 5.53 2.92 +Mg lennard-jones 55.85 2.69 +Mn lennard-jones 6.54 2.64 +Mo lennard-jones 28.18 2.72 +N lennard-jones 34.72 3.26 +Na lennard-jones 15.09 2.66 +Nb lennard-jones 29.69 2.82 +Nd lennard-jones 5.03 3.18 +Ne lennard-jones 21.13 2.66 +Ni lennard-jones 7.55 2.52 +No lennard-jones 5.53 2.89 +Np lennard-jones 9.56 3.05 +O lennard-jones 30.19 3.12 +Os lennard-jones 18.62 2.78 +P lennard-jones 153.46 3.69 +Pa lennard-jones 11.07 3.05 +Pb lennard-jones 333.59 3.83 +Pd lennard-jones 24.15 2.58 +Pm lennard-jones 4.53 3.16 +Po lennard-jones 163.52 4.20 +Pr lennard-jones 5.03 3.21 +Pt lennard-jones 40.25 2.45 +Pu lennard-jones 8.05 3.05 +Ra lennard-jones 203.27 3.28 +Rb lennard-jones 20.13 3.67 +Re lennard-jones 33.21 2.63 +Rh lennard-jones 26.67 2.61 +Rn lennard-jones 124.78 4.25 +Ru lennard-jones 28.18 2.64 +S lennard-jones 137.86 3.59 +Sb lennard-jones 225.91 3.94 +Sc lennard-jones 9.56 2.94 +Se lennard-jones 146.42 3.75 +Si lennard-jones 202.27 3.83 +Sm lennard-jones 4.03 3.14 +Sn lennard-jones 285.28 3.91 +Sr lennard-jones 118.24 3.24 +Ta lennard-jones 40.75 2.82 +Tb lennard-jones 3.52 3.07 +Tc lennard-jones 24.15 2.67 +Te lennard-jones 200.25 3.98 +Th lennard-jones 13.08 3.03 +Ti lennard-jones 8.55 2.83 +Tl lennard-jones 342.14 3.87 +Tm lennard-jones 3.02 3.01 +U lennard-jones 11.07 3.02 +V lennard-jones 8.05 2.80 +W lennard-jones 33.71 2.73 +Xe lennard-jones 221.00 4.10 +Y lennard-jones 36.23 2.98 +Yb lennard-jones 114.72 2.99 +Zn lennard-jones 62.40 2.46 +Zr lennard-jones 34.72 2.783 +C_co2 lennard-jones 27.0 2.80 +O_co2 lennard-jones 79.0 3.05 +N_n2 lennard-jones 36.0 3.31 +N_com lennard-jones 0.0 0.0 +Ow lennard-jones 78.0 3.154 +Hw lennard-jones 0.0 0.0 +Lw lennard-jones 0.0 0.0 +CH4_sp3 lennard-jones 148.0 3.73 +# general mixing rule for Lennard-Jones +Lorentz-Berthelot diff --git a/src/matkit/pygraspa/files/template_mixture/methane.def b/src/matkit/pygraspa/files/template_mixture/methane.def new file mode 100644 index 0000000..5b5aec7 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/methane.def @@ -0,0 +1,18 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +190.564 +4599200.0 +0.01142 +# Number Of Atoms +1 +# Number Of Groups +1 +# Alkane-group +rigid +# number of atoms +1 +# atomic positions +0 CH4_sp3 0.0 0.0 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Bond/Bend Bend/Bend Bond/Torsion Bend/Torsion IntraVDW Intra ch-ch Intra ch-bd Intra bd-bd + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Number of config moves +0 \ No newline at end of file diff --git a/src/matkit/pygraspa/files/template_mixture/pseudo_atoms.def b/src/matkit/pygraspa/files/template_mixture/pseudo_atoms.def new file mode 100644 index 0000000..8d728d1 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/pseudo_atoms.def @@ -0,0 +1,113 @@ +#number of pseudo atoms +110 +#type print as chem oxi mass charge polarization B-factor radii connectivity anisotropic anisotropic-type tinker-type +Ac yes Ac Ac 0 227 0 0 1 1 0 0 relative 0 +Ag yes Ag Ag 0 107.868 0 0 1 1 0 0 relative 0 +Al yes Al Al 0 26.982 0 0 1 1 0 0 relative 0 +Am yes Am Am 0 243 0 0 1 1 0 0 relative 0 +Ar yes Ar Ar 0 39.948 0 0 1 1 0 0 relative 0 +As yes As As 0 74.922 0 0 1 1 0 0 relative 0 +At yes At At 0 210 0 0 1 1 0 0 relative 0 +Au yes Au Au 0 196.967 0 0 1 1 0 0 relative 0 +B yes B B 0 10.811 0 0 1 1 0 0 relative 0 +Ba yes Ba Ba 0 137.327 0 0 1 1 0 0 relative 0 +Be yes Be Be 0 9.012 0 0 1 1 0 0 relative 0 +Bi yes Bi Bi 0 208.98 0 0 1 1 0 0 relative 0 +Bk yes Bk Bk 0 247 0 0 1 1 0 0 relative 0 +Br yes Br Br 0 79.904 0 0 1 1 0 0 relative 0 +C yes C C 0 12.011 0 0 1 1 0 0 relative 0 +Ca yes Ca Ca 0 40.078 0 0 1 1 0 0 relative 0 +Cd yes Cd Cd 0 112.411 0 0 1 1 0 0 relative 0 +Ce yes Ce Ce 0 140.116 0 0 1 1 0 0 relative 0 +Cf yes Cf Cf 0 251 0 0 1 1 0 0 relative 0 +Cl yes Cl Cl 0 35.453 0 0 1 1 0 0 relative 0 +Cm yes Cm Cm 0 247 0 0 1 1 0 0 relative 0 +Co yes Co Co 0 58.933 0 0 1 1 0 0 relative 0 +Cr yes Cr Cr 0 51.996 0 0 1 1 0 0 relative 0 +Cs yes Cs Cs 0 132.905 0 0 1 1 0 0 relative 0 +Cu yes Cu Cu 0 63.546 0 0 1 1 0 0 relative 0 +Dy yes Dy Dy 0 162.5 0 0 1 1 0 0 relative 0 +Er yes Er Er 0 167.26 0 0 1 1 0 0 relative 0 +Es yes Es Es 0 252 0 0 1 1 0 0 relative 0 +Eu yes Eu Eu 0 151.964 0 0 1 1 0 0 relative 0 +F yes F F 0 18.998 0 0 1 1 0 0 relative 0 +Fe yes Fe Fe 0 55.845 0 0 1 1 0 0 relative 0 +Fm yes Fm Fm 0 257 0 0 1 1 0 0 relative 0 +Fr yes Fr Fr 0 223 0 0 1 1 0 0 relative 0 +Ga yes Ga Ga 0 69.723 0 0 1 1 0 0 relative 0 +Gd yes Gd Gd 0 157.25 0 0 1 1 0 0 relative 0 +Ge yes Ge Ge 0 72.61 0 0 1 1 0 0 relative 0 +H yes H H 0 1.008 0 0 1 1 0 0 relative 0 +He yes He He 0 4.003 0 0 1 1 0 0 relative 0 +Hf yes Hf Hf 0 178.49 0 0 1 1 0 0 relative 0 +Hg yes Hg Hg 0 200.59 0 0 1 1 0 0 relative 0 +Ho yes Ho Ho 0 164.93 0 0 1 1 0 0 relative 0 +I yes I I 0 126.904 0 0 1 1 0 0 relative 0 +In yes In In 0 114.818 0 0 1 1 0 0 relative 0 +Ir yes Ir Ir 0 192.217 0 0 1 1 0 0 relative 0 +K yes K K 0 39.098 0 0 1 1 0 0 relative 0 +Kr yes Kr Kr 0 83.8 0 0 1 1 0 0 relative 0 +La yes La La 0 138.906 0 0 1 1 0 0 relative 0 +Li yes Li Li 0 6.941 0 0 1 1 0 0 relative 0 +Lu yes Lu Lu 0 174.967 0 0 1 1 0 0 relative 0 +Md yes Md Md 0 258 0 0 1 1 0 0 relative 0 +Mg yes Mg Mg 0 24.305 0 0 1 1 0 0 relative 0 +Mn yes Mn Mn 0 54.938 0 0 1 1 0 0 relative 0 +Mo yes Mo Mo 0 95.94 0 0 1 1 0 0 relative 0 +N yes N N 0 14.007 0 0 1 1 0 0 relative 0 +Na yes Na Na 0 22.991 0 0 1 1 0 0 relative 0 +Nb yes Nb Nb 0 92.906 0 0 1 1 0 0 relative 0 +Nd yes Nd Nd 0 144.24 0 0 1 1 0 0 relative 0 +Ne yes Ne Ne 0 20.18 0 0 1 1 0 0 relative 0 +Ni yes Ni Ni 0 58.693 0 0 1 1 0 0 relative 0 +No yes No No 0 259 0 0 1 1 0 0 relative 0 +Np yes Np Np 0 237 0 0 1 1 0 0 relative 0 +O yes O O 0 15.999 0 0 1 1 0 0 relative 0 +Os yes Os Os 0 190.23 0 0 1 1 0 0 relative 0 +P yes P P 0 30.974 0 0 1 1 0 0 relative 0 +Pa yes Pa Pa 0 231.036 0 0 1 1 0 0 relative 0 +Pb yes Pb Pb 0 207.2 0 0 1 1 0 0 relative 0 +Pd yes Pd Pd 0 106.42 0 0 1 1 0 0 relative 0 +Pm yes Pm Pm 0 145 0 0 1 1 0 0 relative 0 +Po yes Po Po 0 210 0 0 1 1 0 0 relative 0 +Pr yes Pr Pr 0 140.908 0 0 1 1 0 0 relative 0 +Pt yes Pt Pt 0 195.078 0 0 1 1 0 0 relative 0 +Pu yes Pu Pu 0 244 0 0 1 1 0 0 relative 0 +Ra yes Ra Ra 0 226 0 0 1 1 0 0 relative 0 +Rb yes Rb Rb 0 85.468 0 0 1 1 0 0 relative 0 +Re yes Re Re 0 186.207 0 0 1 1 0 0 relative 0 +Rh yes Rh Rh 0 102.906 0 0 1 1 0 0 relative 0 +Rn yes Rn Rn 0 222 0 0 1 1 0 0 relative 0 +Ru yes Ru Ru 0 101.07 0 0 1 1 0 0 relative 0 +S yes S S 0 32.066 0 0 1 1 0 0 relative 0 +Sb yes Sb Sb 0 121.76 0 0 1 1 0 0 relative 0 +Sc yes Sc Sc 0 44.956 0 0 1 1 0 0 relative 0 +Se yes Se Se 0 78.96 0 0 1 1 0 0 relative 0 +Si yes Si Si 0 28.086 0 0 1 1 0 0 relative 0 +Sm yes Sm Sm 0 150.36 0 0 1 1 0 0 relative 0 +Sn yes Sn Sn 0 118.71 0 0 1 1 0 0 relative 0 +Sr yes Sr Sr 0 87.62 0 0 1 1 0 0 relative 0 +Ta yes Ta Ta 0 180.948 0 0 1 1 0 0 relative 0 +Tb yes Tb Tb 0 158.925 0 0 1 1 0 0 relative 0 +Tc yes Tc Tc 0 98 0 0 1 1 0 0 relative 0 +Te yes Te Te 0 127.6 0 0 1 1 0 0 relative 0 +Th yes Th Th 0 232.038 0 0 1 1 0 0 relative 0 +Ti yes Ti Ti 0 47.867 0 0 1 1 0 0 relative 0 +Tl yes Tl Tl 0 204.383 0 0 1 1 0 0 relative 0 +Tm yes Tm Tm 0 168.934 0 0 1 1 0 0 relative 0 +U yes U U 0 238.029 0 0 1 1 0 0 relative 0 +V yes V V 0 50.942 0 0 1 1 0 0 relative 0 +W yes W W 0 183.84 0 0 1 1 0 0 relative 0 +Xe yes Xe Xe 0 131.29 0 0 1 1 0 0 relative 0 +Y yes Y Y 0 88.906 0 0 1 1 0 0 relative 0 +Yb yes Yb Yb 0 173.04 0 0 1 1 0 0 relative 0 +Zn yes Zn Zn 0 65.39 0 0 1 1 0 0 relative 0 +Zr yes Zr Zr 0 91.224 0 0 1 1 0 0 relative 0 +C_co2 yes C C 0 12.011 0.70 0 1 1 0 0 relative 0 +O_co2 yes O O 0 15.9994 -0.35 0 1 1 0 0 relative 0 +N_n2 yes N N 0 14.0067 -0.482 0 1 1 0 0 relative 0 +N_com no N - 0 0.0 0.964 0 1 1 0 0 relative 0 +Ow yes O O 0 15.9994 0.0 0 1.0 0.5 2 0 relative 0 +Hw yes H H 0 1.008 0.52 0 1.0 0.32 1 0 relative 0 +Lw no Lw - 0 0.0 -1.04 0 1.0 1.00 1 0 relative 0 +CH4_sp3 yes C C 0 16.04246 0.0 0 1.0 1.00 0 0 relative 0 diff --git a/src/matkit/pygraspa/files/template_mixture/simulation.input b/src/matkit/pygraspa/files/template_mixture/simulation.input new file mode 100644 index 0000000..e4e5067 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture/simulation.input @@ -0,0 +1,57 @@ +NumberOfInitializationCycles NCYCLE +NumberOfEquilibrationCycles 0 +NumberOfProductionCycles NCYCLE + +RestartFile no +BMCBiasingMethod LJ_Biasing +NumberOfTrialPositions 10 +NumberOfTrialOrientations 10 +NumberOfBlocks 1 +AdsorbateAllocateSpace 30240 + +NumberOfSimulations 1 +SingleSimulation yes +DifferentFrameworks yes + +InputFileType cif +FrameworkName CIFFILE +UnitCells 0 UC_X UC_Y UC_Z +ChargeMethod Ewald +Temperature TEMPERATURE +Pressure PRESSURE +OverlapCriteria 1e5 +CutOffVDW CUTOFF +CutOffCoulomb CUTOFF +EwaldPrecision 1e-6 +UseChargesFromCIFFile yes + +Component 0 MoleculeName CO2 + IdealGasRosenbluthWeight 1.0 + FugacityCoefficient PR-EOS + MolFraction 0.14 + TranslationProbability 1.0 + RotationProbability 1.0 + ReinsertionProbability 1.0 + IdentityChangeProbability 1.0 + SwapProbability 2.0 + CreateNumberOfMolecules 0 +Component 1 MoleculeName N2 + IdealGasRosenbluthWeight 1.0 + FugacityCoefficient PR-EOS + MolFraction 0.8553 + TranslationProbability 1.0 + RotationProbability 1.0 + ReinsertionProbability 1.0 + IdentityChangeProbability 1.0 + SwapProbability 2.0 + CreateNumberOfMolecules 0 +Component 2 MoleculeName TIP4P + IdealGasRosenbluthWeight 1.0 + FugacityCoefficient PR-EOS + MolFraction 0.0047 + TranslationProbability 1.0 + RotationProbability 1.0 + ReinsertionProbability 1.0 + IdentityChangeProbability 1.0 + SwapProbability 2.0 + CreateNumberOfMolecules 0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/CO2.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/CO2.def new file mode 100644 index 0000000..fa1ff0a --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/CO2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +304.1282 +7377300.0 +0.22394 +#Number Of Atoms + 3 +# Number of groups +1 +# CO2-group +rigid +# number of atoms +3 +# atomic positions +0 O_co2 0.0 0.0 1.16 +1 C_co2 0.0 0.0 0.0 +2 O_co2 0.0 0.0 -1.16 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/H2.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/H2.def new file mode 100644 index 0000000..89598c2 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/H2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +33.19 +1315000.0 +-0.214 +#Number Of Atoms +3 +# Number of groups +1 +# H2-group +rigid +# number of atoms +3 +# atomic positions +0 H_h2 0.0 0.0 0.37 +1 H_com 0.0 0.0 0.0 +2 H_h2 0.0 0.0 -0.37 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/N2.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/N2.def new file mode 100644 index 0000000..53068c5 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/N2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +126.192 +3395800.0 +0.0372 +#Number Of Atoms +3 +# Number of groups +1 +# N2-group +rigid +# number of atoms +3 +# atomic positions +0 N_n2 0.0 0.0 0.55 +1 N_com 0.0 0.0 0.0 +2 N_n2 0.0 0.0 -0.55 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/SO2.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/SO2.def new file mode 100644 index 0000000..f9b97bf --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/SO2.def @@ -0,0 +1,23 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +430.8 +7883100.0 +0.256 +#Number Of Atoms + 3 +# Number of groups +1 +# SO2-group +rigid +# number of atoms +3 +# atomic positions +0 O_so2 0.0 0.0 1.432 +1 S_so2 0.0 0.0 0.0 +2 O_so2 0.0 0.0 -1.432 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +1 2 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/TIP4P.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/TIP4P.def new file mode 100644 index 0000000..b6b4fd6 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/TIP4P.def @@ -0,0 +1,25 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +647.14 +22064000.0 +-0.217000 +# total number Of atoms +4 +# Number of groups +1 +# water-group +rigid +# number of atoms +4 +# atomic positions +0 Ow 0.0 0.0 0.0 +1 Lw 0.0 0.15 0.0 +2 Hw 0.75695 0.58588 0.0 +3 Hw -0.75695 0.58588 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Stretch/Bend Bend/Bend Stretch/Torsion Bend/Torsion IntraVDW IntraCoulomb + 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Bond stretch: atom n1-n2, type, parameters +0 1 RIGID_BOND +0 2 RIGID_BOND +0 3 RIGID_BOND +# Number of config moves +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field.def new file mode 100644 index 0000000..0e7bb78 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field.def @@ -0,0 +1,6 @@ +# rules to overwrite +0 +# number of defined interactions +0 +# mixing rules to overwrite +0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field_mixing_rules.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field_mixing_rules.def new file mode 100644 index 0000000..16df1b9 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/force_field_mixing_rules.def @@ -0,0 +1,123 @@ +# general rule for shifted vs truncated +truncated +# general rule tailcorrections +no +# number of defined interactions +114 +# type interaction +Ac lennard-jones 16.60 3.10 +Ag lennard-jones 18.11 2.80 +Al lennard-jones 254.09 4.01 +Am lennard-jones 7.04 3.01 +Ar lennard-jones 93.08 3.45 +As lennard-jones 155.47 3.77 +At lennard-jones 142.89 4.23 +Au lennard-jones 19.62 2.93 +B lennard-jones 90.57 3.64 +Ba lennard-jones 183.15 3.30 +Be lennard-jones 42.77 2.45 +Bi lennard-jones 260.63 3.89 +Bk lennard-jones 6.54 2.97 +Br lennard-jones 126.30 3.73 +C lennard-jones 52.83 3.43 +Ca lennard-jones 119.75 3.03 +Cd lennard-jones 114.72 2.54 +Ce lennard-jones 6.54 3.17 +Cf lennard-jones 6.54 2.95 +Cl lennard-jones 114.23 3.517 +Cm lennard-jones 6.54 2.96 +Co lennard-jones 7.04 2.56 +Cr lennard-jones 7.55 2.69 +Cs lennard-jones 22.64 4.02 +Cu lennard-jones 2.516 3.114 +Dy lennard-jones 3.52 3.05 +Er lennard-jones 3.52 3.02 +Es lennard-jones 6.04 2.94 +Eu lennard-jones 4.03 3.11 +F lennard-jones 25.16 2.997 +Fe lennard-jones 6.54 2.59 +Fm lennard-jones 6.04 2.93 +Fr lennard-jones 25.16 4.37 +Ga lennard-jones 208.81 3.90 +Gd lennard-jones 4.53 3.00 +Ge lennard-jones 190.69 3.81 +H lennard-jones 22.14 2.57 +He lennard-jones 10.9 2.64 +Hf lennard-jones 36.23 2.80 +Hg lennard-jones 193.71 2.41 +Ho lennard-jones 3.52 3.04 +I lennard-jones 170.57 4.01 +In lennard-jones 301.39 3.98 +Ir lennard-jones 36.73 2.53 +K lennard-jones 17.61 3.40 +Kr lennard-jones 166.40 3.636 +La lennard-jones 8.55 3.14 +Li lennard-jones 12.58 2.18 +Lu lennard-jones 20.63 3.24 +Md lennard-jones 5.53 2.92 +Mg lennard-jones 55.85 2.69 +Mn lennard-jones 6.54 2.64 +Mo lennard-jones 28.18 2.72 +N lennard-jones 34.72 3.26 +Na lennard-jones 15.09 2.66 +Nb lennard-jones 29.69 2.82 +Nd lennard-jones 5.03 3.18 +Ne lennard-jones 21.13 2.66 +Ni lennard-jones 7.55 2.52 +No lennard-jones 5.53 2.89 +Np lennard-jones 9.56 3.05 +O lennard-jones 30.19 3.12 +Os lennard-jones 18.62 2.78 +P lennard-jones 153.46 3.69 +Pa lennard-jones 11.07 3.05 +Pb lennard-jones 333.59 3.83 +Pd lennard-jones 24.15 2.58 +Pm lennard-jones 4.53 3.16 +Po lennard-jones 163.52 4.20 +Pr lennard-jones 5.03 3.21 +Pt lennard-jones 40.25 2.45 +Pu lennard-jones 8.05 3.05 +Ra lennard-jones 203.27 3.28 +Rb lennard-jones 20.13 3.67 +Re lennard-jones 33.21 2.63 +Rh lennard-jones 26.67 2.61 +Rn lennard-jones 124.78 4.25 +Ru lennard-jones 28.18 2.64 +S lennard-jones 137.86 3.59 +Sb lennard-jones 225.91 3.94 +Sc lennard-jones 9.56 2.94 +Se lennard-jones 146.42 3.75 +Si lennard-jones 202.27 3.83 +Sm lennard-jones 4.03 3.14 +Sn lennard-jones 285.28 3.91 +Sr lennard-jones 118.24 3.24 +Ta lennard-jones 40.75 2.82 +Tb lennard-jones 3.52 3.07 +Tc lennard-jones 24.15 2.67 +Te lennard-jones 200.25 3.98 +Th lennard-jones 13.08 3.03 +Ti lennard-jones 8.55 2.83 +Tl lennard-jones 342.14 3.87 +Tm lennard-jones 3.02 3.01 +U lennard-jones 11.07 3.02 +V lennard-jones 8.05 2.80 +W lennard-jones 33.71 2.73 +Xe lennard-jones 221.00 4.10 +Y lennard-jones 36.23 2.98 +Yb lennard-jones 114.72 2.99 +Zn lennard-jones 62.40 2.46 +Zr lennard-jones 34.72 2.783 +C_co2 lennard-jones 27.0 2.80 +O_co2 lennard-jones 79.0 3.05 +N_n2 lennard-jones 36.0 3.31 +N_com lennard-jones 0.0 0.0 +Ow lennard-jones 78.0 3.154 +Hw lennard-jones 0.0 0.0 +Lw lennard-jones 0.0 0.0 +CH4_sp3 lennard-jones 148.0 3.73 +S_so2 lennard-jones 73.8 3.39 +O_so2 lennard-jones 79.0 3.05 +H_h2 lennard-jones 0.0 0.0 +H_com lennard-jones 36.7 2.958 +# general mixing rule for Lennard-Jones +Lorentz-Berthelot diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/methane.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/methane.def new file mode 100644 index 0000000..5b5aec7 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/methane.def @@ -0,0 +1,18 @@ +# critical constants: Temperature [T], Pressure [Pa], and Acentric factor [-] +190.564 +4599200.0 +0.01142 +# Number Of Atoms +1 +# Number Of Groups +1 +# Alkane-group +rigid +# number of atoms +1 +# atomic positions +0 CH4_sp3 0.0 0.0 0.0 +# Chiral centers Bond BondDipoles Bend UrayBradley InvBend Torsion Imp. Torsion Bond/Bond Bond/Bend Bend/Bend Bond/Torsion Bend/Torsion IntraVDW Intra ch-ch Intra ch-bd Intra bd-bd + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +# Number of config moves +0 \ No newline at end of file diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/pseudo_atoms.def b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/pseudo_atoms.def new file mode 100644 index 0000000..8ee2cd7 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/pseudo_atoms.def @@ -0,0 +1,117 @@ +#number of pseudo atoms +114 +#type print as chem oxi mass charge polarization B-factor radii connectivity anisotropic anisotropic-type tinker-type +Ac yes Ac Ac 0 227 0 0 1 1 0 0 relative 0 +Ag yes Ag Ag 0 107.868 0 0 1 1 0 0 relative 0 +Al yes Al Al 0 26.982 0 0 1 1 0 0 relative 0 +Am yes Am Am 0 243 0 0 1 1 0 0 relative 0 +Ar yes Ar Ar 0 39.948 0 0 1 1 0 0 relative 0 +As yes As As 0 74.922 0 0 1 1 0 0 relative 0 +At yes At At 0 210 0 0 1 1 0 0 relative 0 +Au yes Au Au 0 196.967 0 0 1 1 0 0 relative 0 +B yes B B 0 10.811 0 0 1 1 0 0 relative 0 +Ba yes Ba Ba 0 137.327 0 0 1 1 0 0 relative 0 +Be yes Be Be 0 9.012 0 0 1 1 0 0 relative 0 +Bi yes Bi Bi 0 208.98 0 0 1 1 0 0 relative 0 +Bk yes Bk Bk 0 247 0 0 1 1 0 0 relative 0 +Br yes Br Br 0 79.904 0 0 1 1 0 0 relative 0 +C yes C C 0 12.011 0 0 1 1 0 0 relative 0 +Ca yes Ca Ca 0 40.078 0 0 1 1 0 0 relative 0 +Cd yes Cd Cd 0 112.411 0 0 1 1 0 0 relative 0 +Ce yes Ce Ce 0 140.116 0 0 1 1 0 0 relative 0 +Cf yes Cf Cf 0 251 0 0 1 1 0 0 relative 0 +Cl yes Cl Cl 0 35.453 0 0 1 1 0 0 relative 0 +Cm yes Cm Cm 0 247 0 0 1 1 0 0 relative 0 +Co yes Co Co 0 58.933 0 0 1 1 0 0 relative 0 +Cr yes Cr Cr 0 51.996 0 0 1 1 0 0 relative 0 +Cs yes Cs Cs 0 132.905 0 0 1 1 0 0 relative 0 +Cu yes Cu Cu 0 63.546 0 0 1 1 0 0 relative 0 +Dy yes Dy Dy 0 162.5 0 0 1 1 0 0 relative 0 +Er yes Er Er 0 167.26 0 0 1 1 0 0 relative 0 +Es yes Es Es 0 252 0 0 1 1 0 0 relative 0 +Eu yes Eu Eu 0 151.964 0 0 1 1 0 0 relative 0 +F yes F F 0 18.998 0 0 1 1 0 0 relative 0 +Fe yes Fe Fe 0 55.845 0 0 1 1 0 0 relative 0 +Fm yes Fm Fm 0 257 0 0 1 1 0 0 relative 0 +Fr yes Fr Fr 0 223 0 0 1 1 0 0 relative 0 +Ga yes Ga Ga 0 69.723 0 0 1 1 0 0 relative 0 +Gd yes Gd Gd 0 157.25 0 0 1 1 0 0 relative 0 +Ge yes Ge Ge 0 72.61 0 0 1 1 0 0 relative 0 +H yes H H 0 1.008 0 0 1 1 0 0 relative 0 +He yes He He 0 4.003 0 0 1 1 0 0 relative 0 +Hf yes Hf Hf 0 178.49 0 0 1 1 0 0 relative 0 +Hg yes Hg Hg 0 200.59 0 0 1 1 0 0 relative 0 +Ho yes Ho Ho 0 164.93 0 0 1 1 0 0 relative 0 +I yes I I 0 126.904 0 0 1 1 0 0 relative 0 +In yes In In 0 114.818 0 0 1 1 0 0 relative 0 +Ir yes Ir Ir 0 192.217 0 0 1 1 0 0 relative 0 +K yes K K 0 39.098 0 0 1 1 0 0 relative 0 +Kr yes Kr Kr 0 83.8 0 0 1 1 0 0 relative 0 +La yes La La 0 138.906 0 0 1 1 0 0 relative 0 +Li yes Li Li 0 6.941 0 0 1 1 0 0 relative 0 +Lu yes Lu Lu 0 174.967 0 0 1 1 0 0 relative 0 +Md yes Md Md 0 258 0 0 1 1 0 0 relative 0 +Mg yes Mg Mg 0 24.305 0 0 1 1 0 0 relative 0 +Mn yes Mn Mn 0 54.938 0 0 1 1 0 0 relative 0 +Mo yes Mo Mo 0 95.94 0 0 1 1 0 0 relative 0 +N yes N N 0 14.007 0 0 1 1 0 0 relative 0 +Na yes Na Na 0 22.991 0 0 1 1 0 0 relative 0 +Nb yes Nb Nb 0 92.906 0 0 1 1 0 0 relative 0 +Nd yes Nd Nd 0 144.24 0 0 1 1 0 0 relative 0 +Ne yes Ne Ne 0 20.18 0 0 1 1 0 0 relative 0 +Ni yes Ni Ni 0 58.693 0 0 1 1 0 0 relative 0 +No yes No No 0 259 0 0 1 1 0 0 relative 0 +Np yes Np Np 0 237 0 0 1 1 0 0 relative 0 +O yes O O 0 15.999 0 0 1 1 0 0 relative 0 +Os yes Os Os 0 190.23 0 0 1 1 0 0 relative 0 +P yes P P 0 30.974 0 0 1 1 0 0 relative 0 +Pa yes Pa Pa 0 231.036 0 0 1 1 0 0 relative 0 +Pb yes Pb Pb 0 207.2 0 0 1 1 0 0 relative 0 +Pd yes Pd Pd 0 106.42 0 0 1 1 0 0 relative 0 +Pm yes Pm Pm 0 145 0 0 1 1 0 0 relative 0 +Po yes Po Po 0 210 0 0 1 1 0 0 relative 0 +Pr yes Pr Pr 0 140.908 0 0 1 1 0 0 relative 0 +Pt yes Pt Pt 0 195.078 0 0 1 1 0 0 relative 0 +Pu yes Pu Pu 0 244 0 0 1 1 0 0 relative 0 +Ra yes Ra Ra 0 226 0 0 1 1 0 0 relative 0 +Rb yes Rb Rb 0 85.468 0 0 1 1 0 0 relative 0 +Re yes Re Re 0 186.207 0 0 1 1 0 0 relative 0 +Rh yes Rh Rh 0 102.906 0 0 1 1 0 0 relative 0 +Rn yes Rn Rn 0 222 0 0 1 1 0 0 relative 0 +Ru yes Ru Ru 0 101.07 0 0 1 1 0 0 relative 0 +S yes S S 0 32.066 0 0 1 1 0 0 relative 0 +Sb yes Sb Sb 0 121.76 0 0 1 1 0 0 relative 0 +Sc yes Sc Sc 0 44.956 0 0 1 1 0 0 relative 0 +Se yes Se Se 0 78.96 0 0 1 1 0 0 relative 0 +Si yes Si Si 0 28.086 0 0 1 1 0 0 relative 0 +Sm yes Sm Sm 0 150.36 0 0 1 1 0 0 relative 0 +Sn yes Sn Sn 0 118.71 0 0 1 1 0 0 relative 0 +Sr yes Sr Sr 0 87.62 0 0 1 1 0 0 relative 0 +Ta yes Ta Ta 0 180.948 0 0 1 1 0 0 relative 0 +Tb yes Tb Tb 0 158.925 0 0 1 1 0 0 relative 0 +Tc yes Tc Tc 0 98 0 0 1 1 0 0 relative 0 +Te yes Te Te 0 127.6 0 0 1 1 0 0 relative 0 +Th yes Th Th 0 232.038 0 0 1 1 0 0 relative 0 +Ti yes Ti Ti 0 47.867 0 0 1 1 0 0 relative 0 +Tl yes Tl Tl 0 204.383 0 0 1 1 0 0 relative 0 +Tm yes Tm Tm 0 168.934 0 0 1 1 0 0 relative 0 +U yes U U 0 238.029 0 0 1 1 0 0 relative 0 +V yes V V 0 50.942 0 0 1 1 0 0 relative 0 +W yes W W 0 183.84 0 0 1 1 0 0 relative 0 +Xe yes Xe Xe 0 131.29 0 0 1 1 0 0 relative 0 +Y yes Y Y 0 88.906 0 0 1 1 0 0 relative 0 +Yb yes Yb Yb 0 173.04 0 0 1 1 0 0 relative 0 +Zn yes Zn Zn 0 65.39 0 0 1 1 0 0 relative 0 +Zr yes Zr Zr 0 91.224 0 0 1 1 0 0 relative 0 +C_co2 yes C C 0 12.011 0.70 0 1 1 0 0 relative 0 +O_co2 yes O O 0 15.9994 -0.35 0 1 1 0 0 relative 0 +N_n2 yes N N 0 14.0067 -0.482 0 1 1 0 0 relative 0 +N_com no N - 0 0.0 0.964 0 1 1 0 0 relative 0 +Ow yes O O 0 15.9994 0.0 0 1.0 0.5 2 0 relative 0 +Hw yes H H 0 1.008 0.52 0 1.0 0.32 1 0 relative 0 +Lw no Lw - 0 0.0 -1.04 0 1.0 1.00 1 0 relative 0 +CH4_sp3 yes C C 0 16.04246 0.0 0 1.0 1.00 0 0 relative 0 +S_so2 yes S S 0 32.066 0.59 0 1 1 0 0 relative 0 +O_so2 yes O O 0 15.9994 -0.295 0 1 1 0 0 relative 0 +H_h2 yes H H 0 1.00794 0.468 0 1 1 0 0 relative 0 +H_com no H H 0 0.0 -0.936 0 1 1 0 0 relative 0 diff --git a/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/simulation.input b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/simulation.input new file mode 100644 index 0000000..b35ca70 --- /dev/null +++ b/src/matkit/pygraspa/files/template_mixture_isotherm/template_mixture/simulation.input @@ -0,0 +1,29 @@ +NumberOfInitializationCycles NCYCLE +NumberOfEquilibrationCycles 0 +NumberOfProductionCycles NCYCLE + +RestartFile no +BMCBiasingMethod LJ_Biasing +NumberOfTrialPositions 10 +NumberOfTrialOrientations 10 +NumberOfBlocks 1 +AdsorbateAllocateSpace 30240 + +NumberOfSimulations 1 +SingleSimulation yes +DifferentFrameworks yes + +InputFileType cif +FrameworkName CIFFILE +UnitCells 0 UC_X UC_Y UC_Z +ChargeMethod Ewald +Temperature TEMPERATURE +Pressure PRESSURE +OverlapCriteria 1e5 +CutOffVDW CUTOFF +CutOffCoulomb CUTOFF +EwaldPrecision 1e-6 +UseChargesFromCIFFile yes + +__COMPONENTS__ + diff --git a/src/matkit/pygraspa/pygraspa.py b/src/matkit/pygraspa/pygraspa.py new file mode 100644 index 0000000..37ed162 --- /dev/null +++ b/src/matkit/pygraspa/pygraspa.py @@ -0,0 +1,661 @@ +"""matkit backend for [pygRASPA](https://github.com/tdpham2/pygRASPA). + +pygRASPA runs gRASPA Monte Carlo adsorption simulations with ML +interatomic potentials (FAIRChem eSEN/UMA/AllScAIP, MACE-Polar, Allegro). +This module produces simulation directories ready to launch with +``pygraspa run`` / ``pygraspa run-auto`` on a GPU machine. matkit itself +does not import pygRASPA at setup time — it only writes input files and a +small ``run.py`` that invokes pygRASPA when executed. + +Mirrors the API of :mod:`matkit.graspa.graspa` so the two backends are +interchangeable apart from the ML-specific keyword arguments. +""" + +from __future__ import annotations + +import json +import shutil +from concurrent.futures import ThreadPoolExecutor +from itertools import product +from pathlib import Path + +from ase.io import read as ase_read + +from matkit.graspa.graspa import generate_component_blocks +from matkit.utils.cif import sanitize_cif_stem +from matkit.utils.template import copy_template, render_template +from matkit.utils.unitcell_calculator import calculate_cell_size + +VALID_MODES = ("run", "run-auto") + + +def _render_run_py( + *, + mode: str, + model_path: str, + model_type: str, + task: str | None, + E_comps: list[float], + save_poscar: bool, + log_file: str = "output.log", + checkpoint_file: str = "checkpoint.json", + stop_flag: str = "STOP_FLAG", +) -> str: + """Return the contents of the per-simulation ``run.py`` launcher.""" + if mode not in VALID_MODES: + raise ValueError(f"mode must be one of {VALID_MODES}, got {mode!r}") + task_repr = repr(task) if task is None else repr(task) + ecomps_repr = "[" + ", ".join(repr(float(e)) for e in E_comps) + "]" + + header = ( + '"""Auto-generated by matkit.pygraspa.\n\n' + "Run with `python run.py` (single-process) or wrap in" + " your scheduler.\nFor run-auto mode, touch STOP_FLAG" + ' to checkpoint and exit gracefully.\n"""\n\n' + "from pygRASPA.utils import read_pseudo_atoms\n" + ) + + common_kwargs = ( + f" model_path={model_path!r},\n" + f" model_type={model_type!r},\n" + f" task={task_repr},\n" + f' simulation_input="simulation.input",\n' + f" E_comps={ecomps_repr},\n" + f" save_trial_poscar={bool(save_poscar)!r},\n" + f' type_index_to_label=read_pseudo_atoms("pseudo_atoms.def"),\n' + ) + + if mode == "run": + body = ( + "from pygRASPA.main import run\n\n" + "run(\n" + f"{common_kwargs}" + f' log_file_path="{log_file}",\n' + ")\n" + ) + else: # run-auto + body = ( + "from pygRASPA.auto.main import run_auto\n\n" + "run_auto(\n" + f"{common_kwargs}" + f' log_file="{log_file}",\n' + f' checkpoint_file="{checkpoint_file}",\n' + f' stop_flag="{stop_flag}",\n' + ")\n" + ) + + return header + body + + +def setup_simulation( + cif: str, + outpath: str, + adsorbates: list[dict], + *, + model_path: str, + model_type: str, + E_comps: list[float], + task: str | None = None, + mode: str = "run-auto", + save_poscar: bool = False, + temperature: float = 298.0, + pressure: float = 1e5, + cutoff: float = 12.8, + n_cycle: int = 1000, + template_dir: str = "template", + cell_size: list[int] | None = None, +) -> bool: + """Set up a pygRASPA (ML-potential GCMC) simulation directory. + + Writes the standard gRASPA input files plus a ``run.py`` launcher that + invokes ``pygRASPA.main.run`` (``mode='run'``) or + ``pygRASPA.auto.main.run_auto`` (``mode='run-auto'``) at run time. + + Args: + cif: Path to the input CIF structure file. + outpath: Directory to write simulation files to. + adsorbates: List of adsorbate dicts with at least ``MoleculeName``. + Length must equal ``len(E_comps)``. + model_path: Path to the ML model checkpoint (``.pt`` or deployed + NequIP bundle). + model_type: Backend identifier; one of pygRASPA's ``MODEL_TYPES`` + (e.g. ``FAIRChem-esen``, ``FAIRChem-uma``, ``mace_polar``). + E_comps: Pre-computed ML energy of each isolated rigid adsorbate + molecule, in eV. One entry per adsorbate. + task: Task name required by FAIRChem-uma / FAIRChem-AllScAIP + backends. Ignored otherwise. + mode: ``'run'`` for plain ML MC, ``'run-auto'`` for checkpointed + cluster-friendly run. + save_poscar: Pass through to pygRASPA's ``save_trial_poscar``. + temperature: Simulation temperature in Kelvin. + pressure: Simulation pressure in Pascals. + cutoff: Van der Waals cutoff radius in Angstrom. + n_cycle: Number of Monte Carlo cycles. + template_dir: Template subdirectory under ``files/``. + cell_size: Pre-computed unit cell dimensions ``[uc_x, uc_y, uc_z]``. + If provided, skips re-reading the CIF. + + Returns: + ``True`` on success. + + Raises: + FileNotFoundError: If the CIF file does not exist. + ValueError: If ``len(adsorbates) != len(E_comps)``, or mode is + invalid, or the FAIRChem task is missing when required. + """ + if len(adsorbates) != len(E_comps): + raise ValueError( + f"adsorbates ({len(adsorbates)}) and E_comps ({len(E_comps)}) " + "must have the same length" + ) + if model_type in ("FAIRChem-uma", "FAIRChem-AllScAIP") and not task: + raise ValueError(f"task is required when model_type is '{model_type}'") + + outdir = Path(outpath) + cifpath = Path(cif) + if not cifpath.exists(): + raise FileNotFoundError(f"CIF file does not exist: {cif}") + safe_stem = sanitize_cif_stem(cifpath.stem) + + template_path = Path(__file__).parent / "files" / template_dir + copy_template(template_path, outdir) + shutil.copy(cifpath, outdir / f"{safe_stem}.cif") + + if cell_size is not None: + uc_x, uc_y, uc_z = cell_size + else: + atoms = ase_read(cifpath) + uc_x, uc_y, uc_z = calculate_cell_size(atoms) + + input_path = outdir / "simulation.input" + render_template( + input_path, + { + "NCYCLE": str(n_cycle), + "TEMPERATURE": str(temperature), + "PRESSURE": str(pressure), + "CUTOFF": str(cutoff), + "CIFFILE": safe_stem, + "UC_X": str(uc_x), + "UC_Y": str(uc_y), + "UC_Z": str(uc_z), + }, + ) + + content = input_path.read_text() + if "__COMPONENTS__" in content: + component_block = generate_component_blocks(adsorbates) + content = content.replace("__COMPONENTS__", component_block) + input_path.write_text(content) + + (outdir / "run.py").write_text( + _render_run_py( + mode=mode, + model_path=model_path, + model_type=model_type, + task=task, + E_comps=E_comps, + save_poscar=save_poscar, + ) + ) + + return True + + +def _setup_single_cif( + cif: Path, + out_path: Path, + adsorbates: list[dict], + temperatures: list[float], + pressures: list[float], + cutoff: float, + n_cycle: int, + template_dir: str, + model_path: str, + model_type: str, + E_comps: list[float], + task: str | None, + mode: str, + save_poscar: bool, +) -> list[dict]: + """Set up all T x P pygRASPA simulations for one CIF file.""" + atoms = ase_read(cif) + cell_size = calculate_cell_size(atoms) + safe_stem = sanitize_cif_stem(cif.stem) + + entries = [] + for temp, pres in product(temperatures, pressures): + sim_dir = out_path / safe_stem / f"T{temp}_P{pres:g}" + setup_simulation( + cif=str(cif), + outpath=str(sim_dir), + adsorbates=adsorbates, + model_path=model_path, + model_type=model_type, + E_comps=E_comps, + task=task, + mode=mode, + save_poscar=save_poscar, + temperature=temp, + pressure=pres, + cutoff=cutoff, + n_cycle=n_cycle, + template_dir=template_dir, + cell_size=cell_size, + ) + entry = { + "sim_dir": str(sim_dir), + "cif": cif.name, + "temperature": temp, + "pressure": pres, + "adsorbates": [ad["MoleculeName"] for ad in adsorbates], + "model_path": model_path, + "model_type": model_type, + "mode": mode, + } + if safe_stem != cif.stem: + entry["original_cif_stem"] = cif.stem + entry["safe_cif_stem"] = safe_stem + entries.append(entry) + return entries + + +def setup_batch( + cif_dir: str, + outpath: str, + adsorbates: list[dict], + *, + model_path: str, + model_type: str, + E_comps: list[float], + temperatures: list[float], + pressures: list[float], + task: str | None = None, + mode: str = "run-auto", + save_poscar: bool = False, + cutoff: float = 12.8, + n_cycle: int = 1000, + template_dir: str = "template", + max_workers: int | None = None, +) -> list[dict]: + """Set up pygRASPA simulations for all CIF x T x P combinations. + + Mirrors :func:`matkit.graspa.graspa.setup_batch`. Writes a + ``simulations.jsonl`` manifest and (if any CIFs were renamed) + a ``cif_mapping.json``. + + See :func:`setup_simulation` for argument semantics. + """ + cif_path = Path(cif_dir) + out_path = Path(outpath) + + if not cif_path.is_dir(): + raise ValueError(f"CIF directory does not exist: {cif_dir}") + + cif_files = sorted(cif_path.glob("*.cif")) + if not cif_files: + raise ValueError(f"No .cif files found in {cif_dir}") + + manifest = [] + with ThreadPoolExecutor(max_workers=max_workers) as pool: + futures = [ + pool.submit( + _setup_single_cif, + cif, + out_path, + adsorbates, + temperatures, + pressures, + cutoff, + n_cycle, + template_dir, + model_path, + model_type, + E_comps, + task, + mode, + save_poscar, + ) + for cif in cif_files + ] + for future in futures: + manifest.extend(future.result()) + + manifest_path = out_path / "simulations.jsonl" + with manifest_path.open("w") as f: + for entry in manifest: + f.write(json.dumps(entry) + "\n") + + rename_map = { + entry["original_cif_stem"]: entry["safe_cif_stem"] + for entry in manifest + if "original_cif_stem" in entry + } + if rename_map: + _write_cif_mapping(out_path, rename_map) + + return manifest + + +def _write_cif_mapping(out_path: Path, mapping: dict[str, str]) -> None: + """Write/merge {original_stem: safe_stem} to ``cif_mapping.json``.""" + mapping_path = out_path / "cif_mapping.json" + existing: dict[str, str] = {} + if mapping_path.exists(): + existing = json.loads(mapping_path.read_text()) + existing.update(mapping) + mapping_path.write_text(json.dumps(existing, indent=2, sort_keys=True)) + + +def compute_ecomp( + adsorbate_def: str, + model_path: str, + model_type: str, + *, + task: str | None = None, + device: str = "cuda", + cache_path: str | None = None, +) -> float: + """Compute E_comp (ML energy of one isolated rigid adsorbate molecule). + + Wraps :func:`pygRASPA.utils.read_molecule` and + :func:`pygRASPA.ml_setup.setup_calculator`. Used to obtain the per- + component reference energy that pygRASPA requires for the Metropolis + criterion. Results may be cached to a JSON file so repeated isotherm + sweeps with the same model/adsorbate avoid recomputation. + + Args: + adsorbate_def: Path to a gRASPA ``*.def`` molecule file. + model_path: ML model checkpoint path. + model_type: pygRASPA model-type identifier. + task: Task name (required for FAIRChem-uma / FAIRChem-AllScAIP). + device: ``"cuda"`` or ``"cpu"``. + cache_path: Optional JSON file for memoizing results keyed by + ``(model_path, model_type, task, adsorbate stem)``. + + Returns: + E_comp in eV. + + Raises: + ImportError: If pygRASPA is not installed. + FileNotFoundError: If ``adsorbate_def`` does not exist. + """ + def_path = Path(adsorbate_def) + if not def_path.exists(): + raise FileNotFoundError( + f"Adsorbate def file not found: {adsorbate_def}" + ) + + cache_key = f"{model_path}|{model_type}|{task or ''}|{def_path.stem}" + if cache_path is not None: + cp = Path(cache_path) + if cp.exists(): + cache = json.loads(cp.read_text()) + if cache_key in cache: + return float(cache[cache_key]) + + try: + from pygRASPA.ml_setup import setup_calculator + from pygRASPA.utils import read_molecule + except ImportError as e: + raise ImportError( + "pygRASPA is required for compute_ecomp. " + "Install pygRASPA on the GPU machine first." + ) from e + + atoms = read_molecule(str(def_path)) + atoms.info = {"charge": 0, "spin": 1} + atoms.calc = setup_calculator( + model_path, device, model_type=model_type, task=task + ) + value = float(atoms.get_potential_energy()) + + if cache_path is not None: + cp = Path(cache_path) + cache = json.loads(cp.read_text()) if cp.exists() else {} + cache[cache_key] = value + cp.write_text(json.dumps(cache, indent=2, sort_keys=True)) + + return value + + +def _convert_uptake( + molecules_per_uc: float, + error_per_uc: float, + cif_path: str, + n_uc: tuple[int, int, int], + unit: str, + molar_mass: float | None = None, +) -> tuple[float, float]: + """Convert molecules-per-supercell to a chemistry-friendly uptake unit. + + Args: + molecules_per_uc: Mean adsorbed molecule count over the full + ``n_uc`` supercell. + error_per_uc: Standard error on the same scale. + cif_path: Path to the framework CIF (used for mass / volume). + n_uc: ``(uc_x, uc_y, uc_z)`` replication factors. + unit: ``"mol/kg"``, ``"mg/g"``, or ``"g/L"``. + molar_mass: Adsorbate molar mass (g/mol), required for ``mg/g`` + and ``g/L``. + + Returns: + ``(uptake, error)`` in the requested unit. + """ + atoms = ase_read(cif_path) + n_cells = n_uc[0] * n_uc[1] * n_uc[2] + framework_mass = sum(atoms.get_masses()) * n_cells # g/mol + if unit == "mol/kg": + return ( + molecules_per_uc / framework_mass * 1000, + error_per_uc / framework_mass * 1000, + ) + if molar_mass is None: + raise ValueError(f"molar_mass is required for unit={unit!r}") + if unit == "mg/g": + return ( + molecules_per_uc * molar_mass / framework_mass * 1000, + error_per_uc * molar_mass / framework_mass * 1000, + ) + if unit == "g/L": + vol_l = atoms.get_volume() * n_cells * 1e-27 + return ( + molecules_per_uc / 6.022e23 * molar_mass / vol_l, + error_per_uc / 6.022e23 * molar_mass / vol_l, + ) + raise ValueError(f"Unsupported unit: {unit!r}") + + +_MOLAR_MASS = { + "CO2": 44.0098, + "N2": 28.0134, + "H2": 2.01588, + "CH4": 16.0425, + "methane": 16.0425, + "H2O": 18.0153, + "TIP4P": 18.0153, + "SO2": 64.066, + "Xe": 131.293, + "Kr": 83.798, +} + + +def get_output_data( + output_path: str, + *, + unit: str = "mol/kg", + output_fname: str = "output.log", + n_skip_cycles: int | None = None, + cifname: str | None = None, +) -> dict: + """Parse a pygRASPA per-cycle JSON-lines log and average uptake. + + pygRASPA writes one JSON object per Monte Carlo cycle to the log file + specified via ``log_file_path`` (``run`` mode) or ``log_file`` + (``run-auto`` mode). Each line is expected to carry the current cycle + index and per-component molecule counts. + + The exact field names vary across pygRASPA versions; this parser is + tolerant and looks for the first plausible key in each category + (``cycle``/``step``/``i_cycle`` for the index; + ``molecules``/``N_components``/``N`` for counts). If the schema in your + log differs, raise an issue and pass the field names through + pygRASPA's own ``readCycles``/utils helpers. + + Args: + output_path: Directory containing the simulation output. + unit: ``"mol/kg"``, ``"mg/g"``, or ``"g/L"``. + output_fname: Log filename (default ``output.log``). + n_skip_cycles: Number of leading cycles to drop. Defaults to + ``NumberOfInitializationCycles + NumberOfEquilibrationCycles`` + from ``simulation.input``. + cifname: Framework CIF name in ``output_path``. Auto-detected if + ``None``. + + Returns: + Dict with keys: ``success``, ``uptake`` (per-component list), + ``error`` (per-component list), ``unit``, ``n_components``, + ``n_production_cycles``. + + Raises: + ValueError: If the log cannot be parsed or required metadata is + missing. + """ + output_dir = Path(output_path) + log_path = output_dir / output_fname + if not log_path.exists(): + raise ValueError(f"Output log not found: {log_path}") + + cycle_keys = ("cycle", "step", "i_cycle", "Cycle") + count_keys = ("molecules", "N_components", "N", "n_molecules") + + cycle_records: list[tuple[int, list[float]]] = [] + with log_path.open() as fh: + for raw in fh: + raw = raw.strip() + if not raw or not raw.startswith("{"): + continue + try: + rec = json.loads(raw) + except json.JSONDecodeError: + continue + cycle = next((rec[k] for k in cycle_keys if k in rec), None) + counts = next((rec[k] for k in count_keys if k in rec), None) + if cycle is None or counts is None: + continue + if isinstance(counts, (int, float)): + counts = [float(counts)] + else: + counts = [float(c) for c in counts] + cycle_records.append((int(cycle), counts)) + + if not cycle_records: + raise ValueError( + f"No parseable cycle records in {log_path}. " + "Confirm the file is the pygRASPA JSON-lines log." + ) + + if n_skip_cycles is None: + sim_input = output_dir / "simulation.input" + n_skip_cycles = ( + _read_skip_cycles(sim_input) if sim_input.exists() else 0 + ) + + production = [c for c in cycle_records if c[0] >= n_skip_cycles] + if not production: + raise ValueError( + f"No cycles >= n_skip_cycles={n_skip_cycles} in {log_path}" + ) + + n_components = len(production[0][1]) + sums = [0.0] * n_components + sq_sums = [0.0] * n_components + for _, counts in production: + if len(counts) != n_components: + raise ValueError( + "Inconsistent component count across cycles " + f"({n_components} vs {len(counts)})" + ) + for i, v in enumerate(counts): + sums[i] += v + sq_sums[i] += v * v + n = len(production) + means = [s / n for s in sums] + variances = [max(0.0, sq / n - m * m) for sq, m in zip(sq_sums, means)] + sems = [(v / n) ** 0.5 for v in variances] + + # Read unit cell replication from simulation.input. + n_uc = _read_unit_cells(output_dir / "simulation.input") + + if cifname is None: + cifs = list(output_dir.glob("*.cif")) + if len(cifs) != 1: + raise ValueError( + f"Found {len(cifs)} CIFs in {output_path}; pass cifname." + ) + cif_path = str(cifs[0]) + else: + cif_path = str(output_dir / cifname) + + uptakes: list[float] = [] + errors: list[float] = [] + for i in range(n_components): + molar_mass = None + if unit in ("mg/g", "g/L"): + ads = _component_name_from_input(output_dir / "simulation.input", i) + molar_mass = _MOLAR_MASS.get(ads) + if molar_mass is None: + raise ValueError( + f"No molar mass for component {i} ({ads!r}); " + "extend _MOLAR_MASS or request mol/kg." + ) + u, e = _convert_uptake( + means[i], sems[i], cif_path, n_uc, unit, molar_mass + ) + uptakes.append(u) + errors.append(e) + + return { + "success": True, + "uptake": uptakes if n_components > 1 else uptakes[0], + "error": errors if n_components > 1 else errors[0], + "unit": unit, + "n_components": n_components, + "n_production_cycles": n, + } + + +def _read_skip_cycles(sim_input: Path) -> int: + init = equil = 0 + for line in sim_input.read_text().splitlines(): + parts = line.split() + if len(parts) >= 2 and parts[0] == "NumberOfInitializationCycles": + init = int(parts[1]) + elif len(parts) >= 2 and parts[0] == "NumberOfEquilibrationCycles": + equil = int(parts[1]) + return init + equil + + +def _read_unit_cells(sim_input: Path) -> tuple[int, int, int]: + for line in sim_input.read_text().splitlines(): + parts = line.split() + if parts and parts[0] == "UnitCells" and len(parts) >= 5: + # "UnitCells 0 ux uy uz" (box_index ux uy uz) + return int(parts[2]), int(parts[3]), int(parts[4]) + return 1, 1, 1 + + +def _component_name_from_input(sim_input: Path, index: int) -> str | None: + """Return the MoleculeName for component *index* in simulation.input.""" + if not sim_input.exists(): + return None + for line in sim_input.read_text().splitlines(): + parts = line.split() + if ( + len(parts) >= 4 + and parts[0] == "Component" + and parts[1] == str(index) + and parts[2] == "MoleculeName" + ): + return parts[3] + return None diff --git a/tests/test_graspa.py b/tests/test_graspa.py index 2b977b1..980c388 100644 --- a/tests/test_graspa.py +++ b/tests/test_graspa.py @@ -192,7 +192,9 @@ def test_batch_writes_mapping_only_for_renamed( mapping = json.loads(mapping_path.read_text()) assert mapping == {"weird.v2": "weird_v2"} - assert (out_dir / "weird_v2" / "T298.0_P100000" / "weird_v2.cif").exists() + assert ( + out_dir / "weird_v2" / "T298.0_P100000" / "weird_v2.cif" + ).exists() assert (out_dir / "clean" / "T298.0_P100000" / "clean.cif").exists() def test_batch_no_mapping_when_all_clean(self, sample_cif, tmp_path): diff --git a/tests/test_pygraspa.py b/tests/test_pygraspa.py new file mode 100644 index 0000000..4025620 --- /dev/null +++ b/tests/test_pygraspa.py @@ -0,0 +1,265 @@ +"""Tests for matkit.pygraspa module.""" + +import json +import shutil +from pathlib import Path + +import pytest + +from matkit.pygraspa.pygraspa import ( + _render_run_py, + get_output_data, + setup_batch, + setup_simulation, +) + + +class TestRenderRunPy: + """Verify the per-simulation run.py launcher is rendered correctly.""" + + def test_run_mode(self): + out = _render_run_py( + mode="run", + model_path="/m/model.pt", + model_type="FAIRChem-esen", + task=None, + E_comps=[-22.97915], + save_poscar=False, + ) + assert "from pygRASPA.main import run" in out + assert "run(" in out + assert "/m/model.pt" in out + assert "FAIRChem-esen" in out + assert "-22.97915" in out + assert "task=None" in out + assert "save_trial_poscar=False" in out + assert 'log_file_path="output.log"' in out + assert "run_auto" not in out + + def test_run_auto_mode(self): + out = _render_run_py( + mode="run-auto", + model_path="/m/model.pt", + model_type="FAIRChem-uma", + task="ads_eng", + E_comps=[-22.0, -17.5], + save_poscar=True, + ) + assert "from pygRASPA.auto.main import run_auto" in out + assert "run_auto(" in out + assert "'ads_eng'" in out + assert "[-22.0, -17.5]" in out + assert "save_trial_poscar=True" in out + assert 'checkpoint_file="checkpoint.json"' in out + assert 'stop_flag="STOP_FLAG"' in out + + def test_invalid_mode_raises(self): + with pytest.raises(ValueError, match="mode must be one of"): + _render_run_py( + mode="invalid", + model_path="x", + model_type="FAIRChem-esen", + task=None, + E_comps=[1.0], + save_poscar=False, + ) + + +class TestSetupSimulation: + """End-to-end setup writes simulation.input, run.py, and CIF.""" + + def _run(self, cif, outdir, **overrides): + defaults = dict( + cif=cif, + outpath=str(outdir), + adsorbates=[{"MoleculeName": "CO2"}], + model_path="/path/to/model.pt", + model_type="FAIRChem-esen", + E_comps=[-22.97915], + ) + defaults.update(overrides) + return setup_simulation(**defaults) + + def test_creates_simulation_input(self, sample_cif, tmp_path): + out = tmp_path / "sim" + self._run(sample_cif, out, temperature=300.0, pressure=2e5, n_cycle=500) + assert (out / "simulation.input").exists() + content = (out / "simulation.input").read_text() + assert "300.0" in content + assert "200000.0" in content + assert "500" in content + assert "NCYCLE" not in content + assert "TEMPERATURE" not in content + assert "__COMPONENTS__" not in content + + def test_writes_run_py_with_baked_args(self, sample_cif, tmp_path): + out = tmp_path / "sim" + self._run(sample_cif, out, mode="run-auto") + run_py = (out / "run.py").read_text() + assert "/path/to/model.pt" in run_py + assert "FAIRChem-esen" in run_py + assert "-22.97915" in run_py + assert "from pygRASPA.auto.main import run_auto" in run_py + + def test_run_py_is_valid_python(self, sample_cif, tmp_path): + out = tmp_path / "sim" + self._run(sample_cif, out) + compile((out / "run.py").read_text(), str(out / "run.py"), "exec") + + def test_copies_cif(self, sample_cif, tmp_path): + out = tmp_path / "sim" + self._run(sample_cif, out) + assert (out / (Path(sample_cif).stem + ".cif")).exists() + + def test_mismatched_ecomps_raises(self, sample_cif, tmp_path): + with pytest.raises(ValueError, match="same length"): + self._run( + sample_cif, + tmp_path / "sim", + adsorbates=[{"MoleculeName": "CO2"}, {"MoleculeName": "N2"}], + E_comps=[-22.0], + ) + + def test_uma_requires_task(self, sample_cif, tmp_path): + with pytest.raises(ValueError, match="task is required"): + self._run( + sample_cif, + tmp_path / "sim", + model_type="FAIRChem-uma", + task=None, + ) + + def test_missing_cif_raises(self, tmp_path): + with pytest.raises(FileNotFoundError): + self._run("/nonexistent.cif", tmp_path / "sim") + + def test_sanitizes_multi_period_cif(self, sample_cif, tmp_path): + weird = tmp_path / "foo.bar.cif" + shutil.copy(sample_cif, weird) + out = tmp_path / "sim" + self._run(str(weird), out) + assert (out / "foo_bar.cif").exists() + content = (out / "simulation.input").read_text() + assert "FrameworkName foo_bar" in content + + +class TestSetupBatch: + """Batch setup mirrors graspa.setup_batch structure.""" + + def test_batch_manifest_and_rename(self, sample_cif, tmp_path): + cif_dir = tmp_path / "cifs" + cif_dir.mkdir() + shutil.copy(sample_cif, cif_dir / "clean.cif") + shutil.copy(sample_cif, cif_dir / "weird.v2.cif") + + out_dir = tmp_path / "batch" + manifest = setup_batch( + cif_dir=str(cif_dir), + outpath=str(out_dir), + adsorbates=[{"MoleculeName": "CO2"}], + model_path="/m.pt", + model_type="FAIRChem-esen", + E_comps=[-22.97915], + temperatures=[298.0], + pressures=[1e5], + n_cycle=10, + ) + + assert len(manifest) == 2 + assert (out_dir / "simulations.jsonl").exists() + assert (out_dir / "cif_mapping.json").exists() + + # Each sim dir should have its own run.py + for entry in manifest: + assert (Path(entry["sim_dir"]) / "run.py").exists() + assert (Path(entry["sim_dir"]) / "simulation.input").exists() + assert entry["model_path"] == "/m.pt" + assert entry["model_type"] == "FAIRChem-esen" + assert entry["mode"] == "run-auto" + + mapping = json.loads((out_dir / "cif_mapping.json").read_text()) + assert mapping == {"weird.v2": "weird_v2"} + + +class TestGetOutputData: + """JSON-lines log parser averages production-cycle counts.""" + + def _write_log(self, dir_, n_init, n_total, n_comp): + """Write a fake JSON-lines log: N=10*cycle for each cycle.""" + lines = [] + for c in range(1, n_total + 1): + counts = [10.0 * c + i for i in range(n_comp)] + lines.append(json.dumps({"cycle": c, "molecules": counts})) + (dir_ / "output.log").write_text("\n".join(lines) + "\n") + + def _write_sim_input(self, dir_, n_init, n_prod, ux=2, uy=2, uz=2): + (dir_ / "simulation.input").write_text( + f"NumberOfInitializationCycles {n_init}\n" + "NumberOfEquilibrationCycles 0\n" + f"NumberOfProductionCycles {n_prod}\n" + f"UnitCells 0 {ux} {uy} {uz}\n" + "Component 0 MoleculeName CO2\n" + ) + + def _write_cif(self, dir_, sample_cif): + shutil.copy(sample_cif, dir_ / "frame.cif") + + def test_parses_single_component(self, sample_cif, tmp_path): + d = tmp_path / "out" + d.mkdir() + self._write_log(d, n_init=2, n_total=5, n_comp=1) + self._write_sim_input(d, n_init=2, n_prod=3) + self._write_cif(d, sample_cif) + + result = get_output_data(str(d), unit="mol/kg") + assert result["success"] is True + assert result["n_components"] == 1 + assert result["n_production_cycles"] == 4 # cycles 2,3,4,5 + assert isinstance(result["uptake"], float) + assert result["unit"] == "mol/kg" + + def test_parses_mixture(self, sample_cif, tmp_path): + d = tmp_path / "out" + d.mkdir() + self._write_log(d, n_init=1, n_total=4, n_comp=2) + self._write_sim_input(d, n_init=1, n_prod=3) + # Need both components in input for mol/kg path (mg/g would need names) + (d / "simulation.input").write_text( + "NumberOfInitializationCycles 1\n" + "NumberOfEquilibrationCycles 0\n" + "NumberOfProductionCycles 3\n" + "UnitCells 0 1 1 1\n" + "Component 0 MoleculeName CO2\n" + "Component 1 MoleculeName N2\n" + ) + self._write_cif(d, sample_cif) + + result = get_output_data(str(d), unit="mol/kg") + assert result["n_components"] == 2 + assert isinstance(result["uptake"], list) + assert len(result["uptake"]) == 2 + + def test_skip_cycles_override(self, sample_cif, tmp_path): + d = tmp_path / "out" + d.mkdir() + self._write_log(d, n_init=0, n_total=5, n_comp=1) + self._write_sim_input(d, n_init=0, n_prod=5) + self._write_cif(d, sample_cif) + + full = get_output_data(str(d), n_skip_cycles=0) + partial = get_output_data(str(d), n_skip_cycles=4) + # With strictly increasing counts, dropping early cycles raises mean. + assert partial["uptake"] > full["uptake"] + + def test_missing_log_raises(self, tmp_path): + d = tmp_path / "out" + d.mkdir() + with pytest.raises(ValueError, match="Output log not found"): + get_output_data(str(d)) + + def test_unparseable_log_raises(self, tmp_path): + d = tmp_path / "out" + d.mkdir() + (d / "output.log").write_text("not json\nalso not json\n") + with pytest.raises(ValueError, match="No parseable cycle records"): + get_output_data(str(d)) diff --git a/tests/test_utils.py b/tests/test_utils.py index 66b3d22..e0123cd 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -4,7 +4,6 @@ import numpy as np from pathlib import Path from ase import Atoms -from ase.io import write as ase_write from matkit.utils.unitcell_calculator import calculate_cell_size from matkit.utils.remove_solvent import remove_solvent @@ -48,7 +47,7 @@ def test_custom_cutoff(self): atoms.set_pbc(True) small_cutoff = calculate_cell_size(atoms, cutoff=4.0) large_cutoff = calculate_cell_size(atoms, cutoff=12.8) - assert all(s <= l for s, l in zip(small_cutoff, large_cutoff)) + assert all(s <= L for s, L in zip(small_cutoff, large_cutoff)) def test_returns_list_of_ints(self): """Return type should be a list of integers."""