Commit 37366422 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'master' into exx-rewrite

parents 015a7dee ed9c0e1a
Pipeline #99932964 passed with stage
in 3 minutes and 27 seconds
......@@ -7,10 +7,12 @@ master:
- pip install git+https://gitlab.com/ase/[email protected] > ase-install.output
- mkdir ~/.gpaw
- echo "scalapack = True; libraries += ['scalapack-openmpi']" > ~/.gpaw/siteconfig.py
- python setup.py install # > gpaw-install.output
- python setup.py install > gpaw-install.output
- useradd -m user
- su user -c 'gpaw install-data --register gpaw-datasets'
- su user -c 'gpaw -P 1 info'
- su user -c 'gpaw test --range linalg/gemm_complex.py,lcao/dos.py'
- flake8 --ignore E722,E222,E701,E731,E114,E402,E502,E262,E266,E203,E305,E128,E221,E127,E202,E303,E302,E201,E241,E225,E251,E261,E265,E231,E226,W291,E501,E129,E741,W293,W504,W503,W605 --exclude "doc/platforms/*" gpaw doc
- flake8 --exit-zero --exclude "doc/platforms/*,gpaw/lrtddft2/*" gpaw doc > f8.out
- cat f8.out | wc -l | tee nerr
- python -c "assert int(open('nerr').read()) <= 2364, 'Please run flake8 on your code'"
- python -We:invalid -m compileall -f -q gpaw/
This diff is collapsed.
This diff is collapsed.
include MANIFEST.in
include COPYING LICENSE CONTRIBUTING.rst CHANGELOG.rst
include config.py siteconfig.py
include LICENSE CONTRIBUTING.rst CHANGELOG.rst
include config.py siteconfig_example.py
include c/*.c
include c/*.h
include c/xc/*.h
......
......@@ -12,7 +12,7 @@
#ifdef __clang__
#define INLINE static
#else
#define INLINE inline
#define INLINE static inline
#endif
// Index notation
......@@ -230,7 +230,7 @@ PyObject* adjust_positions(PyObject *self, PyObject *args)
PyArrayObject* arrayR_niv = 0; // Output: Adjusted positions will be written here.
PyArrayObject* arraynewR_niv = 0; // Input: gives the positions to be adjusted.
if (!PyArg_ParseTuple(args, "OOOO", &arraylen_x, &arraymass_i,
if (!PyArg_ParseTuple(args, "OOOO", &arraylen_x, &arraymass_i,
&arrayR_niv, &arraynewR_niv))
{
return NULL;
......@@ -368,9 +368,9 @@ PyObject* adjust_momenta(PyObject *self, PyObject *args)
{
PyArrayObject* arraymass_i = 0; // Input: the 3 masses
PyArrayObject* arrayR_niv = 0; //
PyArrayObject* arraynewP_niv = 0; //
PyArrayObject* arraynewP_niv = 0; //
if (!PyArg_ParseTuple(args, "OOO", &arraymass_i, &arrayR_niv,
if (!PyArg_ParseTuple(args, "OOO", &arraymass_i, &arrayR_niv,
&arraynewP_niv))
{
return NULL;
......@@ -459,7 +459,7 @@ PyObject* adjust_momenta(PyObject *self, PyObject *args)
vec9_dot(denom_x, d_xv, d_xv); // (R1-R2)^2 etc.
vec3_div(lambda_x, g_x, denom_x); // g_12 / (R1-R2)^2
vec3_imul(lambda_x, mu_x); // lambda /= (m_1^-1 + m_2^-1) etc.
#ifdef FF_DEBUG_M
vec3_print("lambda_x", lambda_x);
#endif
......@@ -496,7 +496,7 @@ PyObject* calculate_forces_H2O(PyObject *self, PyObject *args)
double cutoff=0;
double width=0;
if (!PyArg_ParseTuple(args, "OOddddOOO", &arraypbc, &arraycell, &A, &B, &cutoff, &width, &arrayZ_i, &arrayR_niv, &arrayF_niv))
{
return NULL;
......@@ -514,7 +514,7 @@ PyObject* calculate_forces_H2O(PyObject *self, PyObject *args)
return NULL;
}
if (!(PyArray_NDIM(arraycell) == 2 &&
if (!(PyArray_NDIM(arraycell) == 2 &&
(PyArray_DIM(arraycell,0) == 3) &&
(PyArray_DIM(arraycell,1) == 3))) {
PyErr_SetString(PyExc_TypeError, "Cell should be array with size 3x3.");
......@@ -523,10 +523,10 @@ PyObject* calculate_forces_H2O(PyObject *self, PyObject *args)
double* cell_vc = DOUBLEP(arraycell);
//printf("Cell\n");
for (unsigned int v1=0; v1<3; v1++)
for (unsigned int v1=0; v1<3; v1++)
for (unsigned int v2=0; v2<3; v2++) {
//printf("%20.15f \n", cell_vc[v1+v2*3]);
if (v1 != v2)
if (v1 != v2)
if (fabs(cell_vc[v1+v2*3]) > 1e-10) {
PyErr_SetString(PyExc_TypeError, "Cell array should be diagonal.");
return NULL;
......@@ -556,7 +556,7 @@ PyObject* calculate_forces_H2O(PyObject *self, PyObject *args)
double* F_niv = DOUBLEP(arrayF_niv);
double E = 0.0;
if (cutoff <= 0.0) { // No cutoff
if (cutoff <= 0.0) { // No cutoff
// Double loop of atoms
for (unsigned int n1 = 0; n1 < NW; n1++)
for (unsigned int n2 = n1+1; n2 < NW; n2++) {
......
......@@ -17,7 +17,7 @@ class Box:
self.owns = []
self.position = None
def set_position(self, position):
self.position = np.asarray(position)
......@@ -46,7 +46,7 @@ class MPL:
def plot(self):
a4 = 100 * np.array([2**-1.75, 2**-2.25])
inch = 2.54
self.fig = plt.figure(1, a4 / inch)
self.ax = ax = self.fig.add_axes([0, 0, 1, 1], frameon=False)
ax.set_xlim(0, a4[0])
......@@ -84,7 +84,7 @@ class MPL:
bbox = b.text.get_window_extent()
t = b.text.get_transform()
b.size = t.inverted().transform(bbox.size)
for b in self.boxes:
for other, name, s, style in b.owns:
d = other.position - b.position
......@@ -116,6 +116,7 @@ def box(*args, **kwargs):
boxes.append(b)
return b
atoms = box('Atoms', [''], ['positions, numbers, cell, pbc'],
color='white')
paw = box('PAW', [], [], 'green')
......@@ -153,7 +154,10 @@ fd = box('FDWaveFunctions',
[], 'magenta')
pt = box('LFC', r'$\tilde{p}_i^a(\mathbf{r})$', [], 'red')
lcao = box('LCAOWaveFunctions',
r"$\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})=\sum_{\mu\mathbf{R}} C_{\sigma\mathbf{k}n\mu} \Phi_\mu(\mathbf{r} - \mathbf{R}) \exp(i\mathbf{k}\cdot\mathbf{R})$",
r"$\tilde{\psi}_{\sigma\mathbf{k}n}(\mathbf{r})=" +
r"\sum_{\mu\mathbf{R}} " +
r"C_{\sigma\mathbf{k}n\mu} \Phi_\mu(\mathbf{r} - \mathbf{R}) " +
r"\exp(i\mathbf{k}\cdot\mathbf{R})$",
['S_qMM, T_qMM, P_aqMi'], 'magenta')
atoms0 = box('Atoms', '(copy)', ['positions, numbers, cell, pbc'],
color='grey')
......@@ -162,8 +166,9 @@ forces = box('ForceCalculator')
occupations = box(
'OccupationNumbers',
r'$\epsilon_{\sigma\mathbf{k}n} \rightarrow f_{\sigma\mathbf{k}n}$')
poisson = box('PoissonSolver',
r'$\nabla^2 \tilde{v}_H(\mathbf{r}) = -4\pi \tilde{\rho}(\mathbf{r})$')
poisson = box(
'PoissonSolver',
r'$\nabla^2 \tilde{v}_H(\mathbf{r}) = -4\pi \tilde{\rho}(\mathbf{r})$')
eigensolver = box('EigenSolver')
symmetry = box('Symmetry')
restrictor = box('Transformer', '(fine -> coarse)',
......@@ -175,7 +180,7 @@ kin = box('FDOperator', r'$-\frac{1}{2}\nabla^2$')
hsoperator = box('HSOperator',
[r"$\langle \psi_n | A | \psi_{n'} \rangle$",
r"$\sum_{n'}U_{nn'}|\tilde{\psi}_{n'}\rangle$"])
overlap = box('Overlap')
basisfunctions = box('BasisFunctions', r'$\Phi_\mu(\mathbf{r})$',
color='red')
......
# creates: lines.png
from __future__ import division
import datetime as dt
import os
import subprocess
......@@ -65,13 +64,13 @@ def count_lines():
print(year, month, hash)
subprocess.call(['git', 'checkout', hash])
c = count('c', '\*.[ch]')
py = count('.', '\*.py')
test = count('gpaw/test', '\*.py')
test += count('test', '\*.py')
doc = count('doc', '\*.py')
c = count('c', r'\*.[ch]')
py = count('.', r'\*.py')
test = count('gpaw/test', r'\*.py')
test += count('test', r'\*.py')
doc = count('doc', r'\*.py')
py -= test + doc # avoid double counting
rst = count('.', '\*.rst')
rst = count('.', r'\*.rst')
print(year, month, 0, c, py, test, doc, rst, file=fd)
month += 1
if month == 13:
......
......@@ -6,6 +6,7 @@ machine = os.environ.get('MACHINE', 'TEST')
ncores = os.environ.get('NCORES', 8)
if rank == 0:
os.chdir(machine)
os.system('python ../memory_bandwidth.py --runs=5 --startcores='+str(ncores))
#os.system('python ../memory_bandwidth.py --runs=5') # full memory benchmark
shutil.copy('memory_bandwidth_'+machine+'_py.png', '..')
os.system(
'python ../memory_bandwidth.py --runs=5 --startcores=' +
str(ncores))
shutil.copy('memory_bandwidth_' + machine + '_py.png', '..')
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import print_function
import Numeric as num
from gpaw import GPAW
from gpaw.spherical_harmonics import Y
......@@ -60,10 +59,12 @@ for n in paw.nuclei:
p.plot(X + x, num.dot(P_i, phi_i), C + '-', lw=2,
label=r'$\psi^%s$' % s.symbol)
C = 'r'
p.plot([-d / 2], [0], 'go', ms=2*0.8*4*80/a*1.0*0.53, mfc=None, label='_nolegend_')
p.plot([d / 2], [0], 'ro', ms=2*0.8*4*80/a*1.2*0.53, mfc=None, label='_nolegend_')
p.plot([-d / 2], [0], 'go', ms=2*0.8*4*80/a*1.0*0.53, mfc=None,
label='_nolegend_')
p.plot([d / 2], [0], 'ro', ms=2*0.8*4*80/a*1.2*0.53, mfc=None,
label='_nolegend_')
p.legend(loc='best')
p.xlabel(u'x [Å]')
p.ylabel(r'$\psi$')
#p.show()
# p.show()
p.savefig('co_wavefunctions.png', dpi=dpi)
# Creates: cl_field.ind_Ffe.png, qm_field.ind_Ffe.png, tot_field.ind_Ffe.png
# -*- coding: utf-8 -*-
from gpaw.mpi import world
assert world.size == 1, 'This script should be run in serial mode (with one process).'
import numpy as np
import matplotlib.pyplot as plt
......@@ -9,6 +8,8 @@ import matplotlib.pyplot as plt
from gpaw.inducedfield.inducedfield_base import BaseInducedField
from gpaw.tddft.units import aufrequency_to_eV
assert world.size == 1, 'This script should be run in serial mode.'
# Helper function
def do_plot(d_g, ng, box, atoms):
......@@ -38,6 +39,7 @@ def do_plot(d_g, ng, box, atoms):
plt.ylim([y[0], y[-1]])
ax.set_aspect('equal')
for fname, name in zip(['cl_field.ind', 'qm_field.ind', 'tot_field.ind'],
['Classical subsystem', 'Quantum subsystem',
'Total hybrid system']):
......
# Creates: field.ind_Ffe.png
# -*- coding: utf-8 -*-
from gpaw.mpi import world
assert world.size == 1, 'This script should be run in serial mode (with one process).'
import numpy as np
import matplotlib.pyplot as plt
......@@ -9,6 +8,8 @@ import matplotlib.pyplot as plt
from gpaw.inducedfield.inducedfield_base import BaseInducedField
from gpaw.tddft.units import aufrequency_to_eV
assert world.size == 1, 'This script should be run in serial mode.'
# Helper function
def do_plot(d_g, ng, box, atoms):
......@@ -38,6 +39,7 @@ def do_plot(d_g, ng, box, atoms):
plt.ylim([y[0], y[-1]])
ax.set_aspect('equal')
for fname, name in zip(['field.ind'], ['Classical system']):
# Read InducedField object
ind = BaseInducedField(fname, readmode='all')
......
......@@ -36,7 +36,7 @@ def plot(fname, fiteps):
# plt.ylim(fiteps_e.real.min(), fiteps_e.real.max())
plt.ylim(-70, 0)
plt.xlabel('Energy (eV)')
plt.ylabel('Real($\epsilon$)')
plt.ylabel(r'Real($\epsilon$)')
plt.legend(loc='best')
plt.subplot(1, 2, 2)
......@@ -46,7 +46,7 @@ def plot(fname, fiteps):
# plt.ylim(fiteps_e.imag.min(), fiteps_e.imag.max())
plt.ylim(0, 7)
plt.xlabel('Energy (eV)')
plt.ylabel('Imaginary($\epsilon$)')
plt.ylabel(r'Imaginary($\epsilon$)')
plt.tight_layout()
plt.savefig('%s.png' % fname)
......@@ -55,6 +55,8 @@ def plot(fname, fiteps):
# http://refractiveindex.info/?shelf=main&book=Au&page=Johnson
# Direct download link:
# wget http://refractiveindex.info/database/main/Au/Johnson.yml -O Au.yml
ymlfname = 'Au.yml'
# Fit to the permittivity
......
......@@ -48,4 +48,6 @@ f2 = lr.lr_transitions.get_transition_contributions(0)
for (ip, val) in enumerate(f2):
if (val > 1e-3):
parprint(' %5d => %5d %lf %%\n ' %
(lr.ks_singles.kss_list[ip].occ_ind, lr.ks_singles.kss_list[ip].unocc_ind, val / 2. * 100))
(lr.ks_singles.kss_list[ip].occ_ind,
lr.ks_singles.kss_list[ip].unocc_ind,
val / 2. * 100))
......@@ -7,7 +7,9 @@ manually instead of using gpawtransport, which currently does not work
from ase import Atoms
from gpaw import GPAW, Mixer, FermiDirac
from gpaw.lcao.tools import remove_pbc, get_lcao_hamiltonian, get_lead_lcao_hamiltonian
from gpaw.lcao.tools import (remove_pbc,
get_lcao_hamiltonian,
get_lead_lcao_hamiltonian)
import pickle as pickle
a = 2.41 # Pt binding length
......@@ -48,7 +50,7 @@ H -= Ef * S
remove_pbc(atoms, H, S, 0)
# Dump the Hamiltonian and Scattering matrix to a pickle file
pickle.dump((H.astype(complex), S.astype(complex)),
pickle.dump((H.astype(complex), S.astype(complex)),
open('scat_hs.pickle', 'wb'), 2)
########################
......
......@@ -54,7 +54,7 @@ the directory where you want to install
mkdir -p $MYLIBXCDIR
cd $MYLIBXCDIR
wget http://www.tddft.org/programs/octopus/down.php?file=libxc/4.3.4/libxc-4.3.4.tar.gz -O libxc-4.3.4.tar.gz
wget http://www.tddft.org/programs/libxc/down.php?file=4.3.4/libxc-4.3.4.tar.gz -O libxc-4.3.4.tar.gz
tar xvzf libxc-4.3.4.tar.gz
cd libxc-4.3.4
mkdir install
......@@ -90,7 +90,7 @@ You might want to install a stable version of ASE::
ASE_SOURCE=$PWD/source/ase
mkdir -p $ASE_SOURCE
cd $ASE_SOURCE
git clone -b 3.17.0 https://gitlab.com/ase/ase.git 3.17.0
git clone -b 3.18.1 https://gitlab.com/ase/ase.git 3.18.1
We add our installation to the module environment::
......@@ -98,7 +98,7 @@ We add our installation to the module environment::
mkdir -p modulefiles/ase
cd modulefiles/ase
Edit the module file :file:`3.17.0` that should read::
Edit the module file :file:`3.18.1` that should read::
#%Module1.0
......@@ -166,7 +166,12 @@ A specific tag can be loaded by::
# load version 1.2.0
git checkout 1.2.0
To build GPAW use::
To build the current trunk version of GPAW we need to create
a file :file:`setup.py` that reads
.. literalinclude:: nemo_siteconfig.py
Then we build the executable::
module purge
module load libxc
......@@ -175,17 +180,16 @@ To build GPAW use::
module load ase
cd $GPAW_SOURCE/trunk
mkdir install
python3 setup.py install --prefix=$PWD/install
python3 setup.py build
which installs GPAW to ``$GPAW_SOURCE/trunk/install``.
We create a module that does the necessary things::
which builds GPAW to ``$GPAW_SOURCE/trunk/build``.
We create a module that creates the necessary definitions::
cd
mkdir -p modulefiles/gpaw
cd modulefiles/gpaw
the file :file:`trunk` that should read::
The file :file:`trunk` that should read::
#%Module1.0
......@@ -195,11 +199,13 @@ the file :file:`trunk` that should read::
if {![is-loaded numlib/mkl]} {module load numlib/mkl}
if {![is-loaded gpaw-setups]} {module load gpaw-setups}
# change this to your needs
set gpawhome /home/fr/fr_fr/fr_mw767/source/gpaw/trunk/install
prepend-path PATH $gpawhome/bin
prepend-path PYTHONPATH $gpawhome/lib/python3.6/site-packages/
setenv GPAW_PYTHON $gpawhome/bin/gpaw-python
# change the following directory definition to your needs
set gpawhome /home/fr/fr_fr/fr_mw767/source/gpaw/trunk
# this can stay as is
prepend-path PATH $gpawhome/build/scripts-3.6
prepend-path PYTHONPATH $gpawhome
setenv GPAW_PYTHON $gpawhome/build/bin.linux-x86_64-3.6/gpaw-python
Running GPAW
------------
......
parallel_python_interpreter = True
libraries += ['mkl_intel_lp64', 'mkl_sequential', 'mkl_core']
......@@ -99,7 +99,7 @@ def bondlengths(Ea, dE):
plt.plot(B, E0, 'g.', label='reference')
plt.legend(loc='lower right')
plt.xlabel('Bond length $\mathrm{\AA}$')
plt.xlabel(r'Bond length $\mathrm{\AA}$')
plt.ylabel('Energy [eV]')
plt.savefig('bondlengths.png')
......
......@@ -6,7 +6,8 @@ Lattice="4.12605070262 0.0 0.0 -2.06302535131 3.57326472577 0.0 0.0 0.0 18.13127
V 0.00000000 -0.00000000 9.06563772 1.00000000 23 0.00000000 0.00000000 0.00000003 2.63358805 -0.00000000
I 2.06302535 1.19108824 10.63303664 0.10000000 53 0.00000000 0.00000000 0.00792029 -0.03990743 0.00000000
I -0.00000000 2.38217648 7.49823879 0.10000000 53 -0.00000000 0.00000000 -0.00792046 -0.03990743 -0.00000000
""")
""") # noqa
Path('CrI3.xyz').write_text("""\
8
Lattice="7.00794410138 5.69308307302e-19 0.0 -3.50397205069 6.0690576201 0.0 -7.55612427101e-19 0.0 18.0635365764" Properties=species:S:1:pos:R:3:initial_magmoms:R:1:Z:I:1:forces:R:3:magmoms:R:1 energy=-31.3985425600531 dipole="-0.00239032749259 0.00138005621976 3.42821467339e-07" free_energy=-31.3985446012773 pbc="T T F" magmom=6.00000000876689 stress="-6.01056204961e-05 0.0 -6.48091502334e-21 0.0 -6.01056204961e-05 -2.59236600933e-20 -6.48091502334e-21 -2.59236600933e-20 6.05488486747e-05"
......@@ -18,4 +19,4 @@ I -2.24281648 3.88469895 7.46488009 1.00000000 53
I 2.52228800 -0.00001162 10.59865633 1.00000000 53 -0.00469075 0.00014771 0.00214167 -0.05825936
I 2.24281799 3.88469794 10.59865633 1.00000000 53 0.00247329 0.00398846 0.00214167 -0.05825936
I -1.26113394 2.18437130 10.59865633 1.00000000 53 0.00221746 -0.00413616 0.00214167 -0.05825936
""")
""") # noqa
......@@ -7,12 +7,13 @@ from gpaw.solvation import (
SolvationGPAW, # the solvation calculator
EffectivePotentialCavity, # cavity using an effective potential
Power12Potential, # a specific effective potential
LinearDielectric, # rule to construct permittivity function from the cavity
GradientSurface, # rule to calculate the surface area from the cavity
SurfaceInteraction # rule to calculate non-electrostatic interactions
LinearDielectric, # rule to construct permittivity func from the cavity
GradientSurface, # rule to calculate the surface area from the cavity
SurfaceInteraction # rule to calculate non-electrostatic interactions
)
# all parameters on the user side of the solvation API follow the ASE unit conventions (eV, Angstrom, ...)
# all parameters on the user side of the solvation API follow the ASE
# unit conventions (eV, Angstrom, ...)
# non-solvent related DFT parameters
h = 0.2
......@@ -21,14 +22,20 @@ vac = 5.0
# solvent parameters for water from J. Chem. Phys. 141, 174108 (2014)
u0 = 0.180 # eV
epsinf = 78.36 # dimensionless
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom ** 2
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom**2
T = 298.15 # Kelvin
vdw_radii = vdw_radii.copy()
vdw_radii[1] = 1.09
# atomic_radii expected by gpaw.solvation.Power12Potential have to be a callable
# mapping an Atoms object to an iterable of floats representing the atomic radii of
# the atoms in the same order as in the Atoms object (in Angstrom)
atomic_radii = lambda atoms: [vdw_radii[n] for n in atoms.numbers]
# atomic_radii expected by gpaw.solvation.Power12Potential have to be a
# callable mapping an Atoms object to an iterable of floats representing
# the atomic radii of the atoms in the same order as in the Atoms object
# (in Angstrom)
def atomic_radii(atoms):
return [vdw_radii[n] for n in atoms.numbers]
# create Atoms object for ethanol and add vacuum
atoms = molecule('CH3CH2OH')
......@@ -38,7 +45,8 @@ atoms.center(vacuum=vac)
atoms.calc = GPAW(xc='PBE', h=h, txt='gasphase.txt')
Egasphase = atoms.get_potential_energy()
# perform calculation with continuum solvent model from J. Chem. Phys. 141, 174108 (2014)
# perform calculation with continuum solvent model from
# J. Chem. Phys. 141, 174108 (2014)
atoms.calc = SolvationGPAW(
xc='PBE', h=h, txt='water.txt',
cavity=EffectivePotentialCavity(
......
......@@ -464,6 +464,26 @@ There are few points about the implementation that we emphasize:
dielectric function converges slowly with the upper bound of the frequency
grid. Refer to :ref:`df_theory` for the details on the Hilbert transform.
Drude phenomenological scattering rate
======================================
A phenomenological scattering rate can be introduced using the ``rate``
parameter in the optical limit. By default this is zero but can be set by
using::
df = DielectricFunction(...
rate=0.1,
...)
to set a scattering rate of 0.1 eV. If rate='eta' then the code with use the
specified ``eta`` parameter. Note that the implementation of the rate parameter
differs from some literature by a factor of 2 for consistency with the linear
response formalism. In practice the Drude rate is implemented as
.. math::
\frac{\omega_\mathrm{p}^2}{(\omega + i\gamma)^2}
where `\gamma` is the rate parameter.
Useful tips
===========
......@@ -525,6 +545,13 @@ keyword type default value description
``hilbert`` ``bool`` True Switch for hilbert transform.
``nbands`` ``int`` nbands from gs calc Number of bands from gs calc
to include.
``rate`` ``float`` or
``str`` 0.0 (eV) Phenomenological Drude
scattering rate. If rate="eta" then
use "eta". Note that this may differ
by a factor of 2 for some definitions
of the Drude scattering rate. See the
section on Drude scattering rate.
================= ================= =================== ================================
......
# Creates: silicon_ABS.png
import numpy as np
import matplotlib.pyplot as plt
plt.figure(figsize=(7, 5))
d = np.loadtxt('si_abs.csv', delimiter=',')
plt.plot(d[:, 0], d[:, 3], '-k', label='$\mathrm{Re}\epsilon(\omega)$')
plt.plot(d[:, 0], d[:, 4], '-r', label='$\mathrm{Im}\epsilon(\omega)$')
plt.plot(d[:, 0], d[:, 3], '-k', label=r'$\mathrm{Re}\epsilon(\omega)$')
plt.plot(d[:, 0], d[:, 4], '-r', label=r'$\mathrm{Im}\epsilon(\omega)$')
plt.title('Dielectric function of Si')
plt.legend()
plt.xlabel('Energy (eV)', fontsize=14)
plt.ylabel('$\epsilon$', fontsize=18)
plt.ylabel(r'$\epsilon$', fontsize=18)
ax = plt.gca()
......@@ -20,7 +20,7 @@ def y(e):
i = (x < e).sum()
return d[i, 4]
# data from G.Kresse, PRB 73, 045112 (2006)
for name, e in zip("E_0 E_1 E_2 E_0' E_1'".split(),
[2.53, 2.71, 3.72, 3.08, 4.50]):
......
......@@ -5,14 +5,14 @@ import matplotlib.pyplot as plt
plt.figure(figsize=(5, 7))
Q = np.loadtxt('graphite_q_list')
for i, q in enumerate(Q):
filename = 'graphite_EELS_' + str(i + 1)
filename = 'graphite_EELS_' + str(i + 1)
d = np.loadtxt(filename, delimiter=',')
plt.plot(d[:, 0], d[:, 2] + 0.4 * (4 - i), '-',
label='%.2f Ang$^{-1}$' % q)
plt.xlabel('Energy [eV]')
plt.ylabel('Loss Function')
plt.title('EELS spectra of graphene: $\Gamma-\mathrm{M}$')
plt.title(r'EELS spectra of graphene: $\Gamma-\mathrm{M}$')
plt.legend(loc='best')
plt.savefig('graphite_EELS.png')
plt.show()
......@@ -18,8 +18,8 @@ for omega2 in [2.5, 5, 10, 15, 20, np.inf]:
label = '$\\omega_2 = %.1f\\, \\mathrm{eV}$' % omega2
plt.plot(x, y, '.', label=label)
plt.ylabel('Freq. no')