Commit 0c641c34 authored by Ask Hjorth Larsen's avatar Ask Hjorth Larsen

merge from master resolving conflict

parents 41516066 40fe01f9
Pipeline #21312980 passed with stage
in 2 minutes and 40 seconds
......@@ -5,9 +5,9 @@ master:
- pip install flake8
- pip install git+https://gitlab.com/ase/ase.git@master > ase-install.output
- sed -i "s/scalapack = False/scalapack = True/" customize.py
- python setup.py install --user > gpaw-install.output
- export PATH=~/.local/bin:$PATH
- gpaw install-data --register gpaw-datasets
- gpaw info
- ase build H2 -V1 | gpaw -P 2 run -p mode=pw
- python setup.py install > gpaw-install.output
- useradd -m user
- su user -c 'gpaw install-data --register gpaw-datasets'
- su user -c 'gpaw info'
- su user -c 'ase build H2 -V1 | gpaw -P 2 run -p mode=pw'
- python -m gpaw.test.pyflakes_check
......@@ -12,9 +12,10 @@ Webpage: http://wiki.fysik.dtu.dk/gpaw
Requirements
------------
* Python_ 2.6-3.5
* Python_ 2.7, 3.4-
* ASE_ (atomic simulation environment)
* NumPy_ (base N-dimensional array package)
* SciPy_ (library for scientific computing)
* LibXC
* BLAS
* LAPACK
......@@ -23,28 +24,27 @@ Optional:
* MPI
* ScaLAPACK
* SciPy_ (library for scientific computing)
Installation
------------
Do this::
$ python setup.py install --user
and make sure you have ``~/.local/bin`` in your $PATH.
For more details, please see:
https://wiki.fysik.dtu.dk/gpaw/install.html
Testing
-------
Please run the tests::
$ gpaw test -j 4 # takes 1 hour!
and send us the output if there are failing tests.
......@@ -63,7 +63,7 @@ Example
-------
Geometry optimization of hydrogen molecule:
>>> from ase import Atoms
>>> from ase.optimize import BFGS
>>> from ase.io import write
......@@ -92,7 +92,7 @@ Getting started
Once you have familiarized yourself with ASE_ and NumPy_, you should take a
look at the GPAW exercises_ and tutorials_.
.. _Python: http://www.python.org/
.. _ASE: http://wiki.fysik.dtu.dk/ase
.. _NumPy: http://docs.scipy.org/doc/numpy/reference/
......
......@@ -77,6 +77,7 @@ PyObject* lxcXCFuncNum(PyObject *self, PyObject *args);
PyObject* exterior_electron_density_region(PyObject *self, PyObject *args);
PyObject* plane_wave_grid(PyObject *self, PyObject *args);
PyObject* tci_overlap(PyObject *self, PyObject *args);
PyObject *pwlfc_expand(PyObject *self, PyObject *args);
PyObject* overlap(PyObject *self, PyObject *args);
PyObject* vdw(PyObject *self, PyObject *args);
PyObject* vdw2(PyObject *self, PyObject *args);
......@@ -127,6 +128,19 @@ PyObject* libvdwxc_init_mpi(PyObject* self, PyObject* args);
PyObject* libvdwxc_init_pfft(PyObject* self, PyObject* args);
#endif // GPAW_WITH_LIBVDWXC
#ifdef GPAW_GITHASH
// For converting contents of a macro to a string, see
// https://en.wikipedia.org/wiki/C_preprocessor#Token_stringification
#define STR(s) #s
#define XSTR(s) STR(s)
PyObject* githash(PyObject* self, PyObject* args)
{
return Py_BuildValue("s", XSTR(GPAW_GITHASH));
}
#undef XSTR
#undef STR
#endif // GPAW_GITHASH
// Moving least squares interpolation
PyObject* mlsqr(PyObject *self, PyObject *args);
......@@ -171,6 +185,7 @@ static PyMethodDef functions[] = {
{"utilities_vdot_self", utilities_vdot_self, METH_VARARGS, 0},
{"eed_region", exterior_electron_density_region, METH_VARARGS, 0},
{"plane_wave_grid", plane_wave_grid, METH_VARARGS, 0},
{"pwlfc_expand", pwlfc_expand, METH_VARARGS, 0},
{"erf", errorfunction, METH_VARARGS, 0},
{"cerf", cerf, METH_VARARGS, 0},
{"pack", pack, METH_VARARGS, 0},
......@@ -245,6 +260,9 @@ static PyMethodDef functions[] = {
{"libvdwxc_init_pfft", libvdwxc_init_pfft, METH_VARARGS, 0},
#endif // GPAW_WITH_LIBVDWXC
{"mlsqr", mlsqr, METH_VARARGS, 0},
#ifdef GPAW_GITHASH
{"githash", githash, METH_VARARGS, 0},
#endif // GPAW_GITHASH
{0, 0, 0, 0}
};
......
......@@ -40,4 +40,82 @@ PyObject *plane_wave_grid(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
PyObject *pwlfc_expand(PyObject *self, PyObject *args)
{
PyArrayObject *G_Qv_obj, *pos_av_obj; //*aji1i2I1I2_obj;
PyObject *lf_aj_obj;
PyArrayObject *Y_LG_obj;
int q, G1, G2;
PyArrayObject *emiGR_G_obj, *f_IG_obj;
if (!PyArg_ParseTuple(args, "OOOOiiiOO", &G_Qv_obj, &pos_av_obj,
&lf_aj_obj, &Y_LG_obj,
&q, &G1, &G2, &f_IG_obj, &emiGR_G_obj))
return NULL;
if(q == -1) {
q = 0;
}
double *G_Qv = PyArray_DATA(G_Qv_obj);
double *pos_av = PyArray_DATA(pos_av_obj);
double complex *emiGR_G = PyArray_DATA(emiGR_G_obj);
double complex *f_IG = PyArray_DATA(f_IG_obj);
double *Y_LG = PyArray_DATA(Y_LG_obj);
int natoms = PyArray_SIZE(pos_av_obj) / 3;
int nG = PyArray_SIZE(G_Qv_obj) / 3;
if(nG != G2 - G1) {
return NULL;
}
int G;
int v;
int L;
//double complex minus_i = 1.0;//-1.0 * I;
double complex imag_powers[4] = {1.0, -I, -1.0, I};
int I1 = 0;
int a;
npy_intp* Ystrides = PyArray_STRIDES(Y_LG_obj);
int nGtotal = Ystrides[0] / sizeof(double);
for(a=0; a < natoms; a++) {
for(G=0; G < nG; G++) {
double GR = 0;
for(v=0; v < 3; v++) {
GR += G_Qv[3 * G + v] * pos_av[3 * a + v];
}
emiGR_G[G] = cexp(-I * GR);
}
PyObject *lf_j_obj = PyList_GET_ITEM(lf_aj_obj, a);
int nj = PyList_GET_SIZE(lf_j_obj);
int j;
for(j=0; j < nj; j++) {
PyObject *l_and_f_qG = PyList_GET_ITEM(lf_j_obj, j);
PyObject *l_obj = PyTuple_GET_ITEM(l_and_f_qG, 0);
PyObject *f_qG_obj = PyTuple_GET_ITEM(l_and_f_qG, 1);
int l = PyLong_AsLong(l_obj);
int nL = 2 * l + 1;
int L2 = l * l;
double complex ilpow = imag_powers[l % 4];
PyObject *f_G_obj = PyList_GET_ITEM(f_qG_obj, q);
double *f_G = PyArray_DATA((PyArrayObject *)f_G_obj);
for(G=0; G < nG; G++) {
double complex emiGR_f = emiGR_G[G] * f_G[G1 + G] * ilpow;
for(L=0; L < nL; L++) {
// Terrible memory locality!
f_IG[(I1 + L) * nG + G] = emiGR_f * Y_LG[(L2 + L) * nGtotal + G1 + G];
}
}
I1 += nL;
}
}
Py_RETURN_NONE;
}
......@@ -50,17 +50,19 @@ if scalapack:
# - static linking:
if 0:
include_dirs += ['/home/user/libxc-4.0.1/include']
extra_link_args += ['/home/user/libxc-4.0.1/lib/libxc.a']
xc = '/home/user/libxc-4.0.4/'
include_dirs += [xc + 'include']
extra_link_args += [xc + 'lib/libxc.a']
if 'xc' in libraries:
libraries.remove('xc')
# - dynamic linking (requires rpath or setting LD_LIBRARY_PATH at runtime):
if 0:
include_dirs += ['/home/user/libxc-4.0.1/include']
library_dirs += ['/home/user/libxc-4.0.1/lib']
xc = '/home/user/libxc-4.0.4/'
include_dirs += [xc + 'include']
library_dirs += [xc + 'lib']
# You can use rpath to avoid changing LD_LIBRARY_PATH:
extra_link_args += ['-Wl,-rpath=/home/user/libxc-4.0.1/lib']
extra_link_args += ['-Wl,-rpath={xc}/lib'.format(xc=xc)]
if 'xc' not in libraries:
libraries.append('xc')
......
......@@ -48,7 +48,7 @@ Common sources of bugs
return a
x = np.ones(5)
print double(x) # x[:] is now 2.
print(double(x)) # x[:] is now 2.
- Forgetting a ``n += 1`` statement in a for loop::
......
......@@ -7,7 +7,7 @@ Debugging
Python debugging
================
Even though some debugging can done just with print statements, a real
Even though some debugging can done just with the print function, a real
debugger offers several advantages. It is possible, for example, to set
breakpoints in certain files or functions, execute the code step by step,
examine and change values of variables. Python contains a standard debugger
......
......@@ -43,18 +43,18 @@ A GPAW instance has the following attributes: ``atoms``,
``occupations``, ``initialized``, ``world`` and ``observers``.
The :class:`gpaw.calculator.GPAW` inherits from:
* :class:`ase.calculators.calculator.Calculator`
This implements the ASE :ref:`ase:calculator interface`
* :class:`gpaw.paw.PAW`
Mixin that adds a number of extra methods without bloating the
:git:`gpaw/calculator.py` file too much.
.. note::
GPAW uses atomic units internally (`\hbar=e=m=1`) and ASE uses
Angstrom and eV (:mod:`~ase.units`).
......@@ -68,12 +68,12 @@ When a GPAW instance is created from scratch::
the GPAW object is almost empty. In order to start a calculation, one
will have to do something like::
atoms = Atoms(...)
atoms.calc = calc
atoms.get_potential_energy()
ASE will then arrange to call the :meth:`~gpaw.calculator.GPAW.calculate`
method with the correct arguments. This will trigger:
......@@ -89,7 +89,7 @@ method with the correct arguments. This will trigger:
a) Pass on the atomic positions to the wave functions, hamiltonian
and density objects (call their ``set_positions()`` methods).
b) Make sure the wave functions are initialized.
c) Reset the :class:`~gpaw.scf.SCFLoop`.
......@@ -171,7 +171,7 @@ Attributes of the wave function object: ``gd``, ``nspins``,
``ibzk_kc``, ``weight_k``, ``symmetry``, ``kpt_comm``,
``rank_a``, ``nibzkpts``, ``kpt_u``, ``setups``,
``ibzk_qc``, ``eigensolver``, and ``timer``.
.. _overview_xc:
......@@ -345,14 +345,14 @@ To investigate the way GPAW distributes calculated quantities across the
various MPI groups, simulating an MPI run can be done using ``gpaw-mpisim``::
$ gpaw-mpisim -v --dry-run=4 --spins=2 --kpoints=4 --bands=3 --domain-decomposition=2,1,1
Simulating: world.size = 4
parsize_c = (2, 1, 1)
parsize_bands = 1
nspins = 2
nibzkpts = 4
nbands = 3
world: rank=0, ranks=None
kpt_comm : rank=0, ranks=[0 2], mynks=4, kpt_u=[0^,1^,2^,3^]
band_comm : rank=0, ranks=[0], mynbands=3, mybands=[0, 1, 2]
......@@ -378,4 +378,5 @@ which in this case is equal to *mpi.world*.
.. [1] J J. Mortensen and L. B. Hansen and K. W. Jacobsen,
Phys. Rev. B 71 (2005) 035109.
.. [2] C. Rostgaard, `The Projector Augmented Wave Method <../paw_note.pdf>`_.
.. [2] C. Rostgaard, :download:`The Projector Augmented Wave Method
<../documentation/paw_note.pdf>`.
......@@ -11,34 +11,49 @@ small and quick tests and by a weekly set of larger test.
Quick test suite
================
Use the :program:`gpaw-test` command to run the tests::
Use the :program:`gpaw` command to run the tests::
$ gpaw-test --help
Usage: gpaw-test [options] [tests]
Options:
--version show program's version number and exit
$ gpaw test --help
usage: gpaw test [-h] [-x test1.py,test2.py,...] [-f] [--from TESTFILE]
[--after TESTFILE] [--range test_i.py,test_j.py] [-j JOBS]
[--reverse] [-k] [-d DIRECTORY] [-s] [--list]
[tests [tests ...]]
Run the GPAW test suite. The test suite can be run in parallel with MPI
through gpaw-python. The test suite supports 1, 2, 4 or 8 CPUs although some
tests are skipped for some parallelizations. If no TESTs are given, run all
tests supporting the parallelization.
positional arguments:
tests
optional arguments:
-h, --help show this help message and exit
-x test1.py,test2.py,..., --exclude=test1.py,test2.py,...
-x test1.py,test2.py,..., --exclude test1.py,test2.py,...
Exclude tests (comma separated list of tests).
-f, --run-failed-tests-only
Run failed tests only.
--from=TESTFILE Run remaining tests, starting from TESTFILE
--after=TESTFILE Run remaining tests, starting after TESTFILE
--range=test_i.py,test_j.py
--from TESTFILE Run remaining tests, starting from TESTFILE
--after TESTFILE Run remaining tests, starting after TESTFILE
--range test_i.py,test_j.py
Run tests in range test_i.py to test_j.py (inclusive)
-j JOBS, --jobs=JOBS Run JOBS threads.
-j JOBS, --jobs JOBS Run JOBS threads. Each test will be executed in serial
by one thread. This option cannot be used for
parallelization together with MPI.
--reverse Run tests in reverse order (less overhead with
multiple jobs)
-k, --keep-temp-dir Do not delete temporary files.
-d DIRECTORY, --directory DIRECTORY
Run test in this directory
-s, --show-output Show standard output from tests.
--list list the full list of tests, then exit
A temporary directory will be made and the tests will run in that
directory. If all tests pass, the directory is removed.
The test suite consists of a large number of small and quick tests
found in the :git:`gpaw/test/` directory. The tests run nightly in serial
and in parallel. Here are the results from `BuildBot
<https://gpaw-buildbot.fysik.dtu.dk/waterfall>`_.
and in parallel.
......@@ -66,7 +81,7 @@ Use this function to check results:
.. _big-test:
.. _agts:
Big tests
=========
......
......@@ -9,8 +9,8 @@ at the :ref:`tutorials` and consulting the :ref:`manual` pages.
If you are just looking for an introduction to the theory of PAW, you
should have a look at the page :ref:`introduction_to_paw`, read this
`pdf file on PAW theory`_, or maybe read some of the original
:ref:`PAW_papers`.
:download:`pdf file on PAW theory <paw_note.pdf>`, or maybe read some
of the original :ref:`PAW_papers`.
If you want some general information on the numerical techniques
utilized by the GPAW code, you should read the page
......@@ -30,8 +30,7 @@ Developers should consult the :ref:`devel` pages.
advanced_topics
literature
introduction_to_paw
utilities/utilities
If you can not find what you are looking for in any of the above,
please consult the :ref:`mail list`.
.. _pdf file on PAW theory: paw_note.pdf
from ase.build import molecule
from gpaw import GPAW
from gpaw import dscf
# Ground state calculation
calc = GPAW(nbands=8,
h=0.2,
xc='PBE',
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO = molecule('CO')
CO.center(vacuum=3)
CO.set_calculator(calc)
E_gs = CO.get_potential_energy()
# Obtain the pseudowavefunctions and projector overlaps of the
# state which is to be occupied. n=5,6 is the 2pix and 2piy orbitals
n = 5
molecule = [0, 1]
wf_u = [kpt.psit_nG[n] for kpt in calc.wfs.kpt_u]
p_uai = [dict([(molecule[a], P_ni[n]) for a, P_ni in kpt.P_ani.items()])
for kpt in calc.wfs.kpt_u]
# Excited state calculation
calc_es = GPAW(nbands=8,
h=0.2,
xc='PBE',
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO.set_calculator(calc_es)
lumo = dscf.AEOrbital(calc_es, wf_u, p_uai)
# lumo = dscf.MolecularOrbital(calc, weights={0: [0, 0, 0, 1],
# 1: [0, 0, 0, -1]})
dscf.dscf_calculation(calc_es, [[1.0, lumo, 1]], CO)
E_es = CO.get_potential_energy()
print('Excitation energy: ', E_es - E_gs)
......@@ -8,9 +8,9 @@ Delta Self-Consistent Field
Linear expansion Delta Self-Consistent Field
--------------------------------------------
The method of linear expansion Delta Self-Consistent Field \ [#delscf]_
adds the density of a specified orbital `\varphi_a(r)` to the
total density in each step of the self-consistency cycle. The extra charge
The method of linear expansion Delta Self-Consistent Field \ [#delscf]_
adds the density of a specified orbital `\varphi_a(r)` to the
total density in each step of the self-consistency cycle. The extra charge
is usually taken from the fermi level to keep the system neutral:
.. math::
......@@ -18,12 +18,12 @@ is usually taken from the fermi level to keep the system neutral:
n(r) = \sum_nf_{N-1}(T,\varepsilon_n)|\varphi_n(r)|^2 + |\varphi_a(r)|^2.
with `N` being the total number of electrons and `f_{N-1}(T,\varepsilon_n)`
is the Fermi-Dirac distribution of the `N-1` electron system . To get the
is the Fermi-Dirac distribution of the `N-1` electron system . To get the
band energy right `\varphi_a(r)` needs to be expanded in Kohn-Sham orbitals:
.. math::
|\varphi_a\rangle = \sum_nc_{na}|\varphi_n\rangle,
|\varphi_a\rangle = \sum_nc_{na}|\varphi_n\rangle,
\qquad c_{na} = \langle\varphi_n|\varphi_a\rangle
and the band energy of the orbital becomes
......@@ -33,119 +33,75 @@ and the band energy of the orbital becomes
\varepsilon_a = \sum_n|c_{na}|^2\varepsilon_n.
The method is a generalization of traditional Delta Self-Consistent Field
where only the occupation numbers are modified and it will reduce to that,
where only the occupation numbers are modified and it will reduce to that,
if only one (normalized) term is included in the expansion of `\varphi_a(r)`.
----------------
Simple molecules
----------------
The example below calculates the excitation energy of the
`5\sigma\rightarrow2\pi` transition in CO. We only specify that the
`2\pi` orbital should be occupied ([[1.0, lumo, 1]] means 1.0 electrons
in lumo with spin 1) and the method will take the electron from highest
The example below calculates the excitation energy of the
`5\sigma\rightarrow2\pi` transition in CO. We only specify that the
`2\pi` orbital should be occupied ([[1.0, lumo, 1]] means 1.0 electrons
in lumo with spin 1) and the method will take the electron from highest
occupied orbital which in this case is `5\sigma`.
The lumo is an instance of the class AEOrbital which calculates the
The lumo is an instance of the class AEOrbital which calculates the
expansion of the saved `2\pi` state in each iteration step.
In order to obtain the all-electron overlaps `\langle\varphi_n|2\pi\rangle`
we need to supply the projector overlaps in addition to the
In order to obtain the all-electron overlaps `\langle\varphi_n|2\pi\rangle`
we need to supply the projector overlaps in addition to the
pseudowavefunction.
Exciting the LUMO in CO::
Exciting the LUMO in CO (:download:`co.py`):
from ase.build import molecule
from gpaw import GPAW
from gpaw import dscf
.. literalinclude:: co.py
# Ground state calculation
#------------------------------------------------------------------
calc = GPAW(nbands=8, h=0.2, xc='PBE', spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO = molecule('CO')
CO.center(vacuum=3)
CO.set_calculator(calc)
E_gs = CO.get_potential_energy()
# Obtain the pseudowavefunctions and projector overlaps of the
# state which is to be occupied. n=5,6 is the 2pix and 2piy orbitals
n = 5
molecule = [0, 1]
wf_u = [kpt.psit_nG[n] for kpt in calc.wfs.kpt_u]
p_uai = [dict([(molecule[a], P_ni[n]) for a, P_ni in kpt.P_ani.items()])
for kpt in calc.wfs.kpt_u]
# Excited state calculation
#--------------------------------------------
calc_es = GPAW(nbands=8, h=0.2, xc='PBE', spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO.set_calculator(calc_es)
lumo = dscf.AEOrbital(calc_es, wf_u, p_uai)
#lumo = dscf.MolecularOrbital(calc, weights={0: [0, 0, 0, 1],
# 1: [0, 0, 0, -1]})
dscf.dscf_calculation(calc_es, [[1.0, lumo, 1]], CO)
E_es = CO.get_potential_energy()
print 'Excitation energy: ', E_es-E_gs
The commented line ``lumo = dscf.Molecular...``
uses another class to specify the `2\pi` orbital of CO which does not require
a ground state calculation of the molecule. In the simple example above the
two methods give identical results, but for more complicated systems the
AEOrbital class should be used \ [#des]_. When using the AEOrbital class
The commented lines ``lumo = dscf.Molecular...``
uses another class to specify the `2\pi` orbital of CO which does not require
a ground state calculation of the molecule. In the simple example above the
two methods give identical results, but for more complicated systems the
AEOrbital class should be used \ [#des]_. When using the AEOrbital class
a new calculator object must be constructed for the dscf calculation.
In the example above we only specify a single state, but the function
``dscf.dscf_calculation`` takes a list of orbitals as input and we could for
example have given the argument [[1.0, lumo, 1], [-1.0, pi, 0]] which would
force the electron to be taken from the `\pi` orbital with spin 0. The pi
In the example above we only specify a single state, but the function
``dscf.dscf_calculation`` takes a list of orbitals as input and we could for
example have given the argument [[1.0, lumo, 1], [-1.0, pi, 0]] which would
force the electron to be taken from the `\pi` orbital with spin 0. The pi
should of course be another instance of the AEOrbital class.
---------------------
Exciting an adsorbate
---------------------
The method of linear expansion Delta Self-Consistent Field was designed
for calculations with strongly hybridized orbitals. For example molecules
chemisorbed on transition metals. In such cases the
for calculations with strongly hybridized orbitals. For example molecules
chemisorbed on transition metals. In such cases the
traditional Delta Self-Consistent Field breaks down since the orbital
to be occupied is no longer well described by a single Kohn-Sham state.
The script :git:`~doc/documentation/dscf/homo.py` calculates
The script :git:`~doc/documentation/dscf/homo.py` calculates
the HOMO energy of CO adsorbed on-top Pt(111). The script starts
from scratch, but usually one would start from an optimized configuration
saved in a file ``gs.gpw``. The script only calculates the total energy of
the excited state so the excitation energy is obtained as the difference
saved in a file ``gs.gpw``. The script only calculates the total energy of
the excited state so the excitation energy is obtained as the difference
between ground and excited state energies.
First a calculation of gas-phase CO is performed and the
HOMO pseudo-wavefunctions and the projector overlaps are saved. The
First a calculation of gas-phase CO is performed and the
HOMO pseudo-wavefunctions and the projector overlaps are saved. The
energy range [-100.0, 0.0] means we only include states below the Fermi
level (default is states above).
The script :git:`~doc/documentation/dscf/lumo.py` calculates
the LUMO energy of the same system, but is slightly more complicated due to
the degeneracy of the `2\pi` orbital. We would like to occupy the `2\pi_y`
orbital and we need to figure out which band (5 or 6) this orbital
the LUMO energy of the same system, but is slightly more complicated due to
the degeneracy of the `2\pi` orbital. We would like to occupy the `2\pi_y`
orbital and we need to figure out which band (5 or 6) this orbital
corresponds to in each k-point before we start the slab calculation.
.. [#delscf] J. Gavnholt, T. Olsen, M. Engelund and J. Schiøtz,
Delta Self-Consistent Field as a method to obtain potential
energy surfaces of excited molecules on surfaces,
energy surfaces of excited molecules on surfaces,
*Phys. Rev. B* **78**, 075441 (2008)
.. [#des] T. Olsen, J. Gavnholt and J. Schiøtz,
Hot electron mediated desorption rates calculated from excited
state potential energy surfaces,
state potential energy surfaces,
*Phys. Rev. B* **79**, 035403 (2009)
#!/usr/bin/env python
from ase.visualize import view
from ase.build import fcc111, add_adsorbate
from gpaw import GPAW
from gpaw.mixer import MixerSum
import gpaw.dscf as dscf
filename='homo'
#-------------------------------------------
filename = 'homo'
c_mol = GPAW(nbands=9, h=0.2, xc='RPBE', kpts=(8,6,1),
c_mol = GPAW(nbands=9,