Commit 23bc70c0 authored by Morten Gjerding's avatar Morten Gjerding

Merge branch 'master' into fix-gpaw-python-command-extra-arguments

parents 702351fd 194fd824
Pipeline #100275374 passed with stage
in 3 minutes and 23 seconds
......@@ -14,5 +14,5 @@ master:
- su user -c 'gpaw test --range linalg/gemm_complex.py,lcao/dos.py'
- 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()) <= 2361, 'Please run flake8 on your code'"
- python -c "assert int(open('nerr').read()) <= 2364, 'Please run flake8 on your code'"
- python -We:invalid -m compileall -f -q gpaw/
......@@ -43,6 +43,7 @@ There are several ways to install GPAW:
* :ref:`siteconfig`
* Using :ref:`homebrew` on MacOSX
* Using :ref:`anaconda`
* This `docker image`_
* Tips and tricks for installation on many :ref:`platforms and
architectures`
* :ref:`troubleshooting`
......@@ -79,6 +80,7 @@ Optional, but highly recommended:
.. _PIP: https://pip.pypa.io/en/stable/
.. _ASE: https://wiki.fysik.dtu.dk/ase
.. _FFTW: http://www.fftw.org/
.. _docker image: https://hub.docker.com/r/marcindulak/gpaw-openmpi
.. _installation using pip:
......@@ -243,6 +245,7 @@ this to your ``siteconfig.py``::
libraries += ['fftw3']
.. _libxc installation:
Libxc Installation
......
......@@ -12,6 +12,9 @@ Git master branch
* Corresponding ASE release: ASE-3.18.1b1
* Self-consistent calculations with hybrid functionals are now possible in
plane-wave mode.
* We are now using setuptools_ instead of :mod:`distutils`.
This means that installation with pip works much better.
......
......@@ -22,8 +22,9 @@ class HydrogenAllElectronSetup(BaseSetup):
self.nao = None
self.pt_j = []
self.ni = 0
self.l_j = []
self.n_j = []
self.l_j = [0]
self.l_orb_j = [0]
self.n_j = [1]
self.nct = Spline(0, 0.5, [0.0, 0.0, 0.0])
self.Nct = 0.0
self.N0_p = []
......@@ -37,9 +38,10 @@ class HydrogenAllElectronSetup(BaseSetup):
v_g[0] = 4 * (alpha1**0.5 - alpha2**0.5)
self.vbar = Spline(0, rc, v_g)
self.Delta_pL = np.zeros((0, 1))
self.Delta_iiL = np.zeros((0, 0, 1))
self.Delta0 = -1 / (4 * np.pi)**0.5
self.lmax = 0
self.K_p = self.M_p = self.MB_p = np.zeros(0)
self.K_p = self.M_p = self.MB_p = self.X_p = np.zeros(0)
self.M_pp = np.zeros((0, 0))
self.Kc = 0.0
self.MB = 0.0
......@@ -50,6 +52,7 @@ class HydrogenAllElectronSetup(BaseSetup):
self.type = 'all-electron'
self.fingerprint = None
self.symbol = 'H'
self.ExxC = 0.0
def get_default_nbands(self):
return 1
......
......@@ -42,9 +42,11 @@ class AppelbaumHamann(BaseSetup):
self.ghat_l = [Spline(0, rc, 4 * alpha**1.5 / np.pi**0.5 * x_g)]
self.vbar = Spline(0, rc, 2 * np.pi**0.5 * (v1 + v2 * r2_g) * x_g)
self.Delta_pL = np.zeros((1, 1))
self.Delta_iiL = np.zeros((1, 1, 1))
self.Delta0 = -4 / (4 * np.pi)**0.5
self.ExxC = 0.0
self.lmax = 0
self.K_p = self.M_p = self.MB_p = np.zeros(1)
self.K_p = self.M_p = self.MB_p = self.X_p = np.zeros(1)
self.M_pp = np.zeros((1, 1))
self.Kc = 0.0
self.MB = 0.0
......
......@@ -647,12 +647,12 @@ class AllElectronAtom:
states = []
for ch in self.channels:
for n, f in enumerate(ch.f_n):
states.append((ch.e_n[n], ch, n))
states.append((ch.e_n[n], n, ch.s, ch))
states.sort()
for e, ch, n in states:
for e, n, s, ch in states:
name = str(n + ch.l + 1) + ch.name
if self.nspins == 2:
name += '(%s)' % '+-'[ch.s]
name += '(%s)' % '+-'[s]
n_g = ch.calculate_density(n)
rave = self.rgd.integrate(n_g, 1)
self.log(' %-7s %6.3f %13.6f %13.5f %6.3f' %
......
from .hybrid import HybridXC
__all__ = ['HybridXC']
import numpy as np
from math import pi
class ShortRangeCoulomb:
def __init__(self, omega):
self.omega = omega
def get_potential(self, pd):
G2_G = pd.G2_qG[0]
x_G = 1 - np.exp(-G2_G / (4 * self.omega**2))
with np.errstate(invalid='ignore'):
v_G = 4 * pi * x_G / G2_G
if pd.kd.gamma:
v_G[0] = pi / self.omega**2
return v_G
This diff is collapsed.
This diff is collapsed.
......@@ -19,7 +19,7 @@ def matrix_matrix_multiply(alpha, a, opa, b, opb, beta=0.0, c=None,
equivalent PBLAS functions for distributed matrices.
The coefficients alpha and beta are of type float. Matrices a, b and c
must have same type (float or complex). The strings apa and opb must be
must have same type (float or complex). The strings opa and opb must be
'N', 'T', or 'C' . For opa='N' and opb='N', the operation performed is
equivalent to::
......
......@@ -2,11 +2,12 @@ import numpy as np
from gpaw.matrix import Matrix
from gpaw.mpi import serial_comm
from gpaw.utilities.partition import AtomPartition
class Projections:
def __init__(self, nbands, nproj_a, atom_partition, bcomm,
collinear=True, spin=0, dtype=float):
collinear=True, spin=0, dtype=float, data=None):
self.nproj_a = np.asarray(nproj_a)
self.atom_partition = atom_partition
self.bcomm = bcomm
......@@ -28,7 +29,8 @@ class Projections:
if not collinear:
I1 *= 2
self.matrix = Matrix(nbands, I1, dtype, dist=(bcomm, bcomm.size, 1))
self.matrix = Matrix(nbands, I1, dtype, data,
dist=(bcomm, bcomm.size, 1))
if collinear:
self.myshape = self.matrix.array.shape
......@@ -53,6 +55,12 @@ class Projections:
self.atom_partition if atom_partition is None else atom_partition,
bcomm, self.collinear, self.spin, self.matrix.dtype)
def view(self, n1: int, n2: int) -> 'Projections':
return Projections(n2 - n1, self.nproj_a,
self.atom_partition,
self.bcomm, self.collinear, self.spin,
self.matrix.dtype, self.matrix.array[n1:n2])
def items(self):
for a, I1, I2 in self.indices:
yield a, self.array[..., I1:I2]
......@@ -64,8 +72,19 @@ class Projections:
def __contains__(self, a):
return a in self.map
def todicttttt(self):
return dict(self.items())
def broadcast(self) -> 'Projections':
ap = AtomPartition(serial_comm, np.zeros(len(self.nproj_a), int))
P = self.new(atom_partition=ap)
comm = self.atom_partition.comm
for a, rank in enumerate(self.atom_partition.rank_a):
P1_ni = P[a]
if comm.rank == rank:
P_ni = self[a].copy()
else:
P_ni = np.empty_like(P1_ni)
comm.broadcast(P_ni, rank)
P1_ni[:] = P_ni
return P
def redist(self, atom_partition):
P = self.new(atom_partition=atom_partition)
......
from __future__ import division, print_function
import functools
import itertools
import os
......
......@@ -571,11 +571,13 @@ class PWSymmetryAnalyzer:
reds = s % self.nU
if self.timereversal(s):
TR = lambda x: x.conj()
TR = np.conj
sign = -1
else:
sign = 1
TR = lambda x: x
def TR(x):
return x
return U_scc[reds], sign, TR, self.shift_sc[s], ft_sc[reds]
......@@ -1209,7 +1211,6 @@ class PairDensity:
n_R = ut1cc_R * ut_R
with self.timer('fft'):
n_G[:] = pd.fft(n_R, 0, Q_G) * dv
# PAW corrections:
with self.timer('gemm'):
for C1_Gi, P2_mi in zip(C1_aGi, kpt2.P_ani):
......@@ -1289,8 +1290,8 @@ class PairDensity:
vel_nv = np.zeros((nb - na, 3), dtype=complex)
if n_n is None:
n_n = np.arange(na, nb)
assert np.max(n_n) < nb, print('This is too many bands')
assert np.min(n_n) >= na, print('This is too few bands')
assert np.max(n_n) < nb, 'This is too many bands'
assert np.min(n_n) >= na, 'This is too few bands'
# Load kpoints
kd = self.calc.wfs.kd
......@@ -1431,16 +1432,21 @@ class PairDensity:
shift_c = shift_c.round().astype(int)
if (U_cc == np.eye(3)).all():
T = lambda f_R: f_R
def T(f_R):
return f_R
else:
N_c = self.calc.wfs.gd.N_c
i_cr = np.dot(U_cc.T, np.indices(N_c).reshape((3, -1)))
i = np.ravel_multi_index(i_cr, N_c, 'wrap')
T = lambda f_R: f_R.ravel()[i].reshape(N_c)
def T(f_R):
return f_R.ravel()[i].reshape(N_c)
if time_reversal:
T0 = T
T = lambda f_R: T0(f_R).conj()
def T(f_R):
return T0(f_R).conj()
shift_c *= -1
a_a = []
......
......@@ -9,7 +9,6 @@ See:
Wigner-Seitz truncated interactions: Towards chemical accuracy
in nontrivial systems
"""
from __future__ import print_function
import sys
from math import pi
......
......@@ -110,7 +110,9 @@ class SCFLoop:
errors['density'] = 0.0
if len(self.old_energies) >= 3:
errors['energy'] = np.ptp(self.old_energies[-3:])
energies = self.old_energies[-3:]
if np.isfinite(energies).all():
errors['energy'] = np.ptp(energies)
# We only want to calculate the (expensive) forces if we have to:
check_forces = (self.max_errors['force'] < np.inf and
......@@ -190,10 +192,13 @@ class SCFLoop:
else:
log(' ', end='')
log('%11.6f %-5s %-7s' %
(Ha * ham.e_total_extrapolated,
niterocc,
niterpoisson), end='')
if np.isfinite(ham.e_total_extrapolated):
energy = '{:11.6f}'.format(Ha * ham.e_total_extrapolated)
else:
energy = ' ' * 11
log('%s %-5s %-7s' %
(energy, niterocc, niterpoisson), end='')
if wfs.nspins == 2:
log(' %+.4f' % occ.magmom, end='')
......
......@@ -294,7 +294,7 @@ class BaseSetup:
D_sii = np.zeros((nspins, ni, ni))
D_sp = np.zeros((nspins, ni * (ni + 1) // 2))
nj = len(self.l_j)
nj = len(self.pt_j)
j = 0
i = 0
ib = 0
......
......@@ -337,6 +337,8 @@ tests = [
'tpss.py', # ~18s
'tddft/td_na2.py', # ~18s
'exx/coarse.py', # ~18s
'exx/double_cell.py',
'exx/derivs.py',
'corehole/si.py', # ~18s
'mgga/mgga_sc.py', # ~19s
'Hubbard_U_Zn.py', # ~20s
......
from myqueue.task import task
def create_tasks():
return [task('kpts.py')]
import numpy as np
from ase import Atoms
from ase.units import Ha
from gpaw import GPAW, PW
from gpaw.hybrids import HybridXC
from gpaw.xc.exx import EXX
def test(kpts, setup, spinpol, symmetry):
a = Atoms('H2', cell=(3, 3, 3), pbc=1, positions=[[0, 0, 0], [0, 0, 0.75]])
a.calc = GPAW(mode=PW(100, force_complex_dtype=True),
setups=setup,
kpts=kpts,
spinpol=spinpol,
# nbands=1,
symmetry=symmetry,
txt=None,
xc='PBE')
a.get_potential_energy()
return a
def check(atoms, xc):
xc1 = HybridXC(xc)
c = atoms.calc
xc1.initialize(c.density, c.hamiltonian, c.wfs, c.occupations)
xc1.set_positions(c.spos_ac)
e = xc1.calculate_energy()
xc1.calculate_eigenvalues(0, 2, None)
xc2 = EXX(c, xc=xc, bands=(0, 2), txt=None)
xc2.calculate()
e0 = xc2.get_exx_energy()
eps0 = xc2.get_eigenvalue_contributions()
assert np.allclose(eps0, xc1.e_skn * Ha)
assert np.allclose(e0, e[0] + e[1])
ecv, evv, v_skn = xc1.test()
assert np.allclose(e0, ecv + evv)
assert np.allclose(v_skn, eps0)
def main():
for spinpol in [False, True]:
for setup in ['ae', 'paw']:
for symmetry in ['off', {}]:
for kpts in [(1, 1, 1),
(1, 1, 2),
(1, 1, 3),
(1, 1, 4),
(2, 2, 1),
[(0, 0, 0.5)],
[(0, 0, 0), (0, 0, 0.5)]]:
atoms = test(kpts, setup, spinpol, symmetry)
for xc in ['EXX', 'PBE0', 'HSE06']:
print(spinpol, setup, symmetry, kpts, xc,
len(atoms.calc.wfs.mykpts))
check(atoms, xc)
main()
from __future__ import print_function
import ase.db
from ase.units import kcal, mol
from ase.data.g2_1_ref import ex_atomization, atomization
......@@ -10,7 +9,7 @@ c = ase.db.connect('results.db')
atoms = {}
for d in c.select(natoms=1):
atoms[d.numbers[0]] = [d.energy, d.exx]
maepbe = 0.0
maeexx = 0.0
print(' PBE EXX')
......
......@@ -3,5 +3,5 @@ from myqueue.task import task
def create_tasks():
return [
task('molecules.py@1:1h'),
task('molecules.py', tmax='1h'),
task('check.py', deps='molecules.py')]
......@@ -67,9 +67,6 @@ for name in ['Si', 'C', 'GaAs', 'MgO', 'NaCl', 'Ar']:
kpts=kpts,
txt='%s.txt' % name)
if name in ['MgO', 'NaCl']:
atoms.calc.set(eigensolver='dav')
atoms.get_potential_energy()
atoms.calc.write(name + '.gpw', mode='all')
......
import numpy as np
from ase import Atoms
from gpaw import GPAW, PW, RMMDIIS
from gpaw.hybrids import HybridXC
from gpaw.mpi import world
def run(symbol, d0, M, ecut, L):
a = Atoms(2 * symbol, cell=[L, L, L])
a.positions[1, 0] = d0
a.center()
D = np.linspace(d0 * 0.98, d0 * 1.02, 5)
E = []
for d in D:
a.set_distance(0, 1, d)
a.calc = GPAW(
mode=PW(ecut, force_complex_dtype=True),
xc=HybridXC('PBE0'),
nbands=0,
# eigensolver='rmm-diis',
txt=f'{symbol}2-{d / d0:.2f}-{ecut}-{L}.txt')
e = a.get_potential_energy()
E.append(e)
p0 = np.polyfit(D, E, 3)
p1 = np.polyder(p0)
p2 = np.polyder(p1)
d = sorted([(np.polyval(p2, d), d) for d in np.roots(p1)])[1][1]
e2 = np.polyval(p0, d)
a = Atoms(symbol, cell=[L, L, L], magmoms=[M])
a.center()
a.calc = GPAW(
mode=PW(ecut, force_complex_dtype=True),
xc=HybridXC('PBE0', mix_all=False),
eigensolver=RMMDIIS(niter=1),
parallel={'band': 1, 'kpt': 1},
txt=f'{symbol}-{ecut}-{L}.txt')
e1 = a.get_potential_energy()
if world.rank == 0:
print(symbol, ecut, L, 2 * e1 - e2, d)
return 2 * e1 - e2, d
if __name__ == '__main__':
for L in [8, 10, 12]:
run('H', 0.75, 1, 500, L)
run('N', 1.089, 3, 500, L)
for ecut in [400, 500, 600, 700]:
run('H', 0.75, 1, ecut, 8)
run('N', 1.089, 3, ecut, 8)
def create_tasks():
from myqueue.task import task
return [task('gaps.py@16:5h'),
task('submit.agts.py', deps='gaps.py')]
task('submit.agts.py', deps='gaps.py'),
task('molecule.py', cores=8, tmax='5h')]
if __name__ == '__main__':
......
from __future__ import print_function
import sys
from ase import Atoms
......
import numpy as np
from ase import Atoms
from gpaw.kpt_descriptor import KPointDescriptor
from gpaw.grid_descriptor import GridDescriptor
from gpaw.response.wstc import WignerSeitzTruncatedCoulomb as WSTC
from gpaw.hybrids.exx import EXX, KPoint
from gpaw.symmetry import Symmetry
from gpaw.wavefunctions.arrays import PlaneWaveExpansionWaveFunctions
from gpaw.wavefunctions.pw import PWDescriptor, PWLFC
from gpaw.projections import Projections
from gpaw.mpi import world
from gpaw.spline import Spline
if world.size > 1:
from unittest import SkipTest
raise SkipTest
N = 20
L = 2.5
nb = 2
r2 = np.linspace(0, 1, 51)**2
spos_ac = np.zeros((1, 3)) + 0.25
class AP:
my_indices = [0]
comm = world
rank_a = [0]
class Setup:
Delta_iiL = np.zeros((1, 1, 1)) + 0.1
X_p = np.zeros(1) + 0.3
ExxC = -10.0
ghat_l = [Spline(0, 1.0, 1 - r2 * (1 - 2 * r2))]
gd = GridDescriptor([N, N, N], np.eye(3) * L)
sym = Symmetry([], gd.cell_cv)
kd = KPointDescriptor(None)
kd.set_symmetry(Atoms(pbc=True), sym)
coulomb = WSTC(gd.cell_cv, kd.N_c)
pd = PWDescriptor(10, gd, complex, kd)
data = pd.zeros(nb)
data[0, 1] = 3.0
data[1, 2] = -2.5
psit = PlaneWaveExpansionWaveFunctions(nb, pd, data=data)
proj = Projections(nb, [1], AP(), world, dtype=complex)
pt = PWLFC([[Spline(0, 1.0, 1 - r2 * (1 - 2 * r2))]], pd)
pt.set_positions(spos_ac)
f_n = np.array([1.0, 0.5])
kpt = KPoint(psit, proj, f_n, np.array([0.0, 0.0, 0.0]), 1.0)
xx = EXX(kd, [Setup()], pt, coulomb, spos_ac)
psit.matrix_elements(pt, out=proj)
C = 0.79
VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C
for a, P_ni in proj.items()}
x = xx.calculate([kpt], [kpt], VV_aii, derivatives=True)
v = x[3][0][0, 1]
eps = 0.00001
data[0, 1] = 3 + eps
psit.matrix_elements(pt, out=proj)
VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C
for a, P_ni in proj.items()}
xp = xx.calculate([kpt], [kpt], VV_aii)
data[0, 1] = 3 - eps
psit.matrix_elements(pt, out=proj)
VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C
for a, P_ni in proj.items()}
xm = xx.calculate([kpt], [kpt], VV_aii)
d = (xp[0] + xp[1] - xm[0] - xm[1]) / (2 * eps) * N**6 / L**3 / 2
assert abs(v - d) < 1e-10, (v, d)
from ase import Atoms
from gpaw import GPAW, PW, Davidson
L = 4.0
a = Atoms('H',
cell=[L, L, 1],
pbc=1)
a.center()
a.calc = GPAW(
mode=PW(400, force_complex_dtype=True),
parallel={'kpt': 1, 'band': 1},
eigensolver=Davidson(1),
kpts={'size': (1, 1, 2), 'gamma': True},
txt='H.txt',
xc='HSE06')
e1 = a.get_potential_energy()
eps1 = a.calc.get_eigenvalues(0)[0]
a *= (1, 1, 2)
es = Davidson(1)
es.keep_htpsit = True
a.calc = GPAW(
mode=PW(400, force_complex_dtype=True),
txt='H2.txt',
xc='HSE06')
e2 = a.get_potential_energy()
eps2 = a.calc.get_eigenvalues(0)[0]
assert abs(e2 - 2 * e1) < 0.002
assert abs(eps1 - eps2) < 0.001
from __future__ import print_function
from ase import Atom, Atoms
from gpaw import GPAW
from gpaw.test import equal
......@@ -7,11 +6,11 @@ a = 4.05
d = a / 2**0.5
bulk = Atoms([Atom('Al', (0, 0, 0)),
Atom('Al', (0, 0, d))],
cell=(4*d, 4*d, 2*d),
cell=(4 * d, 4 * d, 2 * d),
pbc=1)
n = 16
calc = GPAW(gpts=(2*n, 2*n, 1*n),
nbands=1*8,
calc = GPAW(gpts=(2 * n, 2 * n, 1 * n),
nbands=1 * 8,
kpts=(1, 1, 4),
convergence={'eigenstates': 2.3e-9},
xc='LDA')
......@@ -21,7 +20,9 @@ niter2 = calc.get_number_of_iterations()
bulk = bulk.repeat((1, 1, 2))
bulk.set_calculator(calc)
calc.set(nbands=16, kpts=(1, 1, 2), gpts=(2*n, 2*n, 2*n))
calc.set(nbands=16,
kpts=(1, 1, 2),
gpts=(2 * n, 2 * n, 2 * n))
e4 = bulk.get_potential_energy()
niter4 = calc.get_number_of_iterations()
......
from __future__ import print_function
from ase import Atom, Atoms
from ase import Atoms
from ase.calculators.test import numeric_force
from gpaw import GPAW, Mixer, FermiDirac, Davidson
from gpaw.test import equal
a = 4.0
n = 16
atoms = Atoms([Atom('H', [1.234, 2.345, 3.456])],
cell=(a, a, a), pbc=True)
atoms = Atoms('H',
positions=[[1.234, 2.345, 3.456]],
cell=(a, a, a),
pbc=True)
calc = GPAW(nbands=1,
gpts=(n, n, n),
txt=None,
......
......@@ -43,8 +43,8 @@ for l in range(4):
lfcr.set_positions(spos_ac)
lfcr.add(br, cr_axi)
a = pd.ifft(b)
ar = pdr.ifft(br)
a = pd.ifft(b[0])
ar = pdr.ifft(br[0])
equal(abs(a - ar).max(), 0, 1e-14)
if l == 0:
......
from __future__ import print_function
from ase import Atom, Atoms
from ase import Atoms
from gpaw import GPAW
from gpaw.test import equal
from gpaw.xc.tools import vxc
a = 5.0
d = 1.0
x = d / 3**0.5
atoms = Atoms([Atom('C', (0.0, 0.0, 0.0)),
Atom('H', (x, x, x)),
Atom('H', (-x, -x, x)),
Atom('H', (x, -x, -x)),
Atom('H', (-x, x, -x))],
cell=(a, a, a),
pbc=False)
atoms = Atoms('CH4',
[(0.0, 0.0, 0.0),
(x, x, x),
(-x, -x, x),
(x, -x, -x),
(-x, x, -x)],
cell=(a, a, a),
pbc=False)
atoms.positions[:] += a / 2
calc = GPAW(h=0.25, nbands=4, convergence={'eigenstates': 7.8e-10})
......@@ -30,7 +31,6 @@ niter_tolerance = 0
equal(energy, -23.631, energy_tolerance)