Lec_22.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physics 300 \n",
    "## Computational Physics I (Fall 2017)\n",
    "## BPB-248, Tues/Thurs 10:00 - 11:15 am \n",
    "\n",
    "|Instructor| Prof. Qiang Zhu|\n",
    "|--|-------------------------------|\n",
    "|Email | qiang.zhu@unlv.edu|\n",
    "|Website|http://www.physics.unlv.edu/~qzhu/|\n",
    "|Office| BPB 232|\n",
    "|Office hours | Tues/Thurs 8:30 - 10:00 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 22 Global Optmization (III)\n",
    "\n",
    "## 22.1 Genetic Algorithm\n",
    "\n",
    "A genetic algorithm is a search heuristic that is inspired by Charles Darwin’s theory of natural evolution. This algorithm reflects the process of natural selection where the fittest individuals are selected for reproduction in order to produce offspring of the next generation.\n",
    "\n",
    "- 1, Initialization\n",
    "- 2, Selection\n",
    "- 3, Genetic operators\n",
    "- 4, Termination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 22.2 Application to LJ clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def LJ(r):\n",
    "    r6 = r**6\n",
    "    r12 = r6*r6\n",
    "    return 4*(1/r12 - 1/r6)\n",
    "\n",
    "def total_energy(positions):\n",
    "    \"\"\"\n",
    "    Calculate the total energy\n",
    "    input:\n",
    "    positions: 3*N array which represents the atomic positions\n",
    "    output\n",
    "    E: the total energy\n",
    "    \"\"\"\n",
    "    E = 0\n",
    "    N_atom = int(len(positions)/3)\n",
    "\n",
    "    #positions = [x0, y0, z0, x1, y1, z1, .....  , xn, yn, zn]\n",
    "    for i in range(N_atom-1):\n",
    "        for j in range(i+1, N_atom):\n",
    "            pos1 = positions[i*3:(i+1)*3]\n",
    "            pos2 = positions[j*3:(j+1)*3]\n",
    "            #print('pos1:  ', pos1)\n",
    "            #print('pos2:  ', pos2)\n",
    "            dist = np.linalg.norm(pos1-pos2)\n",
    "            #print(i,j, dist)\n",
    "            E += LJ(dist)\n",
    "    return E\n",
    "            \n",
    "def init_pos(N, L=5):\n",
    "    return L*np.random.random_sample((N*3,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "import numpy as np\n",
    "\n",
    "def init_population(population_size, N_atom):\n",
    "    cluster = []\n",
    "    for i in range(population_size):\n",
    "        cluster.append(init_pos(N_atom))\n",
    "    return np.array(cluster)\n",
    "\n",
    "def local_optimize(cluster):\n",
    "    fitness = []\n",
    "    for i, cluster0 in enumerate(cluster):\n",
    "        res = minimize(total_energy, cluster0, method='CG', tol=1e-3) \n",
    "        cluster0 = res.x\n",
    "        fitness.append(res.fun)\n",
    "    return cluster, np.array(fitness)\n",
    "        \n",
    "def ranking(fitness):\n",
    "    return np.argsort(fitness)\n",
    "\n",
    "def mutation(cluster, rank, kT=1):\n",
    "    \n",
    "    id = int(len(rank)*np.random.random()/2)\n",
    "    cluster0 = cluster[id]\n",
    "    return cluster0 + kT*np.random.random_sample((len(cluster0),))\n",
    "\n",
    "def crossover(cluster, rank):\n",
    "    id1 = int(len(rank)*np.random.random()/2)\n",
    "    while True:\n",
    "        id2 = int(len(rank)*np.random.random()/2)\n",
    "        if id2 != id1:\n",
    "            break\n",
    "    frac = np.random.random()\n",
    "    return cluster[id1]*frac + cluster[id2]*(1-frac)\n",
    "\n",
    "\n",
    "def GA(generation_num=10, population_num=10, atom_num=10):\n",
    "    for i in range(generation_num):\n",
    "        if i == 0:\n",
    "            cluster = init_population(population_num, atom_num)\n",
    "            \n",
    "        cluster, fitness = local_optimize(cluster)\n",
    "        #print(fitness)\n",
    "        print('best fitness in generation ', i, ':  ', min(fitness))\n",
    "        \n",
    "        rank = ranking(fitness)\n",
    "        new_cluster = []\n",
    "        for j in range(population_num):\n",
    "            if j < int(0.7*population_num):\n",
    "                new_cluster.append(crossover(cluster, rank))\n",
    "            else:\n",
    "                new_cluster.append(mutation(cluster, rank))\n",
    "        cluster = new_cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best fitness in generation  0 :   -26.5580108882\n",
      "best fitness in generation  1 :   -28.4225318863\n",
      "best fitness in generation  2 :   -27.4468289166\n",
      "best fitness in generation  3 :   -27.4797388754\n",
      "best fitness in generation  4 :   -26.5955206605\n",
      "best fitness in generation  5 :   -27.545206841\n",
      "best fitness in generation  6 :   -25.5040078792\n",
      "best fitness in generation  7 :   -27.5223300827\n",
      "best fitness in generation  8 :   -28.4225318703\n",
      "best fitness in generation  9 :   -27.4468289373\n"
     ]
    }
   ],
   "source": [
    "GA(atom_num=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 22.3 Homework\n",
    "\n",
    "Check the following documentation and use the differential evolution method to perform optimization on LJ clusters\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html#scipy.optimize.differential_evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}