Commit 871c8d27 authored by Ask Hjorth Larsen's avatar Ask Hjorth Larsen

Merge branch 'efield' into 'master'

efield for pw

See merge request !296
parents 9491e710 0c641c34
Pipeline #21313079 passed with stage
in 2 minutes and 35 seconds
...@@ -19,6 +19,7 @@ def create_external_potential(name, **kwargs): ...@@ -19,6 +19,7 @@ def create_external_potential(name, **kwargs):
class ExternalPotential: class ExternalPotential:
vext_g = None vext_g = None
vext_q = None
def get_potential(self, gd): def get_potential(self, gd):
"""Get the potential on a regular 3-d grid. """Get the potential on a regular 3-d grid.
...@@ -29,6 +30,17 @@ class ExternalPotential: ...@@ -29,6 +30,17 @@ class ExternalPotential:
self.calculate_potential(gd) self.calculate_potential(gd)
return self.vext_g return self.vext_g
def get_potentialq(self, gd, pd3):
"""Get the potential on a regular 3-d grid in real space.
Will only call calculate_potential() the first time."""
if self.vext_q is None:
vext_g = self.get_potential(gd)
self.vext_q = pd3.fft(vext_g)
return self.vext_q
def calculate_potential(self, gd): def calculate_potential(self, gd):
raise NotImplementedError raise NotImplementedError
...@@ -69,11 +81,12 @@ class ConstantElectricField(ExternalPotential): ...@@ -69,11 +81,12 @@ class ConstantElectricField(ExternalPotential):
def calculate_potential(self, gd): def calculate_potential(self, gd):
d_v = self.field_v / (self.field_v**2).sum()**0.5 d_v = self.field_v / (self.field_v**2).sum()**0.5
for axis_v in gd.cell_cv[gd.pbc_c]: # Currently skipped, PW mode is periodic in all directions
if abs(np.dot(d_v, axis_v)) > self.tolerance: #for axis_v in gd.cell_cv[gd.pbc_c]:
raise ValueError( # if abs(np.dot(d_v, axis_v)) > self.tolerance:
'Field not perpendicular to periodic axis: {}' # raise ValueError(
.format(axis_v)) # 'Field not perpendicular to periodic axis: {}'
# .format(axis_v))
center_v = 0.5 * gd.cell_cv.sum(0) center_v = 0.5 * gd.cell_cv.sum(0)
r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0)) r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
......
...@@ -194,6 +194,7 @@ tests = [ ...@@ -194,6 +194,7 @@ tests = [
'broydenmixer.py', # ~3s 'broydenmixer.py', # ~3s
'pw/fulldiagk.py', # ~3s 'pw/fulldiagk.py', # ~3s
'ext_potential/external.py', # ~3s 'ext_potential/external.py', # ~3s
'ext_potential/external_pw.py', # ~3s
'lcao/atomic_corrections.py', # ~3s 'lcao/atomic_corrections.py', # ~3s
'vdw/libvdwxc_h2.py', # ~3s 'vdw/libvdwxc_h2.py', # ~3s
'generic/mixer.py', # ~3s 'generic/mixer.py', # ~3s
......
from __future__ import print_function
from ase import Atoms
from gpaw import GPAW
from gpaw.test import equal
from gpaw.wavefunctions.pw import PW
from gpaw.external import ConstantPotential
cp = ConstantPotential()
sc = 2.9
R = 0.7 # approx. experimental bond length
R = 1.0
a = 2 * sc
c = 3 * sc
at = 'H'
H2 = Atoms('HH', [(a / 2, a / 2, (c - R) / 2),
(a / 2, a / 2, (c + R) / 2)],
cell=(a, a, c),
pbc=False)
txt = None
convergence = {'eigenstates': 1.e-6 * 40 * 1.5**3,
'density': 1.e-2,
'energy': 0.1}
# without potential
if True:
if txt:
print('\n################## no potential')
c00 = GPAW(mode=PW(200), nbands=-1,
convergence=convergence,
txt=txt)
c00.calculate(H2)
eps00_n = c00.get_eigenvalues()
# 0 potential
if True:
if txt:
print('\n################## 0 potential')
cp0 = ConstantPotential(0.0)
c01 = GPAW(mode=PW(200), nbands=-2, external=cp0,
convergence=convergence,
txt=txt)
c01.calculate(H2)
# 1 potential
if True:
if txt:
print('################## 1 potential')
cp1 = ConstantPotential(-1.0)
c1 = GPAW(mode=PW(200), nbands=-2, external=cp1,
convergence=convergence,
txt=txt)
c1.calculate(H2)
for i in range(c00.get_number_of_bands()):
f00 = c00.get_occupation_numbers()[i]
if f00 > 0.01:
e00 = c00.get_eigenvalues()[i]
e1 = c1.get_eigenvalues()[i]
print('Eigenvalues no pot, expected, error=',
e00, e1 + 1, e00 - e1 - 1)
equal(e00, e1 + 1., 0.007)
E_c00 = c00.get_potential_energy()
E_c1 = c1.get_potential_energy()
DeltaE = E_c00 - E_c1
equal(DeltaE, 0, 0.002)
...@@ -1581,6 +1581,14 @@ class ReciprocalSpaceHamiltonian(Hamiltonian): ...@@ -1581,6 +1581,14 @@ class ReciprocalSpaceHamiltonian(Hamiltonian):
self.epot = 0.5 * self.pd3.integrate(self.vHt_q, dens.rhot_q) self.epot = 0.5 * self.pd3.integrate(self.vHt_q, dens.rhot_q)
self.vt_Q = self.vbar_Q + self.vHt_q[dens.G3_G] / 8 self.vt_Q = self.vbar_Q + self.vHt_q[dens.G3_G] / 8
self.e_external = 0.0
if self.vext is not None:
gd = self.finegd
vext_q = self.vext.get_potentialq(gd, self.pd3)
self.vt_Q += vext_q[dens.G3_G] / 8
self.e_external = self.pd3.integrate(vext_q, dens.rhot_q)
self.vt_sG[:] = self.pd2.ifft(self.vt_Q) self.vt_sG[:] = self.pd2.ifft(self.vt_Q)
self.timer.start('XC 3D grid') self.timer.start('XC 3D grid')
...@@ -1602,15 +1610,17 @@ class ReciprocalSpaceHamiltonian(Hamiltonian): ...@@ -1602,15 +1610,17 @@ class ReciprocalSpaceHamiltonian(Hamiltonian):
self.timer.stop('XC 3D grid') self.timer.stop('XC 3D grid')
eext = 0.0 return np.array([self.epot, self.ebar, self.e_external, self.exc])
return np.array([self.epot, self.ebar, eext, self.exc])
def calculate_atomic_hamiltonians(self, density): def calculate_atomic_hamiltonians(self, density):
W_aL = {} W_aL = {}
for a in density.D_asp: for a in density.D_asp:
W_aL[a] = np.empty((self.setups[a].lmax + 1)**2) W_aL[a] = np.empty((self.setups[a].lmax + 1)**2)
density.ghat.integrate(self.vHt_q, W_aL) if self.vext:
vext_q = self.vext.get_potentialq(self.finegd, self.pd3)
density.ghat.integrate(self.vHt_q+vext_q, W_aL)
else:
density.ghat.integrate(self.vHt_q, W_aL)
return W_aL return W_aL
def calculate_kinetic_energy(self, density): def calculate_kinetic_energy(self, density):
...@@ -1636,7 +1646,11 @@ class ReciprocalSpaceHamiltonian(Hamiltonian): ...@@ -1636,7 +1646,11 @@ class ReciprocalSpaceHamiltonian(Hamiltonian):
restrict_and_collect = restrict restrict_and_collect = restrict
def calculate_forces2(self, dens, ghat_aLv, nct_av, vbar_av): def calculate_forces2(self, dens, ghat_aLv, nct_av, vbar_av):
dens.ghat.derivative(self.vHt_q, ghat_aLv) if self.vext:
vext_q = self.vext.get_potentialq(self.finegd, self.pd3)
dens.ghat.derivative(self.vHt_q+vext_q, ghat_aLv)
else:
dens.ghat.derivative(self.vHt_q, ghat_aLv)
dens.nct.derivative(self.vt_Q, nct_av) dens.nct.derivative(self.vt_Q, nct_av)
self.vbar.derivative(dens.nt_Q, vbar_av) self.vbar.derivative(dens.nt_Q, vbar_av)
......
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