Commit d81b18a3 authored by ggandus's avatar ggandus
Browse files

Add Exercise-1

parent 6044c5dc
{
"cells": [
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"from scipy.linalg import eigh\n",
"import nglview\n",
"import sys\n",
"from ase import Atoms\n",
"from ase.calculators.emt import EMT\n",
"from IPython import embed"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exercise-1.aux Exercise-1.pdf \u001b[0m\u001b[01;34mlatex-praesentation\u001b[0m/\r\n",
"Exercise1.ipynb \u001b[01;31mExercise-1.synctex.gz\u001b[0m \u001b[01;31mlatex-praesentation.zip\u001b[0m\r\n",
"Exercise-1.log Exercise-1.tex\r\n",
"Exercise-1.out h2.emt.traj\r\n"
]
}
],
"source": [
"ls"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"class BFGS():\n",
" def __init__(self, atoms, trajectory=None, maxstep=0.04):\n",
" self.maxstep = maxstep\n",
" self.nsteps = 0\n",
" self.logfile = sys.stdout\n",
" self.trajectory = []\n",
" self.initialize(atoms)\n",
" \n",
" def initialize(self, atoms):\n",
" self.H = np.eye(3 * len(atoms)) * 70.0\n",
" self.r0 = None\n",
" self.f0 = None\n",
" \n",
" def step(self, atoms, f):\n",
" r = atoms.get_positions()\n",
" f = f.reshape(-1)\n",
" update(self.H, r.flat, f, self.r0, self.f0)\n",
" dr = determine_step(self.H, f, self.maxstep) \n",
" atoms.set_positions(r + dr)\n",
" self.r0 = r.flat.copy()\n",
" self.f0 = f.copy()\n",
" \n",
" def log(self, atoms, forces):\n",
" fmax = np.linalg.norm(atoms.get_forces(),axis=1).max()\n",
" e = atoms.get_potential_energy(force_consistent=True)\n",
" name = self.__class__.__name__\n",
" if self.nsteps == 0:\n",
" self.logfile.write(\n",
" '{} {:4} {:>10} {:>9}\\n'.format\n",
" (' ' * len(name), 'Step', 'Energy', 'fmax'))\n",
" \n",
" self.logfile.write('{}: {:3} {:10.6} {:>9.4}\\n'.format\n",
" (name, self.nsteps, e, fmax))\n",
" self.logfile.flush()\n",
" \n",
" def run(self):\n",
" step = 0 \n",
" max_nstep = 1000\n",
" while not converged(atoms) or step > max_nstep:\n",
" self.trajectory.append(atoms.copy())\n",
" f=atoms.get_forces()\n",
" self.step(atoms, f)\n",
" self.log(atoms, f)\n",
" self.nsteps+=1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Update the Hessian matrix\n",
"\n",
"$$H-=\\frac{df\\cdot df^T}{dr\\cdot df^T}+\\frac{(H\\cdot dr)\\cdot(H\\cdot dr)^T}{dr\\cdot H\\cdot dr^T}$$\n",
"\n",
"Hints: \n",
"> You might find convenient to define: \n",
"$a=dr\\cdot df^T$ <br>\n",
"$dg=H\\cdot dr^T$ <br>\n",
"$b=dr\\cdot dg^T$ "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def update(H, r, f, r0, f0):\n",
" if r0 is None:\n",
" return\n",
" if np.abs(dr).max() < 1e-7:\n",
" # Same configuration again:\n",
" return\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" H -= #..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute the atomic displacements\n",
"$$dr=\\Big( V\\cdot\\Big[\\frac{V\\cdot f^T}{|\\omega|}\\Big] \\Big).reshape(-1,3)$$\n",
"where $\\omega$ and $V$ are the eigevalues and eigenvectors of $H$. <br><br>\n",
"Impose the condition that \n",
"$$\\text{if}\\Big[ max(|dr|)<maxstep \\Big] \\rightarrow |dr|*=\\frac{max(||dr|)}{maxstep}$$\n",
"Hints:\n",
"> Use the **eigh** function from the *scipy.linalg* library"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def determine_step(H, f, maxstep):\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" return dr"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def converged(atoms, fmax=0.05):\n",
" \"\"\"Did the optimization converge?\"\"\"\n",
" forces = atoms.get_forces()\n",
" return (forces**2).sum(axis=1).max() < fmax**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run the simulation\n",
"\n",
"1. Run the your BFGS algorithm and verify that the total energy of the sistem decreases. \n",
"2. Visualize the trajectory of your geometry optimization with the *ngl* interactive viewer."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Step Energy fmax\n",
"BFGS: 0 5.61262 5.88\n",
"BFGS: 1 5.2235 3.234\n",
"BFGS: 2 5.11248 2.588\n",
"BFGS: 3 5.04177 1.534\n",
"BFGS: 4 4.93205 1.842\n",
"BFGS: 5 4.83166 2.273\n",
"BFGS: 6 4.6923 2.64\n",
"BFGS: 7 4.51955 2.985\n",
"BFGS: 8 4.32152 3.311\n",
"BFGS: 9 4.10189 3.598\n",
"BFGS: 10 3.86448 3.805\n",
"BFGS: 11 3.61641 3.868\n",
"BFGS: 12 3.37098 3.698\n",
"BFGS: 13 3.14882 3.172\n",
"BFGS: 14 2.97719 2.205\n",
"BFGS: 15 2.87464 1.069\n",
"BFGS: 16 2.80678 2.139\n",
"BFGS: 17 2.72527 3.124\n",
"BFGS: 18 2.63984 4.048\n",
"BFGS: 19 2.55637 4.904\n",
"BFGS: 20 2.4827 5.643\n",
"BFGS: 21 2.42973 6.215\n",
"BFGS: 22 2.4058 6.552\n",
"BFGS: 23 2.39017 6.514\n",
"BFGS: 24 2.33427 6.249\n",
"BFGS: 25 2.27681 5.945\n",
"BFGS: 26 2.20893 5.408\n",
"BFGS: 27 2.14038 4.745\n",
"BFGS: 28 2.07244 3.972\n",
"BFGS: 29 2.00762 3.116\n",
"BFGS: 30 1.9513 2.229\n",
"BFGS: 31 1.90805 1.341\n",
"BFGS: 32 1.88322 0.4887\n",
"BFGS: 33 1.8789 0.02878\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0c8bf0c2f2f34cb9b518f027976a3b12",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"NGLWidget(count=34)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "cb7097c159a344ab844a00fb8f4344ea",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Tab(children=(Box(children=(Box(children=(Box(children=(Label(value='step'), IntSlider(value=1, min=-100)), la…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"atoms = Atoms('H2O', positions=[[0.0, 0.0, 0.0],\n",
" [0.0, 1.0, 0.0],\n",
" [0.0, 0.0, 2.0]])\n",
"calc = EMT()\n",
"\n",
"atoms.set_calculator(calc)\n",
"opt = BFGS(atoms, trajectory='h2.traj')\n",
"opt.run();nglview.show_asetraj(opt.trajectory,gui=True)"
]
}
],
"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.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
%% Cell type:code id: tags:
``` python
from scipy.linalg import eigh
import nglview
import sys
from ase import Atoms
from ase.calculators.emt import EMT
from IPython import embed
```
%% Cell type:code id: tags:
``` python
ls
```
%%%% Output: stream
Exercise-1.aux Exercise-1.pdf latex-praesentation/
Exercise1.ipynb Exercise-1.synctex.gz latex-praesentation.zip
Exercise-1.log Exercise-1.tex
Exercise-1.out h2.emt.traj
%% Cell type:code id: tags:
``` python
class BFGS():
def __init__(self, atoms, trajectory=None, maxstep=0.04):
self.maxstep = maxstep
self.nsteps = 0
self.logfile = sys.stdout
self.trajectory = []
self.initialize(atoms)
def initialize(self, atoms):
self.H = np.eye(3 * len(atoms)) * 70.0
self.r0 = None
self.f0 = None
def step(self, atoms, f):
r = atoms.get_positions()
f = f.reshape(-1)
update(self.H, r.flat, f, self.r0, self.f0)
dr = determine_step(self.H, f, self.maxstep)
atoms.set_positions(r + dr)
self.r0 = r.flat.copy()
self.f0 = f.copy()
def log(self, atoms, forces):
fmax = np.linalg.norm(atoms.get_forces(),axis=1).max()
e = atoms.get_potential_energy(force_consistent=True)
name = self.__class__.__name__
if self.nsteps == 0:
self.logfile.write(
'{} {:4} {:>10} {:>9}\n'.format
(' ' * len(name), 'Step', 'Energy', 'fmax'))
self.logfile.write('{}: {:3} {:10.6} {:>9.4}\n'.format
(name, self.nsteps, e, fmax))
self.logfile.flush()
def run(self):
step = 0
max_nstep = 1000
while not converged(atoms) or step > max_nstep:
self.trajectory.append(atoms.copy())
f=atoms.get_forces()
self.step(atoms, f)
self.log(atoms, f)
self.nsteps+=1
```
%% Cell type:markdown id: tags:
### Update the Hessian matrix
$$H-=\frac{df\cdot df^T}{dr\cdot df^T}+\frac{(H\cdot dr)\cdot(H\cdot dr)^T}{dr\cdot H\cdot dr^T}$$
Hints:
> You might find convenient to define:
$a=dr\cdot df^T$ <br>
$dg=H\cdot dr^T$ <br>
$b=dr\cdot dg^T$
%% Cell type:code id: tags:
``` python
def update(H, r, f, r0, f0):
if r0 is None:
return
if np.abs(dr).max() < 1e-7:
# Same configuration again:
return
## TODO
## ...
##
H -= #...
```
%% Cell type:markdown id: tags:
### Compute the atomic displacements
$$dr=\Big( V\cdot\Big[\frac{V\cdot f^T}{|\omega|}\Big] \Big).reshape(-1,3)$$
where $\omega$ and $V$ are the eigevalues and eigenvectors of $H$. <br><br>
Impose the condition that
$$\text{if}\Big[ max(|dr|)<maxstep \Big] \rightarrow |dr|*=\frac{max(||dr|)}{maxstep}$$
Hints:
> Use the **eigh** function from the *scipy.linalg* library
%% Cell type:code id: tags:
``` python
def determine_step(H, f, maxstep):
## TODO
## ...
##
return dr
```
%% Cell type:code id: tags:
``` python
def converged(atoms, fmax=0.05):
"""Did the optimization converge?"""
forces = atoms.get_forces()
return (forces**2).sum(axis=1).max() < fmax**2
```
%% Cell type:markdown id: tags:
### Run the simulation
1. Run the your BFGS algorithm and verify that the total energy of the sistem decreases.
2. Visualize the trajectory of your geometry optimization with the *ngl* interactive viewer.
%% Cell type:code id: tags:
``` python
atoms = Atoms('H2O', positions=[[0.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 2.0]])
calc = EMT()
atoms.set_calculator(calc)
opt = BFGS(atoms, trajectory='h2.traj')
opt.run();nglview.show_asetraj(opt.trajectory,gui=True)
```
%%%% Output: stream
Step Energy fmax
BFGS: 0 5.61262 5.88
BFGS: 1 5.2235 3.234
BFGS: 2 5.11248 2.588
BFGS: 3 5.04177 1.534
BFGS: 4 4.93205 1.842
BFGS: 5 4.83166 2.273
BFGS: 6 4.6923 2.64
BFGS: 7 4.51955 2.985
BFGS: 8 4.32152 3.311
BFGS: 9 4.10189 3.598
BFGS: 10 3.86448 3.805
BFGS: 11 3.61641 3.868
BFGS: 12 3.37098 3.698
BFGS: 13 3.14882 3.172
BFGS: 14 2.97719 2.205
BFGS: 15 2.87464 1.069
BFGS: 16 2.80678 2.139
BFGS: 17 2.72527 3.124
BFGS: 18 2.63984 4.048
BFGS: 19 2.55637 4.904
BFGS: 20 2.4827 5.643
BFGS: 21 2.42973 6.215
BFGS: 22 2.4058 6.552
BFGS: 23 2.39017 6.514
BFGS: 24 2.33427 6.249
BFGS: 25 2.27681 5.945
BFGS: 26 2.20893 5.408
BFGS: 27 2.14038 4.745
BFGS: 28 2.07244 3.972
BFGS: 29 2.00762 3.116
BFGS: 30 1.9513 2.229
BFGS: 31 1.90805 1.341
BFGS: 32 1.88322 0.4887
BFGS: 33 1.8789 0.02878
%%%% Output: display_data
%%%% Output: display_data
{
"cells": [
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"from scipy.linalg import eigh\n",
"import nglview\n",
"import sys\n",
"from ase import Atoms\n",
"from ase.calculators.emt import EMT\n",
"from IPython import embed"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exercise-1.aux Exercise-1.pdf \u001b[0m\u001b[01;34mlatex-praesentation\u001b[0m/\r\n",
"Exercise1.ipynb \u001b[01;31mExercise-1.synctex.gz\u001b[0m \u001b[01;31mlatex-praesentation.zip\u001b[0m\r\n",
"Exercise-1.log Exercise-1.tex\r\n",
"Exercise-1.out h2.emt.traj\r\n"
]
}
],
"source": [
"ls"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"class BFGS():\n",
" def __init__(self, atoms, trajectory=None, maxstep=0.04):\n",
" self.maxstep = maxstep\n",
" self.nsteps = 0\n",
" self.logfile = sys.stdout\n",
" self.trajectory = []\n",
" self.initialize(atoms)\n",
" \n",
" def initialize(self, atoms):\n",
" self.H = np.eye(3 * len(atoms)) * 70.0\n",
" self.r0 = None\n",
" self.f0 = None\n",
" \n",
" def step(self, atoms, f):\n",
" r = atoms.get_positions()\n",
" f = f.reshape(-1)\n",
" update(self.H, r.flat, f, self.r0, self.f0)\n",
" dr = determine_step(self.H, f, self.maxstep) \n",
" atoms.set_positions(r + dr)\n",
" self.r0 = r.flat.copy()\n",
" self.f0 = f.copy()\n",
" \n",
" def log(self, atoms, forces):\n",
" fmax = np.linalg.norm(atoms.get_forces(),axis=1).max()\n",
" e = atoms.get_potential_energy(force_consistent=True)\n",
" name = self.__class__.__name__\n",
" if self.nsteps == 0:\n",
" self.logfile.write(\n",
" '{} {:4} {:>10} {:>9}\\n'.format\n",
" (' ' * len(name), 'Step', 'Energy', 'fmax'))\n",
" \n",
" self.logfile.write('{}: {:3} {:10.6} {:>9.4}\\n'.format\n",
" (name, self.nsteps, e, fmax))\n",
" self.logfile.flush()\n",
" \n",
" def run(self):\n",
" step = 0 \n",
" max_nstep = 1000\n",
" while not converged(atoms) or step > max_nstep:\n",
" self.trajectory.append(atoms.copy())\n",
" f=atoms.get_forces()\n",
" self.step(atoms, f)\n",
" self.log(atoms, f)\n",
" self.nsteps+=1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Update the Hessian matrix\n",
"\n",
"$$H-=\\frac{df^T df}{dr\\cdot df^T}+\\frac{(H\\cdot dr)^T(H\\cdot dr)}{dr\\cdot H\\cdot dr^T}$$\n",
"\n",
"Hints: \n",
"> You might find convenient to define: \n",
"$a=dr\\cdot df^T$ <br>\n",
"$dg=H\\cdot dr^T$ <br>\n",
"$b=dr\\cdot dg^T$ <br>\n",
"Use the **outer** and **dot** functions from numpy library"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def update(H, r, f, r0, f0):\n",
" if r0 is None:\n",
" return\n",
" if np.abs(dr).max() < 1e-7:\n",
" # Same configuration again:\n",
" return\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" H -= #..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute the atomic displacements\n",
"$$dr=\\Big( V\\cdot\\Big[\\frac{V\\cdot f^T}{|\\omega|}\\Big] \\Big).reshape(-1,3)$$\n",
"where $\\omega$ and $V$ are the eigevalues and eigenvectors of $H$. <br><br>\n",
"Impose the condition that \n",
"$$\\text{if}\\Big[ max(|dr|)<maxstep \\Big] \\rightarrow |dr|*=\\frac{max(||dr|)}{maxstep}$$\n",
"Hints:\n",
"> Use the **eigh** function from the *scipy.linalg* library"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def determine_step(H, f, maxstep):\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" return dr"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def converged(atoms, fmax=0.05):\n",
" \"\"\"Did the optimization converge?\"\"\"\n",
" forces = atoms.get_forces()\n",
" return (forces**2).sum(axis=1).max() < fmax**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run the simulation\n",
"\n",
"1. Run the your BFGS algorithm and verify that the total energy of the sistem decreases. \n",
"2. Visualize the trajectory of your geometry optimization with the *ngl* interactive viewer."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Step Energy fmax\n",
"BFGS: 0 5.61262 5.88\n",
"BFGS: 1 5.2235 3.234\n",
"BFGS: 2 5.11248 2.588\n",
"BFGS: 3 5.04177 1.534\n",
"BFGS: 4 4.93205 1.842\n",
"BFGS: 5 4.83166 2.273\n",
"BFGS: 6 4.6923 2.64\n",
"BFGS: 7 4.51955 2.985\n",
"BFGS: 8 4.32152 3.311\n",
"BFGS: 9 4.10189 3.598\n",
"BFGS: 10 3.86448 3.805\n",
"BFGS: 11 3.61641 3.868\n",
"BFGS: 12 3.37098 3.698\n",
"BFGS: 13 3.14882 3.172\n",