{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:28.675365Z",
     "start_time": "2019-10-08T13:04:24.575595Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline \n",
    "import json\n",
    "import numpy as np\n",
    "import operator\n",
    "import os\n",
    "import pandas\n",
    "import pylab as plt\n",
    "import random\n",
    "from ase.data import reference_states, atomic_numbers\n",
    "from sklearn.neighbors import KernelDensity\n",
    "from pyiron.project import Project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "mpl.rc('font', family='Times New Roman')\n",
    "plt.rcParams[\"mathtext.fontset\"] = \"stix\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyiron.thermodynamics.interfacemethod import freeze_one_half, remove_selective_dynamics, set_server, create_job_template, fix_iso, fix_z_dir, half_velocity, minimize_pos, minimize_vol, next_calc, npt_solid, npt_liquid, next_step_funct, round_temperature_next, strain_circle, analyse_minimized_structure, get_press, get_center_point, get_strain_lst, get_nve_job_name, plot_solid_liquid_ratio, ratio_selection, plot_equilibration, plot_melting_point_prediction, calc_temp_iteration, get_initial_melting_temperature_guess, validate_convergence, initialise_iterators, get_voronoi_volume, check_for_holes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:29.544782Z",
     "start_time": "2019-10-08T13:04:28.677745Z"
    }
   },
   "outputs": [],
   "source": [
    "pr = Project('melting')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "input_file = os.path.join(pr.path, 'input.json')\n",
    "output_file = os.path.join(pr.path, 'output.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:30.547878Z",
     "start_time": "2019-10-08T13:04:30.051880Z"
    }
   },
   "outputs": [],
   "source": [
    "pr.job_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.exists(input_file) and os.stat(input_file).st_size != 0:\n",
    "    with open(input_file, 'r') as f:\n",
    "        input_dict = json.load(f)\n",
    "else:\n",
    "    input_dict = {\n",
    "        \"config\": [\n",
    "            \"pair_style morse 9.97749\\n\",\n",
    "            \"pair_coeff * * 0.4174 1.3885 2.845\\n\"\n",
    "        ],\n",
    "        \"species\": [\"Al\"],\n",
    "        \"element\": \"Al\"\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot_dict = input_dict.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'model' not in pot_dict.keys():\n",
    "    pot_dict['model'] = 'Lammps'\n",
    "if 'name' not in pot_dict.keys():\n",
    "    pot_dict['name'] = 'CustomPotential'\n",
    "if 'filename' not in pot_dict.keys():\n",
    "    pot_path = []\n",
    "else:\n",
    "    pot_path = [os.path.abspath(pot_dict['filename'])]\n",
    "potential = pandas.DataFrame({'Config': [pot_dict['config']],\n",
    "                             'Filename': [pot_path],\n",
    "                             'Model': [pot_dict['model']],\n",
    "                             'Name': [pot_dict['name']],\n",
    "                             'Species': [pot_dict['species']]\n",
    "                             })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:30.557204Z",
     "start_time": "2019-10-08T13:04:30.549236Z"
    }
   },
   "outputs": [],
   "source": [
    "project_parameter = {\n",
    "    'project': pr,\n",
    "    'run_time_steps': 50000,\n",
    "    'nvt_run_time_steps': 10000,\n",
    "    'nve_run_time_steps': 10000,\n",
    "    'temperature_left': 0,\n",
    "    'temperature_right': 1000,\n",
    "    'strain_run_time_steps': 1000,\n",
    "    'convergence_criterion': 1,\n",
    "    'potential': potential,\n",
    "    'cpu_cores': 1, \n",
    "    'job_type': pr.job_type.Lammps,\n",
    "    'enable_h5md': False,\n",
    "    'points': 21,\n",
    "    'boundary_value': 0.25,\n",
    "    'ratio_boundary': 0.25,\n",
    "    'timestep_lst': [2, 2, 1],\n",
    "    'fit_range_lst': [0.05, 0.01, 0.01],\n",
    "    'nve_run_time_steps_lst': [25000, 20000, 50000],\n",
    "    'number_of_atoms': 8000, \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in input_dict.keys():\n",
    "    project_parameter[k] = input_dict[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'crystalstructure' not in project_parameter.keys():\n",
    "    project_parameter['crystalstructure'] = reference_states[atomic_numbers[project_parameter['element']]]['symmetry']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'seed' not in project_parameter.keys():\n",
    "    project_parameter['seed'] = random.randint(0,99999)\n",
    "project_parameter['seed']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:30.580662Z",
     "start_time": "2019-10-08T13:04:30.571110Z"
    }
   },
   "outputs": [],
   "source": [
    "# Values from a previous calculation can be inserted here to reproduce the results \n",
    "step_dict = {}\n",
    "if os.path.exists(output_file):\n",
    "    with open(output_file, 'r') as f:\n",
    "        step_dict_str = json.load(f)\n",
    "    for k,v in step_dict_str.items():\n",
    "        step_dict[int(k)] = v\n",
    "step_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From here on the notebook is automated - no change required !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.021413Z",
     "start_time": "2019-10-08T13:04:31.012101Z"
    }
   },
   "outputs": [],
   "source": [
    "step_count = 0\n",
    "temperature_next = None\n",
    "enable_iteration = True\n",
    "convergence_goal_achieved = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.033447Z",
     "start_time": "2019-10-08T13:04:31.022716Z"
    }
   },
   "outputs": [],
   "source": [
    "pr = project_parameter['project']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.288827Z",
     "start_time": "2019-10-08T13:04:31.034782Z"
    }
   },
   "outputs": [],
   "source": [
    "if 'lattice_constant' in project_parameter.keys():\n",
    "    a = project_parameter['lattice_constant']\n",
    "else:\n",
    "    a = None\n",
    "if project_parameter['crystalstructure'] == 'hcp':\n",
    "    basis = pr.create_ase_bulk(name=project_parameter['element'], crystalstructure=project_parameter['crystalstructure'].lower(), a=a, orthorhombic=True)\n",
    "else:\n",
    "    basis = pr.create_ase_bulk(name=project_parameter['element'], crystalstructure=project_parameter['crystalstructure'].lower(), a=a, cubic=True)\n",
    "basis_lst = [basis.repeat([i, i, i]) for i in range(5,30)]\n",
    "basis = basis_lst[np.argmin([np.abs(len(b)-project_parameter['number_of_atoms']/2) for b in basis_lst])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis.analyse_ovito_cna_adaptive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.515297Z",
     "start_time": "2019-10-08T13:04:31.290502Z"
    }
   },
   "outputs": [],
   "source": [
    "basis.plot3d()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.518804Z",
     "start_time": "2019-10-08T13:04:31.516776Z"
    }
   },
   "outputs": [],
   "source": [
    "# pr.remove_jobs(recursive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:31.531151Z",
     "start_time": "2019-10-08T13:04:31.520073Z"
    }
   },
   "outputs": [],
   "source": [
    "timestep_iter, fit_range_iter, nve_run_time_steps_iter = initialise_iterators(project_parameter)\n",
    "timestep = next(timestep_iter)\n",
    "fit_range = next(fit_range_iter)\n",
    "nve_run_time_steps = next(nve_run_time_steps_iter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:04:32.021404Z",
     "start_time": "2019-10-08T13:04:31.532487Z"
    }
   },
   "outputs": [],
   "source": [
    "pr.job_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 1: set up the solid sample and roughly estimate a melting point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:06:12.707863Z",
     "start_time": "2019-10-08T13:04:32.022971Z"
    }
   },
   "outputs": [],
   "source": [
    "ham_minimize_pos= minimize_pos(\n",
    "    structure=basis, \n",
    "    project_parameter=project_parameter\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:07:54.491128Z",
     "start_time": "2019-10-08T13:06:13.279478Z"
    }
   },
   "outputs": [],
   "source": [
    "ham_minimize_vol = minimize_vol(\n",
    "    structure=ham_minimize_pos.get_structure(), \n",
    "    project_parameter=project_parameter\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:21:25.585304Z",
     "start_time": "2019-10-08T13:07:54.595319Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_next, structure_left = get_initial_melting_temperature_guess(\n",
    "    project_parameter=project_parameter, \n",
    "    ham_minimize_vol=ham_minimize_vol, \n",
    "    temperature_next=temperature_next\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:21:25.634302Z",
     "start_time": "2019-10-08T13:21:25.630469Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_next   # +/- 100K "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:21:25.798131Z",
     "start_time": "2019-10-08T13:21:25.788508Z"
    }
   },
   "outputs": [],
   "source": [
    "if step_count in step_dict.keys():\n",
    "    temperature_next = step_dict[step_count]['temperature_next']\n",
    "else: \n",
    "    step_dict[step_count]= {'temperature_next': temperature_next}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:21:25.824814Z",
     "start_time": "2019-10-08T13:21:25.799495Z"
    }
   },
   "outputs": [],
   "source": [
    "step_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: set up the interface structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:21:26.283490Z",
     "start_time": "2019-10-08T13:21:25.826153Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_next, temperature_mean, temperature_left, temperature_right, strain_result_lst, pressure_result_lst = calc_temp_iteration(\n",
    "        basis=ham_minimize_vol.get_structure().repeat([1,1,2]), \n",
    "        temperature_next=temperature_next, \n",
    "        project_parameter=project_parameter,\n",
    "        timestep = timestep,\n",
    "        nve_run_time_steps=nve_run_time_steps, \n",
    "        fit_range=fit_range,  \n",
    "        center=None,\n",
    "        debug_plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:26:27.550398Z",
     "start_time": "2019-10-08T13:21:26.285069Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_next, temperature_mean, temperature_left, temperature_right"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: run NVT and NVE MD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:36:31.161068Z",
     "start_time": "2019-10-08T13:26:27.620826Z"
    }
   },
   "outputs": [],
   "source": [
    "output = validate_convergence(\n",
    "    pr=pr,\n",
    "    temperature_left=temperature_left,\n",
    "    temperature_next=temperature_next,\n",
    "    temperature_right=temperature_right,\n",
    "    enable_iteration=enable_iteration,\n",
    "    timestep_iter=timestep_iter,\n",
    "    timestep_lst=project_parameter['timestep_lst'],\n",
    "    timestep=timestep,\n",
    "    fit_range_iter=fit_range_iter,\n",
    "    fit_range_lst=project_parameter['fit_range_lst'],\n",
    "    fit_range=fit_range,\n",
    "    nve_run_time_steps_iter=nve_run_time_steps_iter,\n",
    "    nve_run_time_steps_lst=project_parameter['nve_run_time_steps_lst'],\n",
    "    nve_run_time_steps=nve_run_time_steps,\n",
    "    strain_result_lst=strain_result_lst,\n",
    "    pressure_result_lst=pressure_result_lst,\n",
    "    step_count=step_count,\n",
    "    step_dict=step_dict,\n",
    "    boundary_value=project_parameter['boundary_value'],\n",
    "    ratio_boundary=project_parameter['ratio_boundary'],\n",
    "    convergence_goal=project_parameter['convergence_criterion'],\n",
    "    output_file=output_file\n",
    ")\n",
    "convergence_goal_achieved, enable_iteration, step_count, step_dict, timestep, fit_range, nve_run_time_steps, boundary_value, ratio_boundary, temperature_next, center = output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:36:31.775786Z",
     "start_time": "2019-10-08T13:36:31.234999Z"
    }
   },
   "outputs": [],
   "source": [
    "step_dict[step_count]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:36:32.211084Z",
     "start_time": "2019-10-08T13:36:31.777491Z"
    }
   },
   "outputs": [],
   "source": [
    "convergence_goal_achieved"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: full cycle, predict the final melting point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:36:32.479555Z",
     "start_time": "2019-10-08T13:36:32.461436Z"
    }
   },
   "outputs": [],
   "source": [
    "temperature_estimate_lst, temperature_calculated_lst = [], []\n",
    "while len(temperature_calculated_lst) < 3 or not convergence_goal_achieved and len(temperature_calculated_lst) < 10:\n",
    "    temperature_estimate_lst.append(temperature_next)\n",
    "    temperature_next, temperature_mean, temperature_left, temperature_right, strain_result_lst, pressure_result_lst = calc_temp_iteration(\n",
    "        basis=ham_minimize_vol.get_structure().repeat([1,1,2]), \n",
    "        temperature_next=temperature_next, \n",
    "        project_parameter=project_parameter,\n",
    "        timestep=timestep,\n",
    "        nve_run_time_steps=nve_run_time_steps, \n",
    "        fit_range=fit_range,  \n",
    "        center=center,\n",
    "        debug_plot=True)\n",
    "    print(temperature_next, temperature_mean, temperature_left, temperature_right)\n",
    "    output = validate_convergence(pr=pr,\n",
    "                                  temperature_left=temperature_left, \n",
    "                                  temperature_next=temperature_next, \n",
    "                                  temperature_right=temperature_right, \n",
    "                                  enable_iteration=enable_iteration,\n",
    "                                  timestep_iter=timestep_iter,\n",
    "                                  timestep_lst=project_parameter['timestep_lst'],\n",
    "                                  timestep=timestep,\n",
    "                                  fit_range_iter=fit_range_iter, \n",
    "                                  fit_range_lst=project_parameter['fit_range_lst'], \n",
    "                                  fit_range=fit_range, \n",
    "                                  nve_run_time_steps_iter=nve_run_time_steps_iter, \n",
    "                                  nve_run_time_steps_lst=project_parameter['nve_run_time_steps_lst'], \n",
    "                                  nve_run_time_steps=nve_run_time_steps,\n",
    "                                  strain_result_lst=strain_result_lst, \n",
    "                                  pressure_result_lst=pressure_result_lst, \n",
    "                                  step_count=step_count, \n",
    "                                  step_dict=step_dict, \n",
    "                                  boundary_value=project_parameter['boundary_value'], \n",
    "                                  ratio_boundary=project_parameter['ratio_boundary'],\n",
    "                                  convergence_goal=project_parameter['convergence_criterion'],\n",
    "                                  output_file=output_file)\n",
    "    convergence_goal_achieved, enable_iteration, step_count, step_dict, timestep, fit_range, nve_run_time_steps, boundary_value, ratio_boundary, temperature_next, center = output\n",
    "    print(step_dict[step_count])\n",
    "    temperature_calculated_lst.append(temperature_next)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-10-08T13:36:32.535574Z",
     "start_time": "2019-10-08T13:36:32.480972Z"
    }
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(output_file):\n",
    "    with open(output_file, 'w') as f:\n",
    "        json.dump(step_dict, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the convergence of loop calculations\n",
    "first_key = list(step_dict.keys())[-1]\n",
    "second_key = 'temperature_next'\n",
    "allt = [step_dict[key][second_key] for key in list(step_dict.keys())]\n",
    "plt.plot(np.arange(1, len(allt)), allt[0:-1], 'ro-', label=r\"$T^e$\")\n",
    "plt.plot(np.arange(1, len(allt)), allt[1:], 'bo-', label=r\"$T^p$\")\n",
    "plt.legend(fontsize=14)\n",
    "plt.tick_params(axis='both', labelsize=14)\n",
    "plt.xlabel('Number of loops', fontsize=14)\n",
    "plt.ylabel('Temperature (K)', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pr.job_table()\n",
    "df"
   ]
  }
 ],
 "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "120px",
    "width": "252px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "166.796875px"
   },
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}