Commit 155a2f19 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'master' into agts7

parents 19b80dd9 0e2844c9
Pipeline #20995045 passed with stage
in 3 minutes and 30 seconds
......@@ -15,6 +15,7 @@ Requirements
* Python_ 2.7, 3.4-
* ASE_ (atomic simulation environment)
* NumPy_ (base N-dimensional array package)
* SciPy_ (library for scientific computing)
* LibXC
* BLAS
* LAPACK
......@@ -23,7 +24,6 @@ Optional:
* MPI
* ScaLAPACK
* SciPy_ (library for scientific computing)
Installation
......
......@@ -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);
......@@ -184,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},
......
......@@ -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')
......
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)
......@@ -52,55 +52,11 @@ 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...``
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
......
#!/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,
h=0.2,
xc='RPBE',
kpts=(8, 6, 1),
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': 'occupied'}, txt='CO_homo.txt')
'bands': 'occupied'},
txt='CO_homo.txt')
calc = GPAW(nbands=45, h=0.2, xc='RPBE', kpts=(8,6,1),
calc = GPAW(nbands=80,
h=0.2,
xc='RPBE',
kpts=(8, 6, 1),
eigensolver='cg',
spinpol=True,
mixer=MixerSum(nmaxold=5, beta=0.1, weight=100),
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-7,
'bands': -10}, txt=filename+'.txt')
'bands': -10},
txt=filename + '.txt')
#----------------------------------------
# Import Slab with relaxed CO
#slab = ('gs.gpw').get_atoms()
# Import Slab with relaxed CO
slab = fcc111('Pt', size=(1, 2, 3), orthogonal=True)
add_adsorbate(slab, 'C', 2.0, 'ontop')
add_adsorbate(slab, 'O', 3.15, 'ontop')
......@@ -39,25 +40,22 @@ view(slab)
molecule = slab.copy()
del molecule [:-2]
del molecule[:-2]
# Molecule
#----------------
# Molecule
molecule.set_calculator(c_mol)
molecule.get_potential_energy()
#Homo wavefunction
# Homo wavefunction
wf_u = [kpt.psit_nG[4] for kpt in c_mol.wfs.kpt_u]
#Homo projector overlaps
# Homo projector overlaps
mol = range(len(slab))[-2:]
p_uai = [dict([(mol[a], P_ni[4]) for a, P_ni in kpt.P_ani.items()])
for kpt in c_mol.wfs.kpt_u]
# Slab with adsorbed molecule
#-----------------------------------
# Slab with adsorbed molecule
slab.set_calculator(calc)
orbital = dscf.AEOrbital(calc, wf_u, p_uai, Estart=-100.0, Eend=0.0)
dscf.dscf_calculation(calc, [[-1.0, orbital, 1]], slab)
slab.get_potential_energy()
#!/usr/bin/env python
from numpy import reshape, dot
from ase.visualize import view
......@@ -8,30 +6,33 @@ from gpaw import GPAW
from gpaw.mixer import MixerSum
import gpaw.dscf as dscf
filename='lumo'
#-------------------------------------------
filename = 'lumo'
c_mol = GPAW(nbands=9, h=0.2, xc='RPBE', kpts=(8,6,1),
c_mol = GPAW(nbands=9,
h=0.2,
xc='RPBE',
kpts=(8, 6, 1),
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -2}, txt='CO_lumo.txt')
'bands': -2},
txt='CO_lumo.txt')
calc = GPAW(nbands=60, h=0.2, xc='RPBE', kpts=(8,6,1),
calc = GPAW(nbands=80,
h=0.2,
xc='RPBE',
kpts=(8, 6, 1),
eigensolver='cg',
spinpol=True,
mixer=MixerSum(nmaxold=5, beta=0.1, weight=100),
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-7,
'bands': -10}, txt=filename+'.txt')
'bands': -10},
txt=filename + '.txt')
#----------------------------------------
# Import Slab with relaxed CO
#slab = Calculator('gs.gpw').get_atoms()
# Import Slab with relaxed CO
slab = fcc111('Pt', size=(1, 2, 3), orthogonal=True)
add_adsorbate(slab, 'C', 2.0, 'ontop')
add_adsorbate(slab, 'O', 3.15, 'ontop')
......@@ -41,25 +42,24 @@ view(slab)
molecule = slab.copy()
del molecule [:-2]
del molecule[:-2]
# Molecule
#----------------
# Molecule
molecule.set_calculator(c_mol)
molecule.get_potential_energy()
#Find band corresponding to lumo
# Find band corresponding to lumo
lumo = c_mol.get_pseudo_wave_function(band=5, kpt=0, spin=1)
lumo = reshape(lumo, -1)
wf1_k = [c_mol.get_pseudo_wave_function(band=5, kpt=k, spin=1)
for k in range(len(c_mol.wfs.weight_k))]
for k in range(c_mol.wfs.kd.nibzkpts)]
wf2_k = [c_mol.get_pseudo_wave_function(band=6, kpt=k, spin=1)
for k in range(len(c_mol.wfs.weight_k))]
for k in range(c_mol.wfs.kd.nibzkpts)]
band_k = []
for k in range(len(c_mol.wfs.weight_k)):
for k in range(c_mol.wfs.kd.nibzkpts):
wf1 = reshape(wf1_k[k], -1)
wf2 = reshape(wf2_k[k], -1)
p1 = abs(dot(wf1, lumo))
......@@ -69,16 +69,15 @@ for k in range(len(c_mol.wfs.weight_k)):
else:
band_k.append(6)
#Lumo wavefunction
# Lumo wavefunction
wf_u = [kpt.psit_nG[band_k[kpt.k]] for kpt in c_mol.wfs.kpt_u]
#Lumo projector overlaps
# Lumo projector overlaps
mol = range(len(slab))[-2:]
p_uai = [dict([(mol[a], P_ni[band_k[kpt.k]]) for a, P_ni in kpt.P_ani.items()])
for kpt in c_mol.wfs.kpt_u]
# Slab with adsorbed molecule
#-----------------------------------
slab.set_calculator(calc)
orbital = dscf.AEOrbital(calc, wf_u, p_uai)
dscf.dscf_calculation(calc, [[1.0, orbital, 1]], slab)
......
......@@ -59,6 +59,7 @@ Requirements
* Python_ 2.7, 3.4-
* NumPy_ 1.6.1 or later (base N-dimensional array package)
* SciPy_ 0.7 or later (library for scientific computing)
* ASE_ 3.15.0 or later (atomic simulation environment)
* a C-compiler
* LibXC_ 2.0.1 or later
......@@ -66,8 +67,6 @@ Requirements
Optional, but highly recommended:
* SciPy_ 0.7 or later (library for scientific computing, requirered for
some features)
* an MPI_ library (required for parallel calculations)
* FFTW_ (for increased performance)
* BLACS_ and ScaLAPACK_
......@@ -175,7 +174,7 @@ Sou can get the source from a tar-file or from Git:
$ tar -xf gpaw-1.3.0.tar.gz
$ ln -s gpaw-1.3.0 gpaw
Here is a `list of tarballs <https://pypi.python.org/simple/gpaw/>`__.
Here is a `list of tarballs <https://pypi.org/simple/gpaw/>`__.
:Git clone:
......
......@@ -10,6 +10,11 @@ Git master branch
:git:`master <>`.
* Corresponding ASE release: ASE-3.16.0
* Improved parallelization of operations with localized functions in
PW mode. This solves the current size bottleneck in PW mode.
* Added QNA XC functional.
* Major refactoring of the LCAOTDDFT code and added Kohn--Sham decomposition
......
......@@ -20,7 +20,7 @@ Next we calculate the dynamical dielectric function using the Bethe-Salpeter equ
.. image:: bse_Si.png
:height: 400 px
The Hamiltonian is stored in ``H_SS.gpw`` by default and its eigenstates and eigenvalues are stored in ``v_TS.gpw``. These are easily read by the get_dielectric_function method, such that if you want to calculate the spectrum at a different energy range or broadening you can simply do::
The ``write_v`` keyword ensures that the eigenstates and eigenvalues are stored in ``v_TS.gpw``. These are easily read by the get_dielectric_function method, such that if you want to calculate the spectrum at a different energy range or broadening you can simply do::
bse.get_dielectric_function(filename='bse_0.1.csv',
readfile='v_TS.gpw',
......@@ -33,20 +33,20 @@ The parameters that needs to be converged in the calculation are the k-points in
For large calculations, it may be useful to write the screened interaction, which is the first quantity that is calculated and may be restarted in a subsequent calculation. This may be done with the keyword ``wfile='W_qGG.pckl'``, where the .pckl file contains the screened interaction matrices at all q-points and a few other variables. It may also be useful to set ``write_h=False`` and ``write_v=False``, since these files may become quite large for big calculations.
Excitons in monolayer MoS2
=======================================
Excitons in monolayer MoS2 with Spin-orbit Coupling
===================================================
Spectrum from the Bethe-Salpeter equation
-----------------------------------------
The screening plays a fundamental role in the Bethe-Salpeter equation and for 2D systems the screening requires a special treatment. In particular we use a truncated Coulomb interaction inorder to decouple the screening between periodic images. We refer to Ref. [#Huser]_ for details on the truncated Coulomb interaction in GPAW. As before, we calculate the ground state of `MoS_2` with the script :download:`gs_MoS2.py`, which takes on the order of one hour on 4 CPUs. Note the large density of k-points, which are required to converge the BSE spectrum of two-dimensional systems.
The screening plays a fundamental role in the Bethe-Salpeter equation and for 2D systems the screening requires a special treatment. In particular we use a truncated Coulomb interaction inorder to decouple the screening between periodic images. We refer to Ref. [#Huser]_ for details on the truncated Coulomb interaction in GPAW. As before, we calculate the ground state of `MoS_2` with the script :download:`gs_MoS2.py`, which takes a few minutes. Note the large density of k-points, which are required to converge the BSE spectrum of two-dimensional systems.
The macroscopic dielectric function is calculated as an average of the microscopic screening over the unit cell. Clearly, for a 2D system this will depend on the unit cell size in the direction orthogonal to the slab and in the converged limit the dielectric function becomes unity. Instead we may calculate the longitudinal part of 2D polarizability which is independent of unit cell size. This is done in RPA as well as BSE with the scripts :download:`pol_MoS2.py`, which takes ~10 hours on 16 CPUs. Note that the BSE polarizability is calculated with and without Coulomb truncation for comparison. The results can be plottet with :download:`plot_MoS2.py` and is shown below.
The macroscopic dielectric function is calculated as an average of the microscopic screening over the unit cell. Clearly, for a 2D system this will depend on the unit cell size in the direction orthogonal to the slab and in the converged limit the dielectric function becomes unity. Instead we may calculate the longitudinal part of 2D polarizability which is independent of unit cell size. This is done in RPA as well as BSE with the scripts :download:`pol_MoS2.py`, which takes ~20 hours on 16 CPUs. Note that the BSE polarizability is calculated with and without Coulomb truncation for comparison. In both case spin-orbit coupling is included through the ``spinors`` keyword. We refer to Ref. [#Olsenspin]_ for details on the spin-orbit implementation. The results can be plottet with :download:`plot_MoS2.py` and is shown below.
.. image:: bse_MoS2.png
:height: 400 px
The excitonic effects are much stronger than in the case of Si due to the reduced screening in 2D. In particular, we can identify a distinct exciton well below the band edge. Note that without Coulomb truncation, the BSE spectrum is shifted upward in energy due the screening of electron-hole interactions from periodic images.
The excitonic effects are much stronger than in the case of Si due to the reduced screening in 2D. In particular, we can identify a distinct spin-orbit split exciton well below the band edge. Note that without Coulomb truncation, the BSE spectrum is shifted upward in energy due the screening of electron-hole interactions from periodic images.
2D screening with and without Coulomb truncation
------------------------------------------------
......@@ -84,5 +84,9 @@ which is in good agreement with the BSE computation above
.. [#Huser] F. Huser, T. Olsen and K. S. Thygesen
*Phys. Rev. B* **88**, 245309 (2013)
.. [#Olsenspin] T. Olsen
*Phys. Rev. B.* **94**, 235106 (2016)
.. [#Olsen] T. Olsen, S. Latini, F. Rasmussen and K. S. Thygesen
*Phys. Rev. Lett.* **116**, 056401 (2016)
......@@ -25,6 +25,7 @@ bse = BSE('gs_Si.gpw',
nbands=50,
eshift=eshift,
mode='BSE',
write_v=True,
integrate_gamma=0,
txt='bse_Si.txt')
......
......@@ -3,16 +3,16 @@ from gpaw import GPAW, PW, FermiDirac
calc = GPAW(mode=PW(400),
xc='PBE',
nbands=70,
convergence={'bands': -20},
setups={'Mo': '6'},
parallel={'band': 1, 'domain': 1},
occupations=FermiDirac(width=0.001),
kpts={'size': (30, 30, 1), 'gamma': True},
kpts={'size': (42, 42, 1), 'gamma': True},
txt='gs_MoS2.txt')
slab = mx2(formula='MoS2', a=3.16, thickness=3.17, vacuum=5.0)
slab.set_calculator(calc)
slab.get_potential_energy()
calc.diagonalize_full_hamiltonian(nbands=50)
calc.write('gs_MoS2.gpw', mode='all')
......@@ -20,13 +20,13 @@ df.get_polarizability(filename='pol_rpa_MoS2.csv',
pbc=[True, True, False])
bse = BSE('gs_MoS2.gpw',
spinors=True,
ecut=ecut,
valence_bands=[8],
conduction_bands=[9],
nbands=50,
eshift=eshift,
mode='BSE',
integrate_gamma=1,
txt='bse_MoS2.txt')
bse.get_polarizability(filename='pol_bse_MoS2.csv',
......@@ -36,6 +36,7 @@ bse.get_polarizability(filename='pol_bse_MoS2.csv',
w_w=np.linspace(0, 5, 5001))
bse = BSE('gs_MoS2.gpw',
spinors=True,
ecut=ecut,
valence_bands=[8],
conduction_bands=[9],
......@@ -43,7 +44,6 @@ bse = BSE('gs_MoS2.gpw',
nbands=50,
eshift=eshift,
mode='BSE',
integrate_gamma=1,
txt='bse_MoS2_trun.txt')
bse.get_polarizability(filename='pol_bse_MoS2_trun.csv',
......
......@@ -8,7 +8,7 @@ a = Atoms('GaAs', cell=cell, pbc=True,
calc = GPAW(mode=PW(600),
xc='LDA',
occupations=FermiDirac(width=0.01, maxiter=10000000),
occupations=FermiDirac(width=0.01),
convergence={'density': 1.e-6},
symmetry='off',
kpts={'size': (2, 2, 2), 'gamma': True},
......
......@@ -44,6 +44,8 @@ class KohnShamConvergenceError(ConvergenceError):
class PoissonConvergenceError(ConvergenceError):
pass
class KPointError(Exception):
pass
def parse_extra_parameters(arg):
return {key.replace('-', '_'): value
......
......@@ -21,6 +21,7 @@ from gpaw.io.logger import GPAWLogger
from gpaw.io import Reader, Writer
from gpaw.jellium import create_background_charge
from gpaw.kpt_descriptor import KPointDescriptor
from gpaw.kpt_refine import create_kpoint_descriptor_with_refinement
from gpaw.kohnsham_layouts import get_KohnSham_layouts
from gpaw.matrix import suggest_blocking
from gpaw.occupations import create_occupation_number_object
......@@ -64,7 +65,9 @@ class GPAW(PAW, Calculator):
'eigensolver': None,
'background_charge': None,
'experimental': {'reuse_wfs_method': None,
'magmoms': None},
'magmoms': None,
'soc': None,
'kpt_refine': None},
'external': None,
'random': False,
'hund': False,
......@@ -400,12 +403,24 @@ class GPAW(PAW, Calculator):
self.wfs.set_eigensolver(None)
if key in ['mixer', 'verbose', 'txt', 'hund', 'random',
'eigensolver', 'idiotproof', 'experimental']:
'eigensolver', 'idiotproof']:
continue
if key in ['convergence', 'fixdensity', 'maxiter']:
continue
# Check nested arguments
if key in ['experimental']:
changed_parameters2 = changed_parameters[key]
for key2 in changed_parameters2:
if key2 in ['kpt_refine', 'magmoms', 'soc']:
self.wfs = None
elif key2 in ['reuse_wfs_method']:
continue
else:
raise TypeError('Unknown keyword argument:', key2)
continue
# More drastic changes:
if self.wfs:
self.wfs.set_orthonormalized(False)
......@@ -894,20 +909,42 @@ class GPAW(PAW, Calculator):
self.hamiltonian.soc = self.parameters.experimental.get('soc')
self.log(self.hamiltonian, '\n')
def create_wave_functions(self, mode, realspace,
nspins, collinear, nbands, nao, nvalence,
setups, cell_cv, pbc_c):
def create_kpoint_descriptor(self, nspins):
par = self.parameters
bzkpts_kc = kpts2ndarray(par.kpts, self.atoms)
kd = KPointDescriptor(bzkpts_kc, nspins)
kpt_refine = par.experimental.get('kpt_refine')
if kpt_refine is None:
kd = KPointDescriptor(bzkpts_kc, nspins)
self.timer.start('Set symmetry')
kd.set_symmetry(self.atoms, self.symmetry, comm=self.world)
self.timer.stop('Set symmetry')
self.timer.start('Set symmetry')
kd.set_symmetry(self.atoms, self.symmetry, comm=self.world)
self.timer.stop('Set symmetry')
else:
self.timer.start('Set k-point refinement')
kd = create_kpoint_descriptor_with_refinement(
kpt_refine,
bzkpts_kc, nspins, self.atoms,
self.symmetry, comm=self.world,
timer=self.timer)
self.timer.stop('Set k-point refinement')
# Update quantities which might have changed, if symmetry
# was changed
self.symmetry = kd.symmetry
self.setups.set_symmetry(kd.symmetry)
self.log(kd)
return kd
def create_wave_functions(self, mode, realspace,
nspins, collinear, nbands, nao, nvalence,
setups, cell_cv, pbc_c):
par = self.parameters
kd = self.create_kpoint_descriptor(nspins)
parallelization = mpi.Parallelization(self.world,
nspins * kd.nibzkpts)
......
......@@ -197,7 +197,9 @@ class Density:
self.ghat.set_positions(spos_ac)
self.mixer.reset()
self.nt_xg = None
self.nt_sg = None
self.nt_vg = None
self.nt_g = None
self.rhot_g = None
......
......@@ -114,7 +114,6 @@ class BlacsOrbitalLayouts(BlacsLayouts):
mcpus, ncpus, blocksize, timer)
nbands = bd.nbands
self.blocksize = blocksize
self.mynbands = mynbands = bd.mynbands
self.orbital_comm = self.bd.comm
self.naoblocksize = naoblocksize = -((-nao) // self.orbital_comm.size)
......@@ -131,7 +130,7 @@ class BlacsOrbitalLayouts(BlacsLayouts):
self.mMdescriptor = self.columngrid.new_descriptor(nao, nao,
naoblocksize, nao)
self.nMdescriptor = self.columngrid.new_descriptor(nbands, nao,