Commit f2fdbc71 authored by Robin Worreby's avatar Robin Worreby
Browse files

Added solved ex 1

parent 6e76c0ad
......@@ -2,13 +2,14 @@
"cells": [
{
"cell_type": "code",
"execution_count": 12,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.linalg import eigh\n",
"import nglview\n",
"import sys\n",
"import numpy as np\n",
"from ase import Atoms\n",
"from ase.calculators.emt import EMT\n",
"from IPython import embed"
......@@ -16,27 +17,16 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": null,
"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"
]
}
],
"outputs": [],
"source": [
"ls"
]
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
......@@ -104,20 +94,27 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def update(H, r, f, r0, f0):\n",
" if r0 is None:\n",
" return\n",
" dr = r - r0\n",
" \n",
" if np.abs(dr).max() < 1e-7:\n",
" # Same configuration again:\n",
" return\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" H -= #..."
" \n",
" # Define helpers\n",
" df = f - f0\n",
" a = np.dot(dr, df)\n",
" dg = np.dot(H, dr)\n",
" b = np.dot(dr, dg)\n",
" \n",
" # Updating Hessian for the next iteration\n",
" H -= np.outer(df, df) / a + np.outer(dg, dg) / b"
]
},
{
......@@ -135,20 +132,25 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def determine_step(H, f, maxstep):\n",
" ## TODO\n",
" ## ...\n",
" ##\n",
" w, v = eigh(H)\n",
" dr = np.dot(v, np.dot(v, f) / np.fabs(w)).reshape(-1, 3)\n",
" \n",
" steplength = (dr**2).sum(1)**0.5\n",
" maxsteplength = np.max(steplength)\n",
" if(maxsteplength >= maxstep):\n",
" dr *= maxstep / maxsteplength\n",
"\n",
" return dr"
]
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
......@@ -164,85 +166,15 @@
"source": [
"### Run the simulation\n",
"\n",
"1. Run the your BFGS algorithm and verify that the total energy of the sistem decreases. \n",
"1. Run the your BFGS algorithm and verify that the total energy of the system decreases. \n",
"2. Visualize the trajectory of your geometry optimization with the *ngl* interactive viewer."
]
},
{
"cell_type": "code",
"execution_count": 21,
"execution_count": null,
"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"
}
],
"outputs": [],
"source": [
"atoms = Atoms('H2O', positions=[[0.0, 0.0, 0.0],\n",
" [0.0, 1.0, 0.0],\n",
......@@ -253,6 +185,20 @@
"opt = BFGS(atoms, trajectory='h2.traj')\n",
"opt.run();nglview.show_asetraj(opt.trajectory,gui=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
......@@ -271,7 +217,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
"version": "3.7.2"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
from scipy.linalg import eigh
import nglview
import sys
import numpy as np
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^T df}{dr\cdot df^T}+\frac{(H\cdot dr)^T(H\cdot dr)}{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$ <br>
Use the **outer** and **dot** functions from numpy library
%% Cell type:code id: tags:
``` python
def update(H, r, f, r0, f0):
if r0 is None:
return
dr = r - r0
if np.abs(dr).max() < 1e-7:
# Same configuration again:
return
## TODO
## ...
##
H -= #...
# Define helpers
df = f - f0
a = np.dot(dr, df)
dg = np.dot(H, dr)
b = np.dot(dr, dg)
# Updating Hessian for the next iteration
H -= np.outer(df, df) / a + np.outer(dg, dg) / b
```
%% 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
## ...
##
w, v = eigh(H)
dr = np.dot(v, np.dot(v, f) / np.fabs(w)).reshape(-1, 3)
steplength = (dr**2).sum(1)**0.5
maxsteplength = np.max(steplength)
if(maxsteplength >= maxstep):
dr *= maxstep / maxsteplength
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.
1. Run the your BFGS algorithm and verify that the total energy of the system 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
%% Cell type:code id: tags:
``` python
```
%%%% Output: display_data
%% Cell type:code id: tags:
``` python
```
......
......@@ -2,13 +2,14 @@
"cells": [
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from scipy.linalg import eigh\n",
"import nglview\n",
"import sys\n",
"import numpy as np\n",
"from ase import Atoms\n",
"from ase.calculators.emt import EMT\n",
"from IPython import embed"
......@@ -16,16 +17,14 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exercise-1.aux Exercise-1.out Exercise-1.tex \u001b[0m\u001b[01;34mlatex-praesentation\u001b[0m/\r\n",
"Exercise1.ipynb Exercise-1.pdf \u001b[01;31mExercise-1.zip\u001b[0m \u001b[01;31mlatex-praesentation.zip\u001b[0m\r\n",
"Exercise-1.log \u001b[01;31mExercise-1.synctex.gz\u001b[0m h2.emt.traj Solution1.ipynb\r\n"
"Exercise1.ipynb Exercise1.pdf Solution1.ipynb\r\n"
]
}
],
......@@ -35,7 +34,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
......@@ -103,7 +102,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
......@@ -138,7 +137,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
......@@ -154,7 +153,7 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
......@@ -176,7 +175,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 7,
"metadata": {},
"outputs": [
{
......@@ -223,7 +222,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "67dca8c993d64960a3c20eb53bed55a2",
"model_id": "46c616eb04ce43dd8df206508e3aed41",
"version_major": 2,
"version_minor": 0
},
......@@ -237,7 +236,7 @@
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "498484ec33f74012b5c397564d9c6957",
"model_id": "870201c6f1394ff3bb33d98bd9ee19da",
"version_major": 2,
"version_minor": 0
},
......@@ -277,7 +276,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
"version": "3.7.2"
}
},
"nbformat": 4,
......
%% Cell type:code id: tags:
``` python
from scipy.linalg import eigh
import nglview
import sys
import numpy as np
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.out Exercise-1.tex latex-praesentation/
Exercise1.ipynb Exercise-1.pdf Exercise-1.zip latex-praesentation.zip
Exercise-1.log Exercise-1.synctex.gz h2.emt.traj Solution1.ipynb
Exercise1.ipynb Exercise1.pdf Solution1.ipynb
%% 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^T df}{dr\cdot df^T}+\frac{(H\cdot dr)^T(H\cdot dr)}{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$ <br>
Use the **outer** and **dot** functions from numpy library
%% Cell type:code id: tags:
``` python
def update(H, r, f, r0, f0):
if r0 is None:
return
dr = r - r0
if np.abs(dr).max() < 1e-7:
# Same configuration again:
return
df = f - f0
a = np.dot(dr, df)
dg = np.dot(H, dr)
b = np.dot(dr, dg)
H -= np.outer(df, df) / a + np.outer(dg, dg) / b
```
%% 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):
omega, V = eigh(H)
dr = np.dot(V, np.dot(f, V) / np.fabs(omega)).reshape((-1, 3))
steplengths = (dr**2).sum(1)**0.5
maxsteplength = np.max(steplengths)
if maxsteplength >= maxstep:
dr *= maxstep / maxsteplength
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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment