1from ase import Atoms
2from ase.units import Ha
3
4from gpaw import GPAW
5from gpaw.test import equal
6from gpaw.poisson import NoInteractionPoissonSolver
7from gpaw.external import ExternalPotential, known_potentials
8
9
10class HarmonicPotential(ExternalPotential):
11    def calculate_potential(self, gd):
12        a = gd.cell_cv[0, 0]
13        r_vg = gd.get_grid_point_coordinates()
14        self.vext_g = 0.5 * ((r_vg - a / 2)**2).sum(0)
15
16    def todict(self):
17        return {'name': 'HarmonicPotential'}
18
19
20def test_ext_potential_harmonic(in_tmp_dir):
21    """Test againts analytic result (no xc, no Coulomb)."""
22    a = 4.0
23    x = Atoms(cell=(a, a, a))  # no atoms
24
25    calc = GPAW(charge=-8,
26                nbands=4,
27                h=0.2,
28                xc={'name': 'null'},
29                external=HarmonicPotential(),
30                poissonsolver=NoInteractionPoissonSolver(),
31                eigensolver='cg')
32
33    x.calc = calc
34    x.get_potential_energy()
35
36    eigs = calc.get_eigenvalues()
37    equal(eigs[0], 1.5 * Ha, 0.002)
38    equal(abs(eigs[1:] - 2.5 * Ha).max(), 0, 0.003)
39
40    # Check write + read:
41    calc.write('harmonic.gpw')
42    known_potentials['HarmonicPotential'] = HarmonicPotential
43    GPAW('harmonic.gpw')
44