{ "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/04_DROPPS_1_0_Phase_Separation_Analysis.ipynb)\n\n# DROPPS 1.0 Phase Separation Analysis\n\nAnalyse the **elongated-box NVT trajectory**, never the compact-box NPT\ntrajectory. The notebook verifies that pressure coupling is disabled before\nrunning the analysis chain emphasized in the manuscript: density profiles;\nresidue contact maps and residue-class statistics; chain size and local angle;\nintra/inter-chain distances; MSD; and spontaneous-condensation assembly size,\ncomposition, radius of gyration, asphericity, and ellipticity. The notebook\nalso exposes the auxiliary 1.0 RMSD command.\n\nSelect groups explicitly to avoid interactive prompts. Default group numbering\nis `System = 0`, followed by molecule types in TOP order. Use the trajectory\npreparation notebook to create and inspect custom NDX groups.\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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Upload or reuse a workflow bundle\ninput_source = \"Upload ZIP or individual files\" #@param [\"Upload ZIP or individual files\", \"Reuse files from an earlier notebook in this runtime\"]\nreuse_directory = \"/content/dropps_trajectory\" #@param {type:\"string\"}\n\nimport shutil\n\nBASE = Path(\"/content\") if Path(\"/content\").is_dir() else Path.cwd()\nINPUT_ROOT = reset_task_directory(BASE / \"dropps_uploaded_inputs\")\n\nif input_source.startswith(\"Upload\"):\n try:\n from google.colab import files\n except ImportError as exc:\n raise RuntimeError(\"Upload this bundle in Colab or choose the reuse option.\") from exc\n uploaded = files.upload()\n for name in uploaded:\n source = Path(name)\n if source.suffix.lower() == \".zip\":\n safe_extract_zip(source, INPUT_ROOT)\n else:\n shutil.copy2(source, INPUT_ROOT / source.name)\nelse:\n source_root = Path(reuse_directory)\n if not source_root.is_dir():\n raise FileNotFoundError(source_root)\n for source in source_root.rglob(\"*\"):\n if source.is_file() and source.suffix.lower() in {\n \".pdb\", \".itp\", \".top\", \".mdp\", \".tpr\", \".xtc\", \".ndx\", \".chk\", \".xml\", \".json\"\n }:\n target = INPUT_ROOT / source.name\n if target.exists() and sha256(target) != sha256(source):\n raise ValueError(f\"Duplicate filename with different content: {source.name}\")\n shutil.copy2(source, target)\n\nprint(\"Available files:\")\nfor path in sorted(INPUT_ROOT.rglob(\"*\")):\n if path.is_file():\n print(\" \", path.relative_to(INPUT_ROOT), f\"({path.stat().st_size / 1024:.1f} KiB)\")\n" }, { "cell_type": "markdown", "metadata": {}, "source": "## 1 Resolve files and analysis selections\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Input filenames and index groups\nrun_input_filename = \"slab_nvt.tpr\" #@param {type:\"string\"}\ntrajectory_filename = \"processed.xtc\" #@param {type:\"string\"}\nindex_filename = \"analysis.ndx\" #@param {type:\"string\"}\nreference_group = 1 #@param {type:\"integer\"}\ntwo_component_system = False #@param {type:\"boolean\"}\nselection_group = 2 #@param {type:\"integer\"}\ndensity_center_group = 0 #@param {type:\"integer\"}\nslab_axis = \"z\"\nstart_time_ns = \"\" #@param {type:\"string\"}\nend_time_ns = \"\" #@param {type:\"string\"}\nsampling_interval_ns = \"\" #@param {type:\"string\"}\n\nTPR = find_unique(INPUT_ROOT, run_input_filename)\nTRAJECTORY = find_unique(INPUT_ROOT, trajectory_filename)\n\nimport zipfile\n\nif not zipfile.is_zipfile(TPR):\n raise ValueError(\"Expected a portable DROPPS 1.0 TPR v2 file.\")\nwith zipfile.ZipFile(TPR) as archive:\n tpr_parameters = json.loads(archive.read(\"parameters.json\"))\nif tpr_parameters.get(\"pcoulp\"):\n raise ValueError(\n \"Phase-coexistence analysis requires slab_nvt.tpr (pcoulp=False), \"\n \"not the compact-box NPT run input.\"\n )\ntry:\n NDX = find_unique(INPUT_ROOT, index_filename)\nexcept FileNotFoundError:\n NDX = None\n\nBASE = Path(\"/content\") if Path(\"/content\").is_dir() else Path.cwd()\nWORKDIR = reset_task_directory(BASE / \"dropps_analysis\")\nif NDX is None:\n NDX = WORKDIR / \"analysis.ndx\"\n result = dps(\n \"make_ndx\", \"-s\", TPR, \"-o\", NDX,\n cwd=WORKDIR, input_text=\"q\\n\", check=False,\n )\n if not NDX.is_file():\n raise RuntimeError(\n f\"make_ndx did not create {NDX}; subprocess status was {result.returncode}.\"\n )\ndps(\"check\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX, cwd=WORKDIR)\nprint(\"Verified analysis ensemble: NVT (pcoulp = False)\")\n\nTIME_ARGS = optional_time_arguments(\n start_time_ns, end_time_ns, sampling_interval_ns\n)\nANALYSIS_STATUS = []\n\ndef run_analysis(label, arguments):\n try:\n dps(*arguments, cwd=WORKDIR)\n except Exception as exc:\n ANALYSIS_STATUS.append({\"analysis\": label, \"status\": \"failed\", \"error\": str(exc)})\n print(f\"[{label}] FAILED: {exc}\")\n return False\n ANALYSIS_STATUS.append({\"analysis\": label, \"status\": \"complete\"})\n return True\n" }, { "cell_type": "markdown", "metadata": {}, "source": "## 2 Phase behavior\n\n`density` can output each component separately while centering on a shared\nreference density. `contact` supports a global cutoff or residue-specific\nσ-scaled cutoff. `cstat` reproduces the residue-class aggregation used for\ninteraction interpretation in the manuscript.\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Required z-axis density profile and PNG figure\ndensity_axis = \"z\"\ndensity_type = \"mass\" #@param [\"mass\", \"charge\"]\ndensity_bin_width_nm = 0.05 #@param {type:\"number\"}\ndensity_center_mode = \"frame\" #@param [\"frame\", \"block\", \"none\"]\ndensity_blocks = 5 #@param {type:\"integer\"}\ndense_phase_threshold = 0.5 #@param {type:\"number\"}\n\nif density_axis != \"z\" or slab_axis != \"z\":\n raise ValueError(\"Phase-coexistence density analysis must use the z axis.\")\ncalculate_groups = [reference_group]\nif two_component_system:\n calculate_groups.append(selection_group)\narguments = [\n \"density\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-o\", \"density.xvg\", \"-x\", \"z\", \"-tp\", density_type,\n \"--bin-width\", density_bin_width_nm,\n \"-selfit\", density_center_group, \"-sel\", *calculate_groups,\n \"--center-mode\", density_center_mode, \"--blocks\", density_blocks,\n \"-t\", dense_phase_threshold, *TIME_ARGS,\n]\ndensity_ok = run_analysis(\"density\", arguments)\nif not density_ok:\n raise RuntimeError(\"Required z-axis density analysis failed.\")\n_, DENSITY_FIGURE = plot_xvg(\n WORKDIR / \"density.xvg\",\n title=f\"{density_type.title()} density along z\",\n output=WORKDIR / \"density_z.png\",\n)\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Contact maps and contact-number time series\nrun_contacts = True #@param {type:\"boolean\"}\ncontact_cutoff_scheme = \"global\" #@param [\"global\", \"residue\"]\nglobal_contact_cutoff_nm = 0.7 #@param {type:\"number\"}\nresidue_cutoff_multiplier = 1.2 #@param {type:\"number\"}\nremove_intra_diagonals = 2 #@param {type:\"integer\"}\ncontact_use_pbc = True #@param {type:\"boolean\"}\naverage_intra_chain_maps = True #@param {type:\"boolean\"}\n\ncontact_ok = False\nif run_contacts:\n arguments = [\n \"contact\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-ref\", reference_group,\n \"-sel\", selection_group if two_component_system else reference_group,\n \"-cs\", contact_cutoff_scheme,\n \"-c\", global_contact_cutoff_nm, \"-cm\", residue_cutoff_multiplier,\n \"-rd\", remove_intra_diagonals, \"-otype\", \"dat\",\n \"-orr\", \"contact_reference_reference.dat\",\n \"-or\", \"contact_reference_intra.dat\",\n \"-otrr\", \"contact_number_reference_reference.xvg\",\n *TIME_ARGS,\n ]\n if two_component_system:\n arguments.extend(\n [\n \"-sel\", selection_group,\n \"-ors\", \"contact_reference_selection.dat\",\n \"-otrs\", \"contact_number_reference_selection.xvg\",\n ]\n )\n if contact_use_pbc:\n arguments.append(\"-pbc\")\n if average_intra_chain_maps:\n arguments.append(\"-intraavg\")\n contact_ok = run_analysis(\"contact\", arguments)\n\n if contact_ok:\n for path in sorted(WORKDIR.glob(\"contact_*.dat\")):\n matrix = np.loadtxt(path)\n fig, ax = plt.subplots(figsize=(5.2, 4.4))\n image = ax.imshow(matrix, origin=\"lower\", aspect=\"auto\", cmap=\"magma\")\n ax.set(title=path.stem, xlabel=\"selection residue\", ylabel=\"reference residue\")\n fig.colorbar(image, ax=ax, label=\"contact value\")\n fig.tight_layout()\n figure_path = WORKDIR / f\"{path.stem}.png\"\n fig.savefig(figure_path, dpi=220, bbox_inches=\"tight\")\n print(\"Saved figure:\", figure_path)\n plt.show()\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Residue-class contact statistics with dps cstat\nrun_contact_statistics = True #@param {type:\"boolean\"}\ncontact_grouping_scheme = \"HPST_SC\" #@param [\"HPST\", \"AHCP\", \"HCP\", \"HPST_SC\", \"AHCP_SC\", \"HCP_SC\"]\ncontact_aggregation = \"sum\" #@param [\"sum\", \"average\"]\n\nif run_contact_statistics:\n if not contact_ok:\n print(\"Contact statistics skipped because no contact map was generated in this run.\")\n else:\n map_name = (\n \"contact_reference_selection.dat\"\n if two_component_system\n else \"contact_reference_reference.dat\"\n )\n selected = selection_group if two_component_system else reference_group\n run_analysis(\n \"cstat\",\n [\n \"cstat\", \"-m\", map_name, \"-s\", TPR, \"-n\", NDX,\n \"-ref\", reference_group, \"-sel\", selected,\n \"-o\", \"contact_statistics.xlsx\", \"-gr\",\n \"-gs\", contact_grouping_scheme, \"-ag\", contact_aggregation,\n ],\n )\n" }, { "cell_type": "markdown", "metadata": {}, "source": "## 3 Chain conformation and dynamics\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Radius of gyration and backbone angles\nrun_gyrate = True #@param {type:\"boolean\"}\nrun_angle = True #@param {type:\"boolean\"}\ngyrate_histogram_bin_width_nm = 0.1 #@param {type:\"number\"}\nconformational_use_pbc = True #@param {type:\"boolean\"}\n\nif run_gyrate:\n arguments = [\n \"gyrate\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-sel\", reference_group, \"-oa\", \"gyrate_average.xvg\",\n \"-ov\", \"gyrate_per_chain.xvg\", \"-oh\", \"gyrate_histogram.xvg\",\n \"-bw\", gyrate_histogram_bin_width_nm, *TIME_ARGS,\n ]\n if conformational_use_pbc:\n arguments.append(\"-pbc\")\n if run_analysis(\"gyrate\", arguments):\n plot_xvg(WORKDIR / \"gyrate_average.xvg\", title=\"Mean chain radius of gyration\")\n\nif run_angle:\n arguments = [\n \"angle\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-sel\", reference_group, \"-ot\", \"angle_time.xvg\",\n \"-or\", \"angle_by_residue.xvg\", \"-ors\", \"angle_statistics.xvg\",\n *TIME_ARGS,\n ]\n if conformational_use_pbc:\n arguments.append(\"-pbc\")\n if run_analysis(\"angle\", arguments):\n plot_xvg(WORKDIR / \"angle_by_residue.xvg\", title=\"Backbone angle by residue\")\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Intra- and inter-chain distances\nrun_intra_chain_distance = False #@param {type:\"boolean\"}\nintra_pair_group_ids = \"\" #@param {type:\"string\"}\nrun_inter_chain_distance = False #@param {type:\"boolean\"}\ninter_reference_group = 0 #@param {type:\"integer\"}\ninter_selection_group = 0 #@param {type:\"integer\"}\ndistance_use_pbc = True #@param {type:\"boolean\"}\n\nif run_intra_chain_distance:\n pair_groups = [int(value) for value in intra_pair_group_ids.replace(\",\", \" \").split()]\n if not pair_groups:\n raise ValueError(\"Provide one or more NDX groups containing one bead pair per chain.\")\n arguments = [\n \"idist\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-sel\", *pair_groups, \"-ot\", \"idist_time.xvg\",\n \"-op\", \"idist_pairs.xvg\", \"-ops\", \"idist_statistics.xvg\",\n *TIME_ARGS,\n ]\n if distance_use_pbc:\n arguments.append(\"-pbc\")\n run_analysis(\"idist\", arguments)\n\nif run_inter_chain_distance:\n arguments = [\n \"odist\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-ref\", inter_reference_group, \"-sel\", inter_selection_group,\n \"-oa\", \"odist_average.xvg\", \"-ov\", \"odist_all_pairs.xvg\",\n *TIME_ARGS,\n ]\n if distance_use_pbc:\n arguments.append(\"-pbc\")\n run_analysis(\"odist\", arguments)\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Mean-square displacement and RMSD\nrun_msd = True #@param {type:\"boolean\"}\nmsd_group = 1 #@param {type:\"integer\"}\nmsd_dimensions = \"xyz\" #@param [\"xyz\", \"xy\", \"yz\", \"xz\", \"x\", \"y\", \"z\"]\nrun_rmsd = False #@param {type:\"boolean\"}\nrmsd_selection_expression = \"(mol SCAFFOLD)\" #@param {type:\"string\"}\nrmsd_fit_mode = \"molecule\" #@param [\"none\", \"molecule\", \"selection\"]\nrmsd_output_mode = \"molecule\" #@param [\"molecule\", \"selection\"]\nrmsd_mass_weighted = True #@param {type:\"boolean\"}\n\nif run_msd:\n arguments = [\n \"msd\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-sel\", msd_group, \"-t\", msd_dimensions,\n \"-o\", \"msd.xvg\", *TIME_ARGS,\n ]\n if run_analysis(\"msd\", arguments):\n plot_xvg(WORKDIR / \"msd.xvg\", title=\"Mean-square displacement\")\n\nif run_rmsd:\n arguments = [\n \"rmsd\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-o\", \"rmsd.xvg\", \"--select\", rmsd_selection_expression,\n \"--fit-mode\", rmsd_fit_mode, \"--output-mode\", rmsd_output_mode,\n *TIME_ARGS,\n ]\n if rmsd_mass_weighted:\n arguments.append(\"--mass-weighted\")\n if run_analysis(\"rmsd\", arguments):\n plot_xvg(WORKDIR / \"rmsd.xvg\", title=\"Molecular RMSD\")\n" }, { "cell_type": "markdown", "metadata": {}, "source": "## 4 Spontaneous-condensation assembly analysis\n\nA cluster is a maximal set of chains connected through intermolecular bead\ncontacts. Tune the bead cutoff and the minimum contact count to the model.\nThe molecule-fraction and shape cutoff excludes very small clusters from\ncomposition and morphology statistics.\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Cluster formation, composition, and shape\nrun_assembly = True #@param {type:\"boolean\"}\nassembly_reference_group = 0 #@param {type:\"integer\"}\nassembly_contact_cutoff_nm = 0.7 #@param {type:\"number\"}\ncontacts_to_connect_chains = 5 #@param {type:\"integer\"}\nlarge_cluster_size_cutoff = 10 #@param {type:\"integer\"}\nassembly_use_pbc = True #@param {type:\"boolean\"}\n\nif run_assembly:\n selection_groups = [reference_group]\n if two_component_system:\n selection_groups.append(selection_group)\n arguments = [\n \"assembly\", \"-s\", TPR, \"-f\", TRAJECTORY, \"-n\", NDX,\n \"-ref\", assembly_reference_group, \"-sel\", *selection_groups,\n \"-c\", assembly_contact_cutoff_nm, \"-t\", contacts_to_connect_chains,\n \"-mfc\", large_cluster_size_cutoff,\n \"-cn\", \"assembly_cluster_number.xvg\",\n \"-cs\", \"assembly_largest_size.xvg\",\n \"-csd\", \"assembly_size_distribution.xvg\",\n \"-mf\", \"assembly_molecule_fraction.xvg\",\n \"-rgl\", \"assembly_largest_rg.xvg\",\n \"-asp\", \"assembly_asphericity.xvg\",\n \"-elp\", \"assembly_ellipticity.xvg\",\n *TIME_ARGS,\n ]\n if assembly_use_pbc:\n arguments.append(\"-pbc\")\n if run_analysis(\"assembly\", arguments):\n plot_xvg(WORKDIR / \"assembly_cluster_number.xvg\", title=\"Assembly count\")\n plot_xvg(WORKDIR / \"assembly_largest_size.xvg\", title=\"Largest assembly size\")\n" }, { "cell_type": "markdown", "metadata": {}, "source": "## 5 Analysis manifest and results archive\n" }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "#@title Summarize and download\nmanifest = {\n \"notebook\": \"04_DROPPS_1_0_Phase_Separation_Analysis.ipynb\",\n \"dropps_version\": \"1.0\",\n \"source_tpr\": {\"name\": TPR.name, \"sha256\": sha256(TPR)},\n \"source_trajectory\": {\"name\": TRAJECTORY.name, \"sha256\": sha256(TRAJECTORY)},\n \"source_index\": {\"name\": NDX.name, \"sha256\": sha256(NDX)},\n \"verified_ensemble\": \"NVT\",\n \"slab_axis\": slab_axis,\n \"groups\": {\n \"reference\": reference_group,\n \"selection\": selection_group if two_component_system else None,\n \"density_center\": density_center_group,\n },\n \"time_window_ns\": {\n \"start\": start_time_ns,\n \"end\": end_time_ns,\n \"interval\": sampling_interval_ns,\n },\n \"analysis_status\": ANALYSIS_STATUS,\n \"commands\": COMMAND_LOG,\n}\nmanifest_path = WORKDIR / \"analysis_manifest.json\"\nmanifest_path.write_text(json.dumps(manifest, indent=2), encoding=\"utf-8\")\nexport_files = [path for path in WORKDIR.iterdir() if path.is_file()]\narchive = make_zip(\n export_files,\n WORKDIR / \"DROPPS_1_0_analysis_results.zip\",\n base=WORKDIR,\n)\nprint(json.dumps(ANALYSIS_STATUS, indent=2))\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 }