{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating Geometries for Geant4\n",
    "\n",
    "This notebook contains code blocks used to generate sample\n",
    "geometries used in Geant4 simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "try:\n",
    "    from fractaldna.dna_models import dnachain\n",
    "except (ImportError, ModuleNotFoundError):\n",
    "    sys.path.append(str(Path.cwd().parent.parent.parent))\n",
    "    from fractaldna.dna_models import dnachain\n",
    "\n",
    "import numpy as np\n",
    "from mayavi import mlab\n",
    "\n",
    "mlab.options.offscreen = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single DNA Segments\n",
    "\n",
    "### Straight and Turned Segments for a 50nm box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "bp_separation = dnachain.BP_SEPARATION  # 3.32Å\n",
    "side_length_nm = 50  # nm\n",
    "num_basepairs_straight = int(side_length_nm / (0.1 * bp_separation))\n",
    "num_basepairs_turned = int((side_length_nm * np.pi / 4.0) / (0.1 * bp_separation))\n",
    "\n",
    "chain_straight = dnachain.DNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_straight))\n",
    ")\n",
    "\n",
    "chain_turned = dnachain.TurnedDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned))\n",
    ")\n",
    "\n",
    "chain_turned_twisted = dnachain.TurnedTwistedDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned))\n",
    ")\n",
    "\n",
    "chain_straight.to_frame().to_csv(\"results/50nm_straight.csv\", sep=\" \", index=False)\n",
    "chain_turned.to_frame().to_csv(\"results/50nm_turn.csv\", sep=\" \", index=False)\n",
    "chain_turned_twisted.to_frame().to_csv(\n",
    "    \"results/50nm_turn_twist.csv\", sep=\" \", index=False\n",
    ")\n",
    "\n",
    "chain_straight.to_plot().savefig(\"results/50nm_straight.png\", sep=\" \", index=False)\n",
    "chain_turned.to_plot().savefig(\"results/50nm_turn.png\", sep=\" \", index=False)\n",
    "chain_turned_twisted.to_plot().savefig(\n",
    "    \"results/50nm_turn_twist.png\", sep=\" \", index=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='results/50nm_straight.png' width='30%'> <img src='results/50nm_turn.png' width='30%'> <img src='results/50nm_turn_twist.png' width='30%'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Multi Strand straight and turned segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "bp_separation = dnachain.BP_SEPARATION  # 3.32Å\n",
    "side_length_nm = 50  # nm\n",
    "num_basepairs_straight = int(side_length_nm / (0.1 * bp_separation))\n",
    "num_basepairs_turned = int((side_length_nm * np.pi / 4.0) / (0.1 * bp_separation))\n",
    "strand_separation = 100  # angstroms\n",
    "\n",
    "chain4_straight = dnachain.FourStrandDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_straight)),\n",
    "    strand_separation,\n",
    ")\n",
    "\n",
    "chain4_turned = dnachain.FourStrandTurnedDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned)),\n",
    "    strand_separation,\n",
    ")\n",
    "\n",
    "chain4_turned_twisted = dnachain.FourStrandTurnedDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned)),\n",
    "    strand_separation,\n",
    "    twist=True,\n",
    ")\n",
    "\n",
    "chain4_straight.to_frame().to_csv(\"results/50nm_4_straight.csv\", sep=\" \", index=False)\n",
    "chain4_turned.to_frame().to_csv(\"results/50nm_4_turn.csv\", sep=\" \", index=False)\n",
    "chain4_turned_twisted.to_frame().to_csv(\n",
    "    \"results/50nm_4_turn_twist.csv\", sep=\" \", index=False\n",
    ")\n",
    "\n",
    "chain4_straight.to_plot().savefig(\"results/50nm_4_straight.png\")\n",
    "chain4_turned.to_plot().savefig(\"results/50nm_4_turn.png\")\n",
    "chain4_turned_twisted.to_plot().savefig(\"results/50nm_4_turn_twist.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='results/50nm_4_straight.png' width='30%'> <img src='results/50nm_4_turn.png' width='30%'> <img src='results/50nm_4_turn_twist.png' width='30%'>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "bp_separation = dnachain.BP_SEPARATION  # 3.32Å\n",
    "side_length_nm = 50  # nm\n",
    "num_basepairs_straight = int(side_length_nm / (0.1 * bp_separation))\n",
    "num_basepairs_turned = int((side_length_nm * np.pi / 4.0) / (0.1 * bp_separation))\n",
    "strand_separation_1 = 100  # angstroms\n",
    "strand_separation_2 = 250  # angstroms\n",
    "\n",
    "chain8_straight = dnachain.EightStrandDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_straight)),\n",
    "    strand_separation_1,\n",
    "    strand_separation_2,\n",
    "    turn=False,\n",
    "    twist=False,\n",
    ")\n",
    "\n",
    "chain8_turned = dnachain.EightStrandDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned)),\n",
    "    strand_separation_1,\n",
    "    strand_separation_2,\n",
    "    turn=True,\n",
    "    twist=False,\n",
    ")\n",
    "\n",
    "chain8_turned_twisted = dnachain.EightStrandDNAChain(\n",
    "    \"\".join(np.random.choice([\"G\", \"A\", \"T\", \"C\"], num_basepairs_turned)),\n",
    "    strand_separation_1,\n",
    "    strand_separation_2,\n",
    "    turn=True,\n",
    "    twist=True,\n",
    ")\n",
    "\n",
    "chain8_straight.to_frame().to_csv(\n",
    "    \"results/50nm_8_straight.csv\",\n",
    "    sep=\" \",\n",
    "    index=False,\n",
    ")\n",
    "chain8_turned.to_frame().to_csv(\"results/50nm_8_turn.csv\", sep=\" \", index=False)\n",
    "chain8_turned_twisted.to_frame().to_csv(\n",
    "    \"results/50nm_8_turn_twist.csv\", sep=\" \", index=False\n",
    ")\n",
    "\n",
    "chain8_straight.to_plot().savefig(\"results/50nm_8_straight.png\")\n",
    "chain8_turned.to_plot().savefig(\"results/50nm_8_turn.png\")\n",
    "chain8_turned_twisted.to_plot().savefig(\"results/50nm_8_turn_twist.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='results/50nm_8_straight.png' width='30%'> <img src='results/50nm_8_turn.png' width='30%'> <img src='results/50nm_8_turn_twist.png' width='30%'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Making Solenoidal DNA\n",
    "\n",
    "### Single Solenoids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "side_length = 750  # angstrom\n",
    "radius_solenoid = 100  # angstrom\n",
    "nhistones = 38  # histones\n",
    "\n",
    "solenoid_straight = dnachain.Solenoid(\n",
    "    voxelheight=side_length, radius=radius_solenoid, nhistones=nhistones\n",
    ")\n",
    "solenoid_turned = dnachain.TurnedSolenoid(\n",
    "    voxelheight=side_length, radius=radius_solenoid, nhistones=nhistones\n",
    ")\n",
    "solenoid_turned_twisted = dnachain.TurnedSolenoid(\n",
    "    voxelheight=side_length, radius=radius_solenoid, nhistones=nhistones, twist=True\n",
    ")\n",
    "\n",
    "# centre around (x,y,z)=(0,0,0)\n",
    "solenoid_straight.translate([0, 0, -side_length / 2.0])\n",
    "solenoid_turned.translate([0, 0, -side_length / 2.0])\n",
    "solenoid_turned_twisted.translate([0, 0, -side_length / 2.0])\n",
    "\n",
    "solenoid_straight.to_frame().to_csv(\n",
    "    \"results/solenoid_straight.csv\", sep=\" \", index=False\n",
    ")\n",
    "solenoid_turned.to_frame().to_csv(\"results/solenoid_turned.csv\", sep=\" \", index=False)\n",
    "solenoid_turned_twisted.to_frame().to_csv(\n",
    "    \"results/solenoid_turned_twisted.csv\", sep=\" \", index=False\n",
    ")\n",
    "\n",
    "plot = solenoid_straight.to_line_plot()\n",
    "plot.scene.save_jpg(\"results/solenoid_straight.jpg\")\n",
    "\n",
    "plot = solenoid_turned.to_line_plot()\n",
    "distance = 1500\n",
    "mlab.view(azimuth=180, elevation=0, distance=distance, focalpoint=[0, 0, 0])\n",
    "mlab.move(up=-distance, forward=distance)\n",
    "mlab.pitch(90)\n",
    "plot.scene.save_jpg(\"results/solenoid_turned.jpg\")\n",
    "\n",
    "plot = solenoid_turned_twisted.to_line_plot()\n",
    "mlab.view(azimuth=180, elevation=0, distance=distance, focalpoint=[0, 0, 0])\n",
    "mlab.move(up=-distance, forward=distance)\n",
    "mlab.pitch(90)\n",
    "plot.scene.save_jpg(\"results/solenoid_turned_twisted.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='results/solenoid_straight.jpg' width='30%'> <img src='results/solenoid_turned.jpg' width='30%'> <img src='results/solenoid_turned_twisted.jpg' width='30%'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating multiple solenoids in a Volume\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "side_length = 1000  # angstrom\n",
    "radius_solenoid = 100  # angstrom\n",
    "nhistones = 51  # histones\n",
    "separation = 250  # angstroms\n",
    "\n",
    "solenoid4_straight = dnachain.MultiSolenoidVolume(\n",
    "    voxelheight=side_length,\n",
    "    separation=separation,\n",
    "    radius=radius_solenoid,\n",
    "    nhistones=nhistones,\n",
    "    chains=[1, 2, 3, 4],\n",
    "    turn=False,\n",
    "    twist=False,\n",
    ")\n",
    "\n",
    "solenoid4_turned = dnachain.MultiSolenoidVolume(\n",
    "    voxelheight=side_length,\n",
    "    separation=separation,\n",
    "    radius=radius_solenoid,\n",
    "    nhistones=nhistones,\n",
    "    chains=[1, 2, 3, 4],\n",
    "    turn=True,\n",
    "    twist=False,\n",
    ")\n",
    "\n",
    "solenoid4_turned_twisted = dnachain.MultiSolenoidVolume(\n",
    "    voxelheight=side_length,\n",
    "    separation=separation,\n",
    "    radius=radius_solenoid,\n",
    "    nhistones=nhistones,\n",
    "    chains=[1, 2, 3, 4],\n",
    "    turn=True,\n",
    "    twist=True,\n",
    ")\n",
    "\n",
    "# centre around (x,y,z)=(0,0,0)\n",
    "solenoid4_straight.translate([0, 0, -side_length / 2.0])\n",
    "solenoid4_turned.translate([0, 0, -side_length / 2.0])\n",
    "solenoid4_turned_twisted.translate([0, 0, -side_length / 2.0])\n",
    "\n",
    "solenoid4_straight.to_frame().to_csv(\n",
    "    \"results/solenoid4_straight.csv\", sep=\" \", index=False\n",
    ")\n",
    "solenoid4_turned.to_frame().to_csv(\"results/solenoid4_turned.csv\", sep=\" \", index=False)\n",
    "solenoid4_turned_twisted.to_frame().to_csv(\n",
    "    \"results/solenoid4_turned_twisted.csv\", sep=\" \", index=False\n",
    ")\n",
    "\n",
    "plot = solenoid4_straight.to_line_plot()\n",
    "plot.scene.save_jpg(\"results/solenoid4_straight.jpg\")\n",
    "\n",
    "plot = solenoid4_turned.to_line_plot()\n",
    "distance = 2500\n",
    "mlab.view(azimuth=180, elevation=0, distance=distance, focalpoint=[0, 0, 0])\n",
    "mlab.move(up=-distance, forward=distance)\n",
    "mlab.pitch(90)\n",
    "plot.scene.save_jpg(\"results/solenoid4_turned.jpg\")\n",
    "\n",
    "plot = solenoid4_turned_twisted.to_line_plot()\n",
    "mlab.view(azimuth=180, elevation=0, distance=distance, focalpoint=[0, 0, 0])\n",
    "mlab.move(up=-distance, forward=distance)\n",
    "mlab.pitch(90)\n",
    "plot.scene.save_jpg(\"results/solenoid4_turned_twisted.jpg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='results/solenoid4_straight.jpg' width='40%'> <img src='results/solenoid4_turned.jpg' width='40%'> "
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}