To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit a20fdbbc authored by ggandus's avatar ggandus
Browse files

No commit message

No commit message
parent b2ba60bf
%% 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
```
%% 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: display_data
%%%% Output: display_data
Markdown is supported
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