{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/tymworld/DROPPS/blob/v1.0.0/examples/colab/01_DROPPS_1_0_Model_and_System_Builder.ipynb)\n\n# DROPPS 1.0 Molecular and System Builder\n\nBuild a compact **cubic** box containing many copies of one protein for the\nphase-coexistence workflow. Multi-component systems remain available as an\nextension. Optional PTMs, angle restraints, elastic networks, residue\nmodifications, and deterministic packing match the construction workflow\ndescribed in the manuscript.\n\nDo not elongate the initial box here. The simulation notebook first condenses\nthis compact system under NPT, then expands one axis tenfold and starts the NVT\nslab-production stage.\n\nAll operations use the public `dps` CLI. Bead IDs in angle and elastic input\nfiles are **1-based**; DROPPS index groups used for trajectory analysis are\n**0-based** and are handled in later notebooks.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Install DROPPS 1.0\n#@markdown The default route installs the immutable wheel attached to the\n#@markdown public GitHub `v1.0.0` release. Upload remains available for offline use.\ninstallation_source = \"Install from the DROPPS v1.0.0 GitHub release\" #@param [\"Install from the DROPPS v1.0.0 GitHub release\", \"Upload the DROPPS 1.0 wheel\", \"Install from another wheel URL\"]\nwheel_url = \"https://github.com/tymworld/DROPPS/releases/download/v1.0.0/dropps-1.0-py3-none-any.whl\" #@param {type:\"string\"}\naccelerator_dependencies = \"auto\" #@param [\"auto\", \"cuda12\", \"cuda13\", \"latest\"]\n\nimport importlib.metadata\nimport re\nimport shutil\nimport subprocess\nimport sys\nfrom pathlib import Path\n\n\ndef detect_install_extra():\n    if accelerator_dependencies != \"auto\":\n        return accelerator_dependencies\n    if shutil.which(\"nvidia-smi\"):\n        probe = subprocess.run(\n            [\"nvidia-smi\"], text=True, capture_output=True, check=False\n        ).stdout\n        match = re.search(r\"CUDA Version:\\s*(\\d+)\", probe)\n        if match and int(match.group(1)) >= 13:\n            return \"cuda13\"\n        return \"cuda12\"\n    return \"latest\"\n\n\nif installation_source.startswith(\"Upload\"):\n    try:\n        from google.colab import files\n    except ImportError as exc:\n        raise RuntimeError(\n            \"This upload form is intended for Google Colab. Set installation_source \"\n            \"to the URL option when running elsewhere.\"\n        ) from exc\n    uploaded = files.upload()\n    wheel_candidates = [Path(name) for name in uploaded if name.endswith(\".whl\")]\n    if len(wheel_candidates) != 1:\n        raise ValueError(\"Upload exactly one DROPPS 1.0 .whl file.\")\n    wheel_target = str(wheel_candidates[0].resolve())\nelse:\n    if not wheel_url.strip():\n        raise ValueError(\"Provide the published DROPPS 1.0 wheel URL.\")\n    wheel_target = wheel_url.strip()\n\nextra = detect_install_extra()\nif wheel_target.startswith((\"https://\", \"http://\")):\n    package_spec = f\"dropps[{extra}] @ {wheel_target}\"\nelse:\n    package_spec = f\"{wheel_target}[{extra}]\"\nprint(f\"Installing {package_spec}\")\nsubprocess.run(\n    [\n        sys.executable,\n        \"-m\",\n        \"pip\",\n        \"install\",\n        \"--quiet\",\n        \"--upgrade\",\n        \"--upgrade-strategy\",\n        \"only-if-needed\",\n        package_spec,\n    ],\n    check=True,\n)\n\nversion = importlib.metadata.version(\"dropps\")\nif version != \"1.0\":\n    raise RuntimeError(f\"Expected DROPPS 1.0, but installed {version}.\")\nsubprocess.run([\"dps\", \"--version\"], check=True)\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Shared notebook helpers\nimport hashlib\nimport json\nimport os\nimport shlex\nimport subprocess\nimport zipfile\nfrom datetime import datetime, timezone\nfrom pathlib import Path\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nCOMMAND_LOG = []\n\n\ndef run_command(arguments, *, cwd=None, input_text=None, check=True):\n    command = [str(value) for value in arguments]\n    print(\"$\", shlex.join(command))\n    result = subprocess.run(\n        command,\n        cwd=None if cwd is None else str(cwd),\n        input=input_text,\n        text=True,\n        capture_output=True,\n        check=False,\n    )\n    if result.stdout:\n        print(result.stdout, end=\"\" if result.stdout.endswith(\"\\n\") else \"\\n\")\n    if result.stderr:\n        print(result.stderr, end=\"\" if result.stderr.endswith(\"\\n\") else \"\\n\")\n    COMMAND_LOG.append(\n        {\n            \"time_utc\": datetime.now(timezone.utc).isoformat(),\n            \"cwd\": str(Path(cwd or Path.cwd()).resolve()),\n            \"command\": command,\n            \"returncode\": result.returncode,\n        }\n    )\n    if check and result.returncode:\n        raise RuntimeError(\n            f\"Command failed with exit code {result.returncode}: {shlex.join(command)}\"\n        )\n    return result\n\n\ndef dps(*arguments, cwd=None, input_text=None, check=True):\n    return run_command(\n        [\"dps\", *arguments], cwd=cwd, input_text=input_text, check=check\n    )\n\n\ndef reset_task_directory(path):\n    path = Path(path).resolve()\n    if not path.name.startswith(\"dropps_\"):\n        raise ValueError(f\"Refusing to reset unexpected directory: {path}\")\n    if path.exists():\n        import shutil\n\n        shutil.rmtree(path)\n    path.mkdir(parents=True)\n    return path\n\n\ndef sha256(path):\n    digest = hashlib.sha256()\n    with Path(path).open(\"rb\") as stream:\n        for block in iter(lambda: stream.read(1024 * 1024), b\"\"):\n            digest.update(block)\n    return digest.hexdigest()\n\n\ndef safe_extract_zip(archive, destination):\n    destination = Path(destination).resolve()\n    destination.mkdir(parents=True, exist_ok=True)\n    with zipfile.ZipFile(archive) as bundle:\n        for member in bundle.infolist():\n            target = (destination / member.filename).resolve()\n            if destination not in target.parents and target != destination:\n                raise ValueError(f\"Unsafe ZIP member: {member.filename}\")\n        bundle.extractall(destination)\n\n\ndef find_unique(root, basename):\n    matches = [path for path in Path(root).rglob(basename) if path.is_file()]\n    if len(matches) != 1:\n        raise FileNotFoundError(\n            f\"Expected one file named {basename!r} below {root}, found {len(matches)}.\"\n        )\n    return matches[0]\n\n\ndef make_zip(paths, output_path, *, base=None):\n    output_path = Path(output_path)\n    base = Path(base or output_path.parent).resolve()\n    with zipfile.ZipFile(output_path, \"w\", compression=zipfile.ZIP_DEFLATED) as bundle:\n        for source in sorted({Path(path).resolve() for path in paths}):\n            if source == output_path.resolve() or not source.is_file():\n                continue\n            try:\n                arcname = source.relative_to(base)\n            except ValueError:\n                arcname = Path(source.name)\n            bundle.write(source, arcname)\n    print(f\"Wrote {output_path} ({output_path.stat().st_size / 1024:.1f} KiB)\")\n    return output_path\n\n\ndef download(path):\n    try:\n        from google.colab import files\n    except ImportError:\n        print(f\"Result available at {Path(path).resolve()}\")\n    else:\n        files.download(str(path))\n\n\ndef optional_time_arguments(start, end, interval):\n    arguments = []\n    for option, value in ((\"-b\", start), (\"-e\", end), (\"-dt\", interval)):\n        if str(value).strip():\n            arguments.extend([option, str(value).strip()])\n    return arguments\n\n\ndef read_xvg(path):\n    legends = {}\n    metadata = {}\n    rows = []\n    with Path(path).open(encoding=\"utf-8\") as stream:\n        for raw in stream:\n            line = raw.strip()\n            if not line or line.startswith(\"#\"):\n                continue\n            if line.startswith(\"@\"):\n                re_module = __import__(\"re\")\n                match = re_module.search(r's(\\d+)\\s+legend\\s+\"(.*)\"', line)\n                if match:\n                    legends[int(match.group(1)) + 1] = match.group(2)\n                for key, pattern in {\n                    \"title\": r'@\\s+title\\s+\"(.*)\"',\n                    \"xlabel\": r'@\\s+xaxis\\s+label\\s+\"(.*)\"',\n                    \"ylabel\": r'@\\s+yaxis\\s+label\\s+\"(.*)\"',\n                }.items():\n                    label_match = re_module.search(pattern, line)\n                    if label_match:\n                        metadata[key] = label_match.group(1)\n                continue\n            rows.append([float(value) for value in line.split()])\n    data = np.asarray(rows, dtype=float)\n    if data.ndim != 2 or data.shape[1] < 2:\n        raise ValueError(f\"No plottable data in {path}\")\n    return data, legends, metadata\n\n\ndef plot_xvg(path, *, title=None, output=None):\n    data, legends, metadata = read_xvg(path)\n    fig, ax = plt.subplots(figsize=(7.2, 4.2))\n    for column in range(1, data.shape[1]):\n        ax.plot(data[:, 0], data[:, column], label=legends.get(column, f\"series {column}\"))\n    ax.set_title(title or metadata.get(\"title\", Path(path).stem))\n    ax.set_xlabel(metadata.get(\"xlabel\", \"x / time\"))\n    ax.set_ylabel(metadata.get(\"ylabel\", \"value\"))\n    if data.shape[1] > 2:\n        ax.legend(frameon=False, fontsize=8)\n    ax.spines[[\"top\", \"right\"]].set_visible(False)\n    fig.tight_layout()\n    output_path = Path(output) if output is not None else Path(path).with_suffix(\".png\")\n    fig.savefig(output_path, dpi=220, bbox_inches=\"tight\")\n    print(f\"Saved figure: {output_path}\")\n    plt.show()\n    return fig, output_path\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 1 Define up to three molecular components\n\nThe default is a lightweight 80-copy, single-protein construction example,\nfollowing the manuscript's 5×5×5 placement grid. Replace the sequence with\nthe target protein for a scientific run. Mapping uses Cα coordinates from the\nfirst model of an uploaded all-atom PDB. Enable Components B/C only for a\nco-phase-separation system and use compatible force-field families.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Component A\ncomponent_a_enabled = True #@param {type:\"boolean\"}\ncomponent_a_name = \"SCAFFOLD\" #@param {type:\"string\"}\ncomponent_a_sequence = \"FWFWFWFWFWFWFWFW\" #@param {type:\"string\"}\ncomponent_a_count = 80 #@param {type:\"integer\"}\ncomponent_a_forcefield = \"HPS\" #@param [\"HPS\", \"HPST\", \"CALVADOS2\", \"HPSRNA\", \"MPiPi\", \"MPiPi_PTM\"]\ncomponent_a_input_pdb = \"\" #@param {type:\"string\"}\ncomponent_a_ptms = \"\" #@param {type:\"string\"}\ncomponent_a_radius_nm = 2.0 #@param {type:\"number\"}\ncomponent_a_extension = 0.5 #@param {type:\"slider\", min:0, max:1, step:0.05}\ncomponent_a_residue_index = 1 #@param {type:\"integer\"}\ncomponent_a_charge_ntd = False #@param {type:\"boolean\"}\ncomponent_a_charge_ctd = False #@param {type:\"boolean\"}\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Component B\ncomponent_b_enabled = False #@param {type:\"boolean\"}\ncomponent_b_name = \"CLIENT\" #@param {type:\"string\"}\ncomponent_b_sequence = \"KKGGEEDD\" #@param {type:\"string\"}\ncomponent_b_count = 4 #@param {type:\"integer\"}\ncomponent_b_forcefield = \"HPS\" #@param [\"HPS\", \"HPST\", \"CALVADOS2\", \"HPSRNA\", \"MPiPi\", \"MPiPi_PTM\"]\ncomponent_b_input_pdb = \"\" #@param {type:\"string\"}\ncomponent_b_ptms = \"\" #@param {type:\"string\"}\ncomponent_b_radius_nm = 2.0 #@param {type:\"number\"}\ncomponent_b_extension = 0.5 #@param {type:\"slider\", min:0, max:1, step:0.05}\ncomponent_b_residue_index = 1 #@param {type:\"integer\"}\ncomponent_b_charge_ntd = False #@param {type:\"boolean\"}\ncomponent_b_charge_ctd = False #@param {type:\"boolean\"}\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Component C\ncomponent_c_enabled = False #@param {type:\"boolean\"}\ncomponent_c_name = \"RNA\" #@param {type:\"string\"}\ncomponent_c_sequence = \"AAAAAAAAAA\" #@param {type:\"string\"}\ncomponent_c_count = 2 #@param {type:\"integer\"}\ncomponent_c_forcefield = \"HPSRNA\" #@param [\"HPS\", \"HPST\", \"CALVADOS2\", \"HPSRNA\", \"MPiPi\", \"MPiPi_PTM\"]\ncomponent_c_input_pdb = \"\" #@param {type:\"string\"}\ncomponent_c_ptms = \"\" #@param {type:\"string\"}\ncomponent_c_radius_nm = 2.0 #@param {type:\"number\"}\ncomponent_c_extension = 0.5 #@param {type:\"slider\", min:0, max:1, step:0.05}\ncomponent_c_residue_index = 1 #@param {type:\"integer\"}\ncomponent_c_charge_ntd = False #@param {type:\"boolean\"}\ncomponent_c_charge_ctd = False #@param {type:\"boolean\"}\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Upload optional all-atom PDB reference structures\nupload_reference_pdbs = False #@param {type:\"boolean\"}\nrandom_seed = 1215 #@param {type:\"integer\"}\n\nBASE = Path(\"/content\") if Path(\"/content\").is_dir() else Path.cwd()\nWORKDIR = reset_task_directory(BASE / \"dropps_model_builder\")\nUPLOAD_DIR = WORKDIR / \"uploads\"\nUPLOAD_DIR.mkdir()\n\nif upload_reference_pdbs:\n    try:\n        from google.colab import files\n    except ImportError as exc:\n        raise RuntimeError(\"Upload reference structures in Google Colab.\") from exc\n    import shutil\n\n    for name in files.upload():\n        source = Path(name)\n        if source.suffix.lower() != \".pdb\":\n            raise ValueError(f\"Expected PDB input, received {source}\")\n        shutil.copy2(source, UPLOAD_DIR / source.name)\n\nprint(\"Working directory:\", WORKDIR)\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Build monomer structures and self-contained ITP topologies\nimport re\n\ndef component(prefix):\n    scope = globals()\n    return {\n        \"enabled\": scope[f\"component_{prefix}_enabled\"],\n        \"name\": scope[f\"component_{prefix}_name\"].strip(),\n        \"sequence\": scope[f\"component_{prefix}_sequence\"].strip(),\n        \"count\": int(scope[f\"component_{prefix}_count\"]),\n        \"forcefield\": scope[f\"component_{prefix}_forcefield\"],\n        \"input_pdb\": scope[f\"component_{prefix}_input_pdb\"].strip(),\n        \"ptms\": scope[f\"component_{prefix}_ptms\"].strip(),\n        \"radius\": float(scope[f\"component_{prefix}_radius_nm\"]),\n        \"extension\": float(scope[f\"component_{prefix}_extension\"]),\n        \"residue_index\": int(scope[f\"component_{prefix}_residue_index\"]),\n        \"charge_ntd\": scope[f\"component_{prefix}_charge_ntd\"],\n        \"charge_ctd\": scope[f\"component_{prefix}_charge_ctd\"],\n    }\n\ncomponents = [item for item in map(component, \"abc\") if item[\"enabled\"]]\nif not components:\n    raise ValueError(\"Enable at least one component.\")\nif len({item[\"name\"] for item in components}) != len(components):\n    raise ValueError(\"Component names must be unique.\")\n\nfor item in components:\n    if not re.fullmatch(r\"[A-Za-z][A-Za-z0-9_]*\", item[\"name\"]):\n        raise ValueError(f\"Use a simple alphanumeric molecule name: {item['name']!r}\")\n    if not item[\"sequence\"]:\n        raise ValueError(f\"Sequence is empty for {item['name']}\")\n    arguments = [\n        \"pdb2dps\", \"-s\", item[\"sequence\"], \"-ff\", item[\"forcefield\"],\n        \"-oc\", f\"{item['name']}.pdb\", \"-op\", f\"{item['name']}.itp\",\n        \"-on\", item[\"name\"], \"-ri\", item[\"residue_index\"],\n        \"-r\", item[\"radius\"], \"-e\", item[\"extension\"],\n        \"--seed\", random_seed,\n    ]\n    if item[\"input_pdb\"]:\n        source = Path(item[\"input_pdb\"])\n        if not source.is_file():\n            source = UPLOAD_DIR / item[\"input_pdb\"]\n        if not source.is_file():\n            raise FileNotFoundError(source)\n        arguments.extend([\"-f\", source])\n    if item[\"ptms\"]:\n        arguments.extend([\"-ptm\", *item[\"ptms\"].split()])\n    if item[\"charge_ntd\"]:\n        arguments.append(\"-cNTD\")\n    if item[\"charge_ctd\"]:\n        arguments.append(\"-cCTD\")\n    dps(*arguments, cwd=WORKDIR)\n    item[\"pdb\"] = WORKDIR / f\"{item['name']}.pdb\"\n    item[\"itp\"] = WORKDIR / f\"{item['name']}.itp\"\n\nprint(\"Built:\", \", \".join(item[\"name\"] for item in components))\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 2 Optional topology refinements\n\nEnter semicolon-separated records. Angle records use\n`CENTER_RESIDUE ANGLE_DEG FORCE_CONSTANT`. Elastic groups use compact\nranges such as `1-77; 106-177; 192-256`; the notebook expands them into\nthe explicit 1-based bead lists required by `genelastic`.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Angle restraints with dps addangle\nangle_target = \"SCAFFOLD\" #@param {type:\"string\"}\nangle_records = \"\" #@param {type:\"string\"}\n\nif angle_records.strip():\n    target = next(item for item in components if item[\"name\"] == angle_target)\n    angle_file = WORKDIR / f\"{angle_target}_angles.dat\"\n    angle_file.write_text(\n        \"\\n\".join(part.strip() for part in angle_records.split(\";\") if part.strip()) + \"\\n\",\n        encoding=\"utf-8\",\n    )\n    output = WORKDIR / f\"{angle_target}_angles.itp\"\n    dps(\n        \"addangle\", \"-ip\", target[\"itp\"], \"-op\", output,\n        \"-al\", angle_file, cwd=WORKDIR,\n    )\n    target[\"itp\"] = output\nelse:\n    print(\"No angle restraints requested.\")\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Elastic networks with dps genelastic\nelastic_target = \"SCAFFOLD\" #@param {type:\"string\"}\nelastic_groups = \"\" #@param {type:\"string\"}\nelastic_force_constant = 5000.0 #@param {type:\"number\"}\nelastic_lower_nm = 0.5 #@param {type:\"number\"}\nelastic_upper_nm = 0.9 #@param {type:\"number\"}\n\ndef expand_integer_ranges(expression):\n    values = []\n    for token in expression.replace(\",\", \" \").split():\n        if \"-\" in token:\n            start, stop = map(int, token.split(\"-\", 1))\n            if stop < start:\n                raise ValueError(f\"Descending range: {token}\")\n            values.extend(range(start, stop + 1))\n        else:\n            values.append(int(token))\n    return values\n\nif elastic_groups.strip():\n    target = next(item for item in components if item[\"name\"] == elastic_target)\n    lines = []\n    for group in elastic_groups.split(\";\"):\n        if group.strip():\n            lines.append(\" \".join(map(str, expand_integer_ranges(group))))\n    elastic_file = WORKDIR / f\"{elastic_target}_elastic.dat\"\n    elastic_file.write_text(\"\\n\".join(lines) + \"\\n\", encoding=\"utf-8\")\n    output = WORKDIR / f\"{elastic_target}_elastic.itp\"\n    dps(\n        \"genelastic\", \"-f\", target[\"pdb\"], \"-p\", target[\"itp\"],\n        \"-o\", output, \"-er\", elastic_file,\n        \"-ef\", elastic_force_constant, \"-el\", elastic_lower_nm,\n        \"-eu\", elastic_upper_nm, cwd=WORKDIR,\n    )\n    target[\"itp\"] = output\nelse:\n    print(\"No elastic network requested.\")\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Optional post-build residue modification with dps modifyres\napply_residue_modifications = False #@param {type:\"boolean\"}\nmodification_target = \"SCAFFOLD\" #@param {type:\"string\"}\nresidue_modifications = \"S3SMP\" #@param {type:\"string\"}\nmodification_forcefield = \"MPiPi_PTM\" #@param [\"HPS\", \"HPST\", \"CALVADOS2\", \"HPSRNA\", \"MPiPi\", \"MPiPi_PTM\"]\n\nif apply_residue_modifications:\n    target = next(item for item in components if item[\"name\"] == modification_target)\n    modified_pdb = WORKDIR / f\"{modification_target}_modified.pdb\"\n    modified_itp = WORKDIR / f\"{modification_target}_modified.itp\"\n    dps(\n        \"modifyres\", \"-ip\", target[\"itp\"], \"-if\", target[\"pdb\"],\n        \"-op\", modified_itp, \"-of\", modified_pdb,\n        \"-ff\", modification_forcefield,\n        \"-m\", *residue_modifications.split(), cwd=WORKDIR,\n    )\n    target[\"pdb\"], target[\"itp\"] = modified_pdb, modified_itp\nelse:\n    print(\"No post-build residue modification requested.\")\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 3 Pack the compact cubic simulation box\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title System packing parameters\nmesh_x = 5 #@param {type:\"integer\"}\nmesh_y = 5 #@param {type:\"integer\"}\nmesh_z = 5 #@param {type:\"integer\"}\nminimum_gap_nm = 1.0 #@param {type:\"number\"}\nbox_type = \"cubic\" #@param [\"cubic\", \"anisotropy\", \"xy\"]\nminimum_x_nm = 5.0 #@param {type:\"number\"}\nminimum_y_nm = 5.0 #@param {type:\"number\"}\nminimum_z_nm = 5.0 #@param {type:\"number\"}\nshuffle_components = True #@param {type:\"boolean\"}\nnon_cubic_molecule = False #@param {type:\"boolean\"}\n\ntotal_count = sum(item[\"count\"] for item in components)\nif mesh_x * mesh_y * mesh_z < total_count:\n    raise ValueError(\"The mesh has fewer sites than requested molecules.\")\n\narguments = [\"genmesh\", \"-f\"]\narguments.extend(item[\"pdb\"] for item in components)\narguments.append(\"-p\")\narguments.extend(item[\"itp\"] for item in components)\narguments.append(\"-n\")\narguments.extend(item[\"count\"] for item in components)\narguments.extend(\n    [\n        \"-mesh\", mesh_x, mesh_y, mesh_z,\n        \"-g\", minimum_gap_nm, \"-bt\", box_type,\n        \"-mx\", minimum_x_nm, \"-my\", minimum_y_nm, \"-mz\", minimum_z_nm,\n        \"--seed\", random_seed, \"-oc\", \"system.pdb\", \"-op\", \"system.top\",\n    ]\n)\nif shuffle_components:\n    arguments.append(\"-s\")\nif non_cubic_molecule:\n    arguments.append(\"-ncm\")\ndps(*arguments, cwd=WORKDIR)\nSYSTEM_PDB = WORKDIR / \"system.pdb\"\nSYSTEM_TOP = WORKDIR / \"system.top\"\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "The exported `system.pdb` is deliberately compact and cubic. Its box must be\nelongated only after the NPT stage has produced a condensed `npt.pdb`.\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 4 Inspect and export the model bundle\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Preview bead coordinates\ndef pdb_coordinates_nm(path):\n    rows = []\n    for line in Path(path).read_text(encoding=\"utf-8\").splitlines():\n        if line.startswith((\"ATOM  \", \"HETATM\")):\n            rows.append([float(line[30:38]), float(line[38:46]), float(line[46:54])])\n    return np.asarray(rows) / 10.0\n\ncoordinates = pdb_coordinates_nm(SYSTEM_PDB)\nfig = plt.figure(figsize=(6.5, 5.5))\nax = fig.add_subplot(111, projection=\"3d\")\nax.scatter(*coordinates.T, s=8, alpha=0.7)\nax.set(xlabel=\"x (nm)\", ylabel=\"y (nm)\", zlabel=\"z (nm)\", title=\"Initial system\")\nfig.tight_layout()\nplt.show()\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "#@title Save manifest and download\nselected_files = [SYSTEM_PDB, SYSTEM_TOP]\nfor item in components:\n    selected_files.extend([item[\"pdb\"], item[\"itp\"]])\nselected_files = sorted({Path(path).resolve() for path in selected_files})\nmanifest = {\n    \"notebook\": \"01_DROPPS_1_0_Model_and_System_Builder.ipynb\",\n    \"dropps_version\": \"1.0\",\n    \"seed\": random_seed,\n    \"components\": [\n        {\n            key: value\n            for key, value in item.items()\n            if key not in {\"pdb\", \"itp\", \"input_pdb\"}\n        }\n        | {\"pdb\": item[\"pdb\"].name, \"itp\": item[\"itp\"].name}\n        for item in components\n    ],\n    \"system_pdb\": SYSTEM_PDB.name,\n    \"system_top\": SYSTEM_TOP.name,\n    \"files\": {path.name: sha256(path) for path in selected_files},\n    \"commands\": COMMAND_LOG,\n}\nmanifest_path = WORKDIR / \"model_manifest.json\"\nmanifest_path.write_text(json.dumps(manifest, indent=2), encoding=\"utf-8\")\nselected_files.append(manifest_path)\narchive = make_zip(\n    selected_files, WORKDIR / \"DROPPS_1_0_model_bundle.zip\", base=WORKDIR\n)\ndownload(archive)\n"
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}