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 c6f287f8 authored by spiasko's avatar spiasko
Browse files

Update everything

parent d4cd2666
%% 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^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: display_data
%%%% Output: display_data
This diff is collapsed.
......@@ -277,7 +277,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
"version": "3.6.7"
}
},
"nbformat": 4,
......
1 -0.0604878400
2 -0.0604174878
3 -0.0603416327
4 -0.0604362379
5 -0.0607472659
6 -0.0612509506
7 -0.0618208462
2 -0.0604066276
3 -0.0603506836
4 -0.0604841253
5 -0.0608128447
6 -0.0613056318
7 -0.0618487420
8 -0.0621047902
This diff is collapsed.
38
i = 41, E = -0.0604878400
Ar 15.3022824059 7.0712009285 13.8866605262
Ar 5.1812404922 11.7766857322 12.6696350355
Ar 6.9109487557 8.8302505552 7.0431911232
Ar 13.4873612337 5.1184805322 6.9368194493
Ar 10.8711489775 11.5073690735 15.9888427368
Ar 8.3087544060 4.8108581758 12.2356906303
Ar 5.4796148406 12.4018711529 7.6892292011
Ar 12.9235042175 8.5391378688 5.4569310898
Ar 8.7440808296 13.0637255213 13.2140085468
Ar 14.9942681676 12.9026279329 10.9995454961
Ar 9.5383452908 8.8333323012 3.9407622799
Ar 7.9340772593 11.9519167770 4.9378127621
Ar 6.6366075251 14.5924440901 10.5335452764
Ar 12.4247407298 13.8611357262 13.5373823026
Ar 16.0131101821 7.7691051074 7.4468939980
Ar 15.9554181575 9.2224608565 10.8087274642
Ar 12.8220826518 14.5326135173 8.2551279252
Ar 11.7751187780 12.4180445931 5.1358752224
Ar 10.3437582211 15.5325494626 10.9179533252
Ar 9.2005195461 13.6385589476 7.9405820833
Ar 14.1120671638 11.0865229021 7.8689251333
Ar 7.2502663901 10.4091250686 15.5333404516
Ar 9.9719000712 6.4185874853 6.7123754043
Ar 4.3678582013 9.4519991101 9.7468390154
Ar 12.2699953125 8.0433134316 16.0404328746
Ar 14.8751313672 5.4431692354 10.4913001770
Ar 6.4381473826 8.1585524756 12.5562721116
Ar 5.9346522567 5.9773571608 9.3941237152
Ar 11.1323478713 4.4577364514 9.7103610174
Ar 11.7826231136 6.5351315174 12.6704780776
Ar 8.7862951311 6.9093851869 15.2923678181
Ar 11.4085434989 11.8335421467 10.5031714086
Ar 13.3852404549 10.2222589211 13.2113596426
Ar 10.2618787101 10.1210027348 7.4847726065
Ar 12.6568297469 8.1128276544 9.3918267920
Ar 7.7947652426 10.9915755859 10.1885168579
Ar 9.8334503988 9.5924393695 12.9937339915
Ar 9.1328000184 7.6314567102 10.0321204304
......@@ -117,11 +117,11 @@
 
 
%% Cell type:code id: tags:
 
``` python
de=0. #(for task 1, no interaction)
de = -0.02 #(for task 1, no interaction)
T=220
kb=physical_constants['Boltzmann constant in eV/K'][0]
beta =1.0/(kb*T)
nouter=100
ninner=10
......@@ -146,11 +146,11 @@
## #Molecules
 
%% Cell type:code id: tags:
 
``` python
coverage= #(look in paper lowest coverage)
coverage= 0.1
nmolecules=int(nx*nx*coverage)
```
 
%% Cell type:markdown id: tags:
 
......@@ -172,10 +172,33 @@
if newmolecule not in molecules:
molecules.append(newmolecule)
i=i+1
#### END Create initial geometry
 
def n(mol, molecules):
n+=0
# check if neighbors are occupied
if mol[2] == 1:
if [mol[0],mol[1],0] in molecules:
n += 1
if [mol[0],mol[1]-1,0] in molecules:
n += 1
if [mol[0]+1,mol[1]-1,0] in molecules:
n += 1
if [mol[0]+1,mol[1],0] in molecules:
n += 1
if mol[2] == 0:
if [mol[0],mol[1],0] in molecules:
n += 1
if [mol[0]-1,mol[1],0] in molecules:
n += 1
if [mol[0]-1,mol[1]+1,0] in molecules:
n += 1
if [mol[0],mol[1]+1,0] in molecules:
n += 1
return n
#### Initial energy
totene=0.0
for m in range(nmolecules):
n,l =neighbors(molecules,m,nx,ny)
totene+=n*de
......@@ -211,10 +234,30 @@
6. Update the total energy, totene, and geometry, molecules,
if necessary.
'''
# 1.
k = np.random.randint(0,len(molecules))
mol = molecules[k]
# 2.
while True:
xnew = np.random.randint(0,nx)
ynew = np.random.randint(0,ny)
snew = np.random.randint(0,2)
mol = [xnew, ynew, snew]
# 3.
if mol in molecules:
pass
else break
# 4.
e = de*n(mol,molecules)
nsteps+=1
if nsteps==maxinner:
print('coverage too big')
exit()
......@@ -17965,7 +17965,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
"version": "3.6.7"
}
},
"nbformat": 4,
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