Normal modes animator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------\n",
    "INSTRUCTIONS\n",
    "-------------\n",
    "\n",
    "This code allows obtaining animations of any of the previously calculate normal modes of vibration of a biomolecule. The \"FIRST EXTREMUM\" and \"SECOND EXTREMUM\" cells allow you to obtain a static representation of a specific normal mode's extrema.\n",
    "\n",
    "Run each cell from top to bottom. The final cell allows you to save an .mp4 movie of any animation you may have generated. Running the final cell is not mandatory, only do it if you generated movies that you want to save it. The movie will be saved to the folder where this code is stored.\n",
    "\n",
    "If this is your first time running the code, type ctrl + F and then type r'C, this will take you to the parts of the code where an instruction to read a specific path in your computer is given. Make sure that the path exists in your computer and that it corresponds to the path where the file being specified can be found.\n",
    "\n",
    "------------\n",
    "DESCRIPTIONS\n",
    "-------------\n",
    "\n",
    "**Coordinates file:** File containing the coordinates of the biomolecule which normal modes you wish to animate.\n",
    "\n",
    "**Normal modes file:** File containing the normal modes of the biomolecule which normal modes you wish to animate.\n",
    "\n",
    "**Available number of modes:** This output tells you how many normal modes can be animated.\n",
    "\n",
    "**What mode do you want to animate?:** Choose the mode you wish to animate. You cannot animate the translational or rotational modes (modes 1 to 6) or a mode above the available number of modes.\n",
    "\n",
    "**Graph scaling factor:** This is related to the size of the graph. A large number will result in a very small representation of the molecule, except when all coordinates of the molecule are positive numbers. Suggested value: 1.5.\n",
    "\n",
    "**Maximum normalized value:** This is related to the arbitrary amplitude of motion in the animation. A small maximum normalized value will result in hardly perceivable motion, while a large maximum normalized value will result in clear motion. Suggested value: 100.\n",
    "\n",
    "**Show axes?:** Type \"yes\" if you want to see the cartesian axes, otherwise type \"no\".\n",
    "\n",
    "**Elevation (degrees):** Graph angle. Suggested value: 0.\n",
    "\n",
    "**Azimuth (degrees):** Graph angle. Suggested value: 0.\n",
    "\n",
    "**Node size:** How big do you want the nodes to be? Suggested value: 20 to 50 for small biomolecules (such as a single protein), 0.5 or 1 for bigger biomolecules (such as an entire viral capsid)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LIBRARIES AND FUNCTIONS\n",
    "\n",
    "# Libraries\n",
    "print('\\033[37m' + 'Importing libraries...'  + '\\033[0m')\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation\n",
    "from pathlib import Path\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "print('\\033[37m' + 'Defining functions...'  + '\\033[0m')\n",
    "\n",
    "# Integer validation\n",
    "def is_number(s):\n",
    "    try:\n",
    "        float(s)\n",
    "        return True\n",
    "    except ValueError:\n",
    "        pass\n",
    "    try:\n",
    "        import unicodedata\n",
    "        unicodedata.numeric(s)\n",
    "        return True\n",
    "    except (TypeError, ValueError):\n",
    "        pass\n",
    "    return False\n",
    "\n",
    "# String substraction\n",
    "def substract(a, b):                              \n",
    "    return \"\".join(a.rsplit(b))\n",
    "\n",
    "print('\\033[1m' + 'Done' + '\\033[0m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FILES\n",
    "\n",
    "# Coordinates file directory\n",
    "CoordinatesFile = 0\n",
    "while CoordinatesFile == 0:\n",
    "    CoordinatesFile = input('\\033[1m' + 'Coordinates file: ' + '\\033[0')\n",
    "    if 'coordinates' not in CoordinatesFile:\n",
    "        CoordinatesFile = 0\n",
    "        print('This is not a coordinates file')\n",
    "    elif 'coordinates' in CoordinatesFile:\n",
    "        CoordinatesFile = CoordinatesFile + '.csv'\n",
    "        CoordinatesPath  = os.path.join(r'C:\\Users\\santi\\Desktop\\2019\\Biophysics\\Data', CoordinatesFile)\n",
    "        if Path(CoordinatesPath).is_file():\n",
    "            MoleculeName = substract(CoordinatesFile, '_coordinates.csv')\n",
    "        else:\n",
    "            CoordinatesFile = 0\n",
    "            print('File not found')\n",
    "print('\\033[37m' + 'Retrieving coordinates...'  + '\\033[0m')\n",
    "df1 = pd.read_csv(CoordinatesPath)\n",
    "df1.index = np.arange(1, len(df1) + 1)\n",
    "df1.index.name = 'Node'\n",
    "del df1['Unnamed: 0']\n",
    "N = int(df1.shape[0]) # Number of atoms from coordinates file\n",
    "# Equilibrium positions\n",
    "x0 = np.zeros((N))\n",
    "y0 = np.zeros((N))\n",
    "z0 = np.zeros((N))\n",
    "for i in range(0, N):\n",
    "    x0[i] = df1['x0'][i+1]\n",
    "    y0[i] = df1['y0'][i+1]\n",
    "    z0[i] = df1['z0'][i+1]\n",
    "            \n",
    "# Normal modes file directory\n",
    "NormalModesFile = 0\n",
    "while NormalModesFile == 0:\n",
    "    NormalModesFile = input('\\033[1m' + 'Normal modes file: ' + '\\033[0')\n",
    "    if 'normalModes' not in NormalModesFile:\n",
    "        NormalModesFile = 0\n",
    "        print('This is not a normal modes file')\n",
    "    elif 'normalModes' in NormalModesFile:\n",
    "        NormalModesFile = NormalModesFile + '.csv'\n",
    "        NormalModesPath  = os.path.join(r'C:\\Users\\santi\\Desktop\\2019\\Biophysics\\Data', NormalModesFile)\n",
    "        if Path(NormalModesPath).is_file() == False:\n",
    "            NormalModesFile = 0\n",
    "            print('File not found')\n",
    "print('\\033[37m' + 'Retrieving normal modes...'  + '\\033[0m')\n",
    "df2 = pd.read_csv(NormalModesPath)\n",
    "df2.index = np.arange(1, len(df2) + 1)\n",
    "df2.index = ['x' + str(i) for i in range(1, N + 1)] + ['y' + str(i) for i in range(1, N + 1)] + ['z' + str(i) for i in range(1, N + 1)]\n",
    "df2.index.name = 'Coordinate'\n",
    "del df2['Coordinate']\n",
    "HowMany = df2.shape[1] # Number of modes from normal modes file\n",
    "Eigenvectors = np.zeros((3*N, HowMany))\n",
    "Eigenvectors = df2.values # Eigenmatrix from normal modes file\n",
    "\n",
    "print('\\033[1m' + 'Done. ' + '\\033[0m' + 'There are ' + '\\033[1m' + str(HowMany) + '\\033[0m' + ' modes available (from the lowest)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NORMAL MODE ANIMATION\n",
    "    \n",
    "# Enable dynamic graphs mode\n",
    "%matplotlib nbagg\n",
    "    \n",
    "# Normal mode to animate\n",
    "Mode = 0\n",
    "while Mode == 0:\n",
    "    Mode = input('\\033[1m' + 'What mode do you want to animate? ' + '\\033[0')\n",
    "    if is_number(Mode) == True:\n",
    "        if float(Mode) <= 0:\n",
    "            Mode = 0\n",
    "            print('Must be between 1 and ' + str(HowMany))\n",
    "        elif float(Mode) > 3*N:\n",
    "            Mode = 0\n",
    "            print('Must be between 1 and ' + str(HowMany))\n",
    "        else:\n",
    "            if (float(Mode)).is_integer() == False:\n",
    "                Mode = 0\n",
    "                print('Must be an integer')\n",
    "            else:\n",
    "                Mode = int(Mode)\n",
    "    else:\n",
    "        Mode = 0\n",
    "        print('Must be a number')\n",
    "Mode = Mode - 1\n",
    "\n",
    "# Axes scaling factor\n",
    "ScalingFactor = 0\n",
    "while ScalingFactor == 0:\n",
    "    ScalingFactor = input('\\033[1m' + 'Graph scaling factor: ' + '\\033[0')\n",
    "    if is_number(ScalingFactor) == True:\n",
    "        if float(ScalingFactor) < 1:\n",
    "            ScalingFactor = 0\n",
    "            print('Scaling factor cannot be less than 1')\n",
    "        else:\n",
    "            ScalingFactor = float(ScalingFactor)\n",
    "    else:\n",
    "        ScalingFactor = 0\n",
    "        print('Must be a number')\n",
    "\n",
    "# Maximum normalized value\n",
    "MaxNormalizedValue = 0\n",
    "while MaxNormalizedValue == 0:\n",
    "    MaxNormalizedValue = input('\\033[1m' + 'Maximum normalized value: ' + '\\033[0')\n",
    "    if is_number(MaxNormalizedValue) == True:\n",
    "        if float(MaxNormalizedValue) < 0:\n",
    "            MaxNormalizedValue = 0\n",
    "            print('Must be positive')\n",
    "        elif float(MaxNormalizedValue) == 0:\n",
    "            print('Maximum normalized value cannot be zero')\n",
    "        else:\n",
    "            MaxNormalizedValue = float(MaxNormalizedValue)\n",
    "    else:\n",
    "        MaxNormalizedValue = 0\n",
    "        print('Must be a number')\n",
    "        \n",
    "# Show axes?\n",
    "ShowAxes = 0\n",
    "while ShowAxes == 0:\n",
    "    ShowAxes = input('\\033[1m' + 'Show axes?: ' + '\\033[0')\n",
    "    if ShowAxes not in ['yes', 'no']:\n",
    "        ShowAxes = 0\n",
    "        print('Answer must be yes or no')\n",
    "\n",
    "# Initial angles\n",
    "Angle1 = -1\n",
    "Angle2 = -1\n",
    "while Angle1 == -1:\n",
    "    Angle1 = input('\\033[1m' + 'Elevation (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle1) == True:\n",
    "        if float(Angle1) < 0:\n",
    "            Angle1 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle1) > 360:\n",
    "            Angle1 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle1 = float(Angle1)\n",
    "    else:\n",
    "        Angle1 = -1\n",
    "        print('Must be a number')\n",
    "while Angle2 == -1:\n",
    "    Angle2 = input('\\033[1m' + 'Azimuth (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle2) == True:\n",
    "        if float(Angle2) < 0:\n",
    "            Angle2 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle2) > 360:\n",
    "            Angle2 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle2 = float(Angle2)\n",
    "    else:\n",
    "        Angle2 = -1\n",
    "        print('Must be a number')\n",
    "\n",
    "# Marker size\n",
    "MarkerSize = 0\n",
    "while MarkerSize == -0:\n",
    "    MarkerSize = input('\\033[1m' + 'Node size: ' + '\\033[0')\n",
    "    if is_number(MarkerSize) == True:\n",
    "        if float(MarkerSize) < 0:\n",
    "            MarkerSize = 0\n",
    "            print('Must be positive')\n",
    "        elif float(MarkerSize) == 0:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be zero')\n",
    "        elif float(MarkerSize) > 50:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be greater than 50')\n",
    "        else:\n",
    "            MarkerSize = float(MarkerSize)\n",
    "    else:\n",
    "        MarkerSize = 0\n",
    "        print('Must be a number')\n",
    "        \n",
    "# Normalized eigenmatrix\n",
    "print('\\033[37m' + 'Defining normalized eigenmatrix...'  + '\\033[0m')\n",
    "NormalizedEigenvectors = np.zeros((3*N, HowMany))\n",
    "Max = np.amax(Eigenvectors)\n",
    "for i in range(0, 3*N):\n",
    "    for j in range(0, HowMany):\n",
    "        NormalizedEigenvectors[i,j] = Eigenvectors[i,j]*(MaxNormalizedValue/Max)\n",
    "df3 = pd.DataFrame(NormalizedEigenvectors)\n",
    "df3.index = ['x' + str(i) for i in range(1, N + 1)] + ['y' + str(i) for i in range(1, N + 1)] + ['z' + str(i) for i in range(1, N + 1)]\n",
    "df3.index.name = 'Coordinate'\n",
    "df3.columns = np.arange(1, HowMany + 1)\n",
    "df3.columns.name = 'Mode number'\n",
    "\n",
    "print('\\033[37m' + 'Defining animation...'  + '\\033[0m')\n",
    "# Auxiliary arrays to avoid repetition of transpose function\n",
    "AuxiliaryX = NormalizedEigenvectors[0:N, Mode] # First mode's normalized relative amplitudes of motion for x\n",
    "AuxiliaryY = NormalizedEigenvectors[N:2*N, Mode] # First mode's normalized relative amplitudes of motion for y\n",
    "AuxiliaryZ = NormalizedEigenvectors[2*N:3*N, Mode] # First mode's normalized relative amplitudes of motion for z\n",
    "# Time array that counts from 0 to N, 5 times each (molecule goes from origin to extreme to origin to extreme to origin)\n",
    "t = np.array([np.ones(N)*i for i in range(8)]).flatten()\n",
    "# Vectors containing all the modified coordinates\n",
    "x0Motion = np.zeros((8*N))\n",
    "y0Motion = np.zeros((8*N))\n",
    "z0Motion = np.zeros((8*N))\n",
    "for i in range(0, N):\n",
    "    x0Motion[i] = x0[i]\n",
    "    y0Motion[i] = y0[i]\n",
    "    z0Motion[i] = z0[i]\n",
    "for i in range(N, 2*N):\n",
    "    x0Motion[i] = x0[i-N] + (AuxiliaryX[i-N])/2\n",
    "    y0Motion[i] = y0[i-N] + (AuxiliaryY[i-N])/2\n",
    "    z0Motion[i] = z0[i-N] + (AuxiliaryZ[i-N])/2\n",
    "for i in range(2*N, 3*N):\n",
    "    x0Motion[i] = x0[i-2*N] + AuxiliaryX[i-2*N]\n",
    "    y0Motion[i] = y0[i-2*N] + AuxiliaryY[i-2*N]\n",
    "    z0Motion[i] = z0[i-2*N] + AuxiliaryZ[i-2*N]\n",
    "for i in range(3*N, 4*N):\n",
    "    x0Motion[i] = x0[i-3*N] + (AuxiliaryX[i-3*N])/2\n",
    "    y0Motion[i] = y0[i-3*N] + (AuxiliaryY[i-3*N])/2\n",
    "    z0Motion[i] = z0[i-3*N] + (AuxiliaryZ[i-3*N])/2\n",
    "for i in range(4*N, 5*N):\n",
    "    x0Motion[i] = x0[i-4*N]\n",
    "    y0Motion[i] = y0[i-4*N]\n",
    "    z0Motion[i] = z0[i-4*N]\n",
    "for i in range(5*N, 6*N):\n",
    "    x0Motion[i] = x0[i-5*N] - (AuxiliaryX[i-5*N])/2\n",
    "    y0Motion[i] = y0[i-5*N] - (AuxiliaryY[i-5*N])/2\n",
    "    z0Motion[i] = z0[i-5*N] - (AuxiliaryZ[i-5*N])/2\n",
    "for i in range(6*N, 7*N):\n",
    "    x0Motion[i] = x0[i-6*N] - AuxiliaryX[i-6*N]\n",
    "    y0Motion[i] = y0[i-6*N] - AuxiliaryY[i-6*N]\n",
    "    z0Motion[i] = z0[i-6*N] - AuxiliaryZ[i-6*N]\n",
    "for i in range(7*N, 8*N):\n",
    "    x0Motion[i] = x0[i-7*N] - (AuxiliaryX[i-7*N])/2\n",
    "    y0Motion[i] = y0[i-7*N] - (AuxiliaryY[i-7*N])/2\n",
    "    z0Motion[i] = z0[i-7*N] - (AuxiliaryZ[i-7*N])/2\n",
    "df4 = pd.DataFrame({'time': t ,'x' : x0Motion, 'y' : y0Motion, 'z' : z0Motion}) # Dataframe containing all coordinates\n",
    "# FuncAnimation function argument\n",
    "def update_graph(num):\n",
    "    data = df4[df4['time'] == num]\n",
    "    graph.set_data (data.x, data.y)\n",
    "    graph.set_3d_properties(data.z)\n",
    "    return graph, \n",
    "# Build graph object\n",
    "fig = plt.figure(figsize = (9, 9))\n",
    "ax = fig.add_subplot(111, projection = '3d')\n",
    "ax.view_init(Angle1, Angle2) # Initial angle\n",
    "# x-axis parameters\n",
    "ax.set_xlabel('x (angstroms)')\n",
    "xmin = np.amin(x0)\n",
    "xmax = np.amax(x0)\n",
    "if (xmin < 0 and xmax < 0) or (xmin < 0 and xmax > 0):\n",
    "    ax.set_xlim(ScalingFactor*xmin, xmax+xmin*(1-ScalingFactor))\n",
    "if xmin > 0 and xmax > 0:\n",
    "    ax.set_xlim(xmin*(1-1/ScalingFactor), xmax+(1/ScalingFactor)*xmin)\n",
    "# y-axis parameters\n",
    "ax.set_ylabel('y (angstroms)')\n",
    "ymin = np.amin(y0)\n",
    "ymax = np.amax(y0)\n",
    "if (ymin < 0 and ymax < 0) or (ymin < 0 and ymax > 0):\n",
    "    ax.set_ylim(ScalingFactor*ymin, ymax+ymin*(1-ScalingFactor))\n",
    "if ymin > 0 and ymax > 0:\n",
    "    ax.set_ylim(ymin*(1-1/ScalingFactor), ymax+(1/ScalingFactor)*ymin)\n",
    "# z-axis parameters\n",
    "ax.set_zlabel('z (angstroms)')\n",
    "zmin = np.amin(z0)\n",
    "zmax = np.amax(z0)\n",
    "if (zmin < 0 and zmax < 0) or (zmin < 0 and zmax > 0):\n",
    "    ax.set_zlim(ScalingFactor*zmin, zmax+zmin*(1-ScalingFactor))\n",
    "if zmin > 0 and zmax > 0:\n",
    "    ax.set_zlim(zmin*(1-1/ScalingFactor), zmax+(1/ScalingFactor)*zmin)\n",
    "# Show or hide axes\n",
    "if ShowAxes == 'no':\n",
    "    ax.set_axis_off()\n",
    "elif ShowAxes == 'yes':\n",
    "    ax.set_axis_on()\n",
    "data = df4[df4['time'] == 0] # Initial dataset for the animated graph\n",
    "graph, = ax.plot(data.x, data.y, data.z, linestyle = '', marker = 'o', color = 'black', markersize = MarkerSize) # Graph's main parameters\n",
    "animation = matplotlib.animation.FuncAnimation(fig, update_graph, frames = 8, blit = True, interval = 200) # Generate animation\n",
    "plt.show() # Show graph\n",
    "\n",
    "print('\\033[1m' + 'Done' + '\\033[0m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIRST EXTREMUM\n",
    "\n",
    "# Axes scaling factor\n",
    "ScalingFactor = 0\n",
    "while ScalingFactor == 0:\n",
    "    ScalingFactor = input('\\033[1m' + 'Graph scaling factor: ' + '\\033[0')\n",
    "    if is_number(ScalingFactor) == True:\n",
    "        if float(ScalingFactor) < 0:\n",
    "            ScalingFactor = 0\n",
    "            print('Must be positive')\n",
    "        elif float(ScalingFactor) == 0:\n",
    "            ScalingFactor = 0\n",
    "            print('Scaling factor cannot be zero')\n",
    "        elif float(ScalingFactor) > 5:\n",
    "            ScalingFactor = 0\n",
    "            print('Scaling factor cannot be greater than 5')\n",
    "        else:\n",
    "            ScalingFactor = float(ScalingFactor)\n",
    "    else:\n",
    "        ScalingFactor = 0\n",
    "        print('Must be a number')\n",
    "\n",
    "# Show axes?\n",
    "ShowAxes = 0\n",
    "while ShowAxes == 0:\n",
    "    ShowAxes = input('\\033[1m' + 'Show axes? ' + '\\033[0')\n",
    "    if ShowAxes not in ['yes', 'no']:\n",
    "        ShowAxes = 0\n",
    "        print('Answer must be yes or no')\n",
    "\n",
    "# Initial angles\n",
    "Angle1 = -1\n",
    "Angle2 = -1\n",
    "while Angle1 == -1:\n",
    "    Angle1 = input('\\033[1m' + 'Angle 1 (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle1) == True:\n",
    "        if float(Angle1) < 0:\n",
    "            Angle1 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle1) > 360:\n",
    "            Angle1 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle1 = float(Angle1)\n",
    "    else:\n",
    "        Angle1 = -1\n",
    "        print('Must be a number')\n",
    "while Angle2 == -1:\n",
    "    Angle2 = input('\\033[1m' + 'Angle 2 (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle2) == True:\n",
    "        if float(Angle2) < 0:\n",
    "            Angle2 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle2) > 360:\n",
    "            Angle2 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle2 = float(Angle2)\n",
    "    else:\n",
    "        Angle2 = -1\n",
    "        print('Must be a number')\n",
    "\n",
    "# Marker size\n",
    "MarkerSize = 0\n",
    "while MarkerSize == -0:\n",
    "    MarkerSize = input('\\033[1m' + 'Node size: ' + '\\033[0')\n",
    "    if is_number(MarkerSize) == True:\n",
    "        if float(MarkerSize) < 0:\n",
    "            MarkerSize = 0\n",
    "            print('Must be positive')\n",
    "        elif float(MarkerSize) == 0:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be zero')\n",
    "        elif float(MarkerSize) > 50:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be greater than 50')\n",
    "        else:\n",
    "            MarkerSize = float(MarkerSize)\n",
    "    else:\n",
    "        MarkerSize = 0\n",
    "        print('Must be a number')\n",
    "\n",
    "# First extremum graph\n",
    "print('\\033[37m' + 'Graphing first extremum of mode ' + str(Mode) + '...'  + '\\033[0m')\n",
    "fig2 = plt.figure(figsize = (9, 9)) # Build figure object\n",
    "ax2 = fig2.add_subplot(111, projection = '3d')\n",
    "ax2.scatter(x0Motion[2*N:3*N], y0Motion[2*N:3*N], z0Motion[2*N:3*N], c = 'black', s = MarkerSize, marker = 'o')\n",
    "ax2.view_init(Angle1, Angle2) # Initial angle\n",
    "# x-axis parameters\n",
    "ax2.set_xlabel('x (angstroms)')\n",
    "xmin = np.amin(x0)\n",
    "xmax = np.amax(x0)\n",
    "if (xmin < 0 and xmax < 0) or (xmin < 0 and xmax > 0):\n",
    "    ax2.set_xlim(ScalingFactor*xmin, xmax+xmin*(1-ScalingFactor))\n",
    "if xmin > 0 and xmax > 0:\n",
    "    ax2.set_xlim(xmin*(1-1/ScalingFactor), xmax+(1/ScalingFactor)*xmin)\n",
    "# y-axis parameters\n",
    "ax2.set_ylabel('y (angstroms)')\n",
    "ymin = np.amin(y0)\n",
    "ymax = np.amax(y0)\n",
    "if (ymin < 0 and ymax < 0) or (ymin < 0 and ymax > 0):\n",
    "    ax2.set_ylim(ScalingFactor*ymin, ymax+ymin*(1-ScalingFactor))\n",
    "if ymin > 0 and ymax > 0:\n",
    "    ax2.set_ylim(ymin*(1-1/ScalingFactor), ymax+(1/ScalingFactor)*ymin)\n",
    "# z-axis parameters\n",
    "ax2.set_zlabel('z (angstroms)')\n",
    "zmin = np.amin(z0)\n",
    "zmax = np.amax(z0)\n",
    "if (zmin < 0 and zmax < 0) or (zmin < 0 and zmax > 0):\n",
    "    ax2.set_zlim(ScalingFactor*zmin, zmax+zmin*(1-ScalingFactor))\n",
    "if zmin > 0 and zmax > 0:\n",
    "    ax2.set_zlim(zmin*(1-1/ScalingFactor), zmax+(1/ScalingFactor)*zmin)\n",
    "# Show or hide axes\n",
    "if ShowAxes == 'no':\n",
    "    ax2.set_axis_off()\n",
    "elif ShowAxes == 'yes':\n",
    "    ax2.set_axis_on()\n",
    "plt.show() # Show graph\n",
    "\n",
    "print('\\033[1m' + 'Done' + '\\033[0m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SECOND EXTREMUM\n",
    "\n",
    "# Axes scaling factor\n",
    "ScalingFactor = 0\n",
    "while ScalingFactor == 0:\n",
    "    ScalingFactor = input('\\033[1m' + 'Graph scaling factor: ' + '\\033[0')\n",
    "    if is_number(ScalingFactor) == True:\n",
    "        if float(ScalingFactor) < 0:\n",
    "            ScalingFactor = 0\n",
    "            print('Must be positive')\n",
    "        elif float(ScalingFactor) == 0:\n",
    "            ScalingFactor = 0\n",
    "            print('Scaling factor cannot be zero')\n",
    "        elif float(ScalingFactor) > 5:\n",
    "            ScalingFactor = 0\n",
    "            print('Scaling factor cannot be greater than 5')\n",
    "        else:\n",
    "            ScalingFactor = float(ScalingFactor)\n",
    "    else:\n",
    "        ScalingFactor = 0\n",
    "        print('Must be a number')\n",
    "\n",
    "# Show axes?\n",
    "ShowAxes = 0\n",
    "while ShowAxes == 0:\n",
    "    ShowAxes = input('\\033[1m' + 'Show axes? ' + '\\033[0')\n",
    "    if ShowAxes not in ['yes', 'no']:\n",
    "        ShowAxes = 0\n",
    "        print('Answer must be yes or no')\n",
    "\n",
    "# Initial angles\n",
    "Angle1 = -1\n",
    "Angle2 = -1\n",
    "while Angle1 == -1:\n",
    "    Angle1 = input('\\033[1m' + 'Angle 1 (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle1) == True:\n",
    "        if float(Angle1) < 0:\n",
    "            Angle1 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle1) > 360:\n",
    "            Angle1 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle1 = float(Angle1)\n",
    "    else:\n",
    "        Angle1 = -1\n",
    "        print('Must be a number')\n",
    "while Angle2 == -1:\n",
    "    Angle2 = input('\\033[1m' + 'Angle 2 (degrees): ' + '\\033[0')\n",
    "    if is_number(Angle2) == True:\n",
    "        if float(Angle2) < 0:\n",
    "            Angle2 = -1\n",
    "            print('Must be positive')\n",
    "        elif float(Angle2) > 360:\n",
    "            Angle2 = -1\n",
    "            print('Cannot be greater than 360°')\n",
    "        else:\n",
    "            Angle2 = float(Angle2)\n",
    "    else:\n",
    "        Angle2 = -1\n",
    "        print('Must be a number')  \n",
    "\n",
    "# Marker size\n",
    "MarkerSize = 0\n",
    "while MarkerSize == -0:\n",
    "    MarkerSize = input('\\033[1m' + 'Node size: ' + '\\033[0')\n",
    "    if is_number(MarkerSize) == True:\n",
    "        if float(MarkerSize) < 0:\n",
    "            MarkerSize = 0\n",
    "            print('Must be positive')\n",
    "        elif float(MarkerSize) == 0:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be zero')\n",
    "        elif float(MarkerSize) > 50:\n",
    "            MarkerSize = 0\n",
    "            print('Cannot be greater than 50')\n",
    "        else:\n",
    "            MarkerSize = float(MarkerSize)\n",
    "    else:\n",
    "        MarkerSize = 0\n",
    "        print('Must be a number')\n",
    "\n",
    "# Second extremum graph\n",
    "print('\\033[37m' + 'Graphing second extremum of mode ' + str(Mode) + '...'  + '\\033[0m')\n",
    "fig3 = plt.figure(figsize = (9, 9)) # Build figure object\n",
    "ax3 = fig3.add_subplot(111, projection = '3d')\n",
    "ax3.scatter(x0Motion[6*N:7*N], y0Motion[6*N:7*N], z0Motion[6*N:7*N], c = 'black', s = MarkerSize, marker = 'o')\n",
    "ax3.view_init(Angle1, Angle2) # Initial angle\n",
    "# x-axis parameters\n",
    "ax3.set_xlabel('x (angstroms)')\n",
    "xmin = np.amin(x0)\n",
    "xmax = np.amax(x0)\n",
    "if (xmin < 0 and xmax < 0) or (xmin < 0 and xmax > 0):\n",
    "    ax3.set_xlim(ScalingFactor*xmin, xmax+xmin*(1-ScalingFactor))\n",
    "if xmin > 0 and xmax > 0:\n",
    "    ax3.set_xlim(xmin*(1-1/ScalingFactor), xmax+(1/ScalingFactor)*xmin)\n",
    "# y-axis parameters\n",
    "ax3.set_ylabel('y (angstroms)')\n",
    "ymin = np.amin(y0)\n",
    "ymax = np.amax(y0)\n",
    "if (ymin < 0 and ymax < 0) or (ymin < 0 and ymax > 0):\n",
    "    ax3.set_ylim(ScalingFactor*ymin, ymax+ymin*(1-ScalingFactor))\n",
    "if ymin > 0 and ymax > 0:\n",
    "    ax3.set_ylim(ymin*(1-1/ScalingFactor), ymax+(1/ScalingFactor)*ymin)\n",
    "# z-axis parameters\n",
    "ax3.set_zlabel('z (angstroms)')\n",
    "zmin = np.amin(z0)\n",
    "zmax = np.amax(z0)\n",
    "if (zmin < 0 and zmax < 0) or (zmin < 0 and zmax > 0):\n",
    "    ax3.set_zlim(ScalingFactor*zmin, zmax+zmin*(1-ScalingFactor))\n",
    "if zmin > 0 and zmax > 0:\n",
    "    ax3.set_zlim(zmin*(1-1/ScalingFactor), zmax+(1/ScalingFactor)*zmin)\n",
    "# Show or hide axes\n",
    "if ShowAxes == 'no':\n",
    "    ax3.set_axis_off()\n",
    "elif ShowAxes == 'yes':\n",
    "    ax3.set_axis_on()\n",
    "plt.show() # Show graph\n",
    "\n",
    "print('\\033[1m' + 'Done' + '\\033[0m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SAVE ANIMATION\n",
    "\n",
    "print('\\033[37m' + 'Saving animation...'  + '\\033[0m')\n",
    "animation.save(MoleculeName + '_mode' + str(Mode + 1) + '_scaled' + str(MaxNormalizedValue) + '.mp4', writer = 'ffmpeg', fps = 10)\n",
    "print('\\033[1m' + 'Done' + '\\033[0m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}