Small doc+tutorial updates

parent 9f2bc849
Pipeline #12095156 passed with stage
in 6 minutes and 23 seconds
......@@ -7,6 +7,12 @@ Wave functions
.. autoclass:: gpaw.wavefunctions.fd.FDWaveFunctions
:members:
.. autoclass:: gpaw.wavefunctions.pw.PWWaveFunctions
:members:
.. autoclass:: gpaw.wavefunctions.pw.PW
:members:
.. autoclass:: gpaw.wavefunctions.lcao.LCAOWaveFunctions
:members:
......@@ -15,4 +21,3 @@ Wave functions
.. autoclass:: gpaw.eigensolvers.eigensolver.Eigensolver
:members:
from __future__ import print_function
from ase import Atoms
from ase.build import bulk
from gpaw import GPAW
from gpaw import PW
from gpaw import GPAW, PW
cell = bulk('Si', 'fcc', a=5.421).get_cell()
a = Atoms('Si2', cell=cell, pbc=True,
scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))
a = 5.421
si = bulk('Si', 'fcc', a=a)
# or equivalently:
# b = a / 2
# from ase import Atoms
# si = Atoms('Si2', cell=[[0, b, b], [b, 0, b], [b, b, 0]], pbc=True,
# scaled_positions=[[0, 0, 0], [0.25, 0.25, 0.25]])
for x in [100, 200, 300, 400, 500, 600, 700, 800]:
# for x in [0.24, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.1]:
......@@ -16,6 +17,6 @@ for x in [100, 200, 300, 400, 500, 600, 700, 800]:
kpts=(4, 4, 4),
txt='convergence_%s.txt' % x)
a.set_calculator(calc)
si.calc = calc
print(x, a.get_potential_energy())
print(x, si.get_potential_energy())
from __future__ import print_function
from ase import Atoms
import numpy as np
from ase.build import bulk
from ase.optimize.bfgs import BFGS
from ase.constraints import UnitCellFilter
from gpaw import GPAW
from gpaw import PW
import numpy as np
cell = bulk('Si', 'fcc', a=6.0).get_cell()
si = bulk('Si', 'fcc', a=6.0)
# Experimental Lattice constant is a=5.421 A
a = Atoms('Si2', cell=cell, pbc=True,
scaled_positions=((0, 0, 0), (0.25, 0.25, 0.25)))
calc = GPAW(xc='PBE',
mode=PW(400),
# mode=PW(400, cell=a.get_cell()), # fix no of planewaves!
kpts=(4, 4, 4),
# convergence={'eigenstates': 1.e-10}, # converge tightly!
txt='stress.txt')
a.set_calculator(calc)
si.calc = GPAW(xc='PBE',
mode=PW(400, dedecut='estimate'),
kpts=(4, 4, 4),
# convergence={'eigenstates': 1.e-10}, # converge tightly!
txt='stress.txt')
uf = UnitCellFilter(a)
uf = UnitCellFilter(si)
relax = BFGS(uf)
relax.run(fmax=0.05) # Consider much tighter fmax!
a = np.dot(a.get_cell()[0], a.get_cell()[0])**0.5 * 2**0.5
print('Relaxed lattice parameter: a = %s A' % a)
a = np.linalg.norm(si.cell[0]) * 2**0.5
print('Relaxed lattice parameter: a = {} Ang'.format(a))
.. _stress:
==========================================================
=================================
Plane wave mode and Stress tensor
==========================================================
=================================
The major advantage of running DFT calculations on a real space grid, is a
very efficient parallelization scheme when dealing with large systems.
However, for small systems it is often faster to use a plane wave basis set
instead. In this case all quantities are represented by their Fourier
transforms on the periodic super cell and periodic boundary conditions are
required. In grid mode, the key convergence parameter is the grid spacing
`h`, whereas in planewave mode the corresponding parameter is
`E_{cut}=G_{cut}^2/2` (in atomic units). `G_{cut}` determines the maximum
size of reciprocal lattice vectors to be included in the plane wave
expansion.
The major advantage of running DFT calculations on a real space grid, is a very efficient parallelization scheme when dealing with large systems. However, for small systems it is often faster to use a plane wave basis set instead. In this case all quantities are represented by their Fourier transforms on the periodic super cell and periodic boundary conditions are required. In grid mode, the key convergence parameter is the grid spacing `h`, whereas in planewave mode the corresponding parameter is `E_{cut}=G_{cut}^2/2` (in atomic units). `G_{cut}` determines the maximum size of reciprocal lattice vectors to be included in the plane wave expansion.
Converging the plane wave cutoff
-------------------------------------
--------------------------------
As a simple example we perform a calculation on bulk Si using the plane wave basis. The following script tests the convergence of the total energy with respect to the cutoff parameter `E_{cut}`.
As a simple example we perform a calculation on bulk Si using the plane wave
basis. The following script tests the convergence of the total energy with
respect to the cutoff parameter `E_{cut}`.
.. literalinclude:: con_pw.py
Note that we use a rather rough k-point sampling and the default Fermi smearing of 0.1 eV. These parameters should of course be converged, but for now we will just keep them fixed in order to speed up the calculations. At `E_{cut}` = 400 eV, the total energy is seen to be converged to within 5 meV. Now edit the script such that the calculation is performed in grid mode with various values of the grid spacing (comment out the lines :samp:`for x in [100, ...` and :samp:`mode=PW(x)` and uncomment the lines :samp:`for h in [0.24, ...` and :samp:`h=x`).
Note that we use a rather rough k-point sampling and the default Fermi
smearing of 0.1 eV. These parameters should of course be converged, but for
now we will just keep them fixed in order to speed up the calculations. At
`E_{cut}` = 400 eV, the total energy is seen to be converged to within 5 meV.
Now edit the script such that the calculation is performed in grid mode with
various values of the grid spacing (comment out the lines :samp:`for x in
[100, ...` and :samp:`mode=PW(x)` and uncomment the lines :samp:`for h in
[0.24, ...` and :samp:`h=x`).
What grid spacing is needed in order to converge the total energy to within 5
meV? What is the timing of this calculation compared to the one at 400 eV in
plane wave mode? Compare the number of grid points and plane waves used in
these two calculations respectively (look in the output file
'convergence_400.txt' and the corresponding one for the grid calculation.)
You can see how long a calculation took by looking at the times recorded in
the .txt file.
What grid spacing is needed in order to converge the total energy to within 5 meV? What is the timing of this calculation compared to the one at 400 eV in plane wave mode? Compare the number of grid points and plane waves used in these two calculations respectively (look in the output file 'convergence_400.txt' and the corresponding one for the grid calculation.) You can see how long a calculation took by looking at the times recorded in the .txt file.
Optimizing the unit cell
------------------------
**Warning**: due to difficulties in optimizing cell and positions simultaneously
:class:`ase.constraints.UnitCellFilter` may produce
incorrect results. Always verify obtained structures by means of
performing separate cell
(see :class:`ase.constraints.StrainFilter`)
and positions optimizations (see :mod:`ase.optimize`).
Consider much more tighter fmax than the one used in this tutorial!
**Warning**: due to difficulties in optimizing cell and positions
simultaneously :class:`ase.constraints.UnitCellFilter` may produce incorrect
results. Always verify obtained structures by means of performing separate
cell (see :class:`ase.constraints.StrainFilter`) and positions optimizations
(see :mod:`ase.optimize`). Consider a much more tight fmax than the one used
in this tutorial!
In the :ref:`aluminium_exercise` exercise the lattice constant of bulk Al was found by calculating the total energy at various lattice distances. A nice feature e of the plane wave mode is that it allows a simple implementation of the stress tensor, which can be used to optimize unit unit cells of periodic systems directly. The following script performs such an optimization for bulk Si.
In the :ref:`aluminium_exercise` exercise the lattice constant of bulk Al was
found by calculating the total energy at various lattice distances. A nice
feature e of the plane wave mode is that it allows a simple implementation of
the stress tensor, which can be used to optimize unit unit cells of periodic
systems directly. The following script performs such an optimization for bulk
Si.
.. literalinclude:: stress.py
The calculation uses 12 iterations to find the optimal lattice constant and the relaxation can be viewed with the command line tool :program:`ase`:
The calculation uses 12 iterations to find the optimal lattice constant and
the relaxation can be viewed with the command line tool :program:`ase`:
.. highlight:: bash
......@@ -40,4 +71,8 @@ The calculation uses 12 iterations to find the optimal lattice constant and the
$ ase gui stress.txt
Since we know the experimental lattice constant, we could probably have calculated the PBE lattice constant faster by fitting a parabola to five points in the vicinity of the expermental lattice constant. However, for complicated unit cells with more than one lattice parameter, the stress tensor becomes a highly valuable tool.
Since we know the experimental lattice constant, we could probably have
calculated the PBE lattice constant faster by fitting a parabola to five
points in the vicinity of the expermental lattice constant. However, for
complicated unit cells with more than one lattice parameter, the stress
tensor becomes a highly valuable tool.
......@@ -37,7 +37,7 @@ class PW(Mode):
ecut: float
Plane-wave cutoff in eV.
dedecut: float or None
dedecut: float or None or 'estimate'
Estimate of derivative of total energy with respect to
plane-wave cutoff. Used to calculate pulay_stress.
pulay_stress: float or None
......
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