1
2import numpy as np
3
4from ase.build import molecule
5from ase.parallel import parprint
6
7from gpaw import GPAW
8from gpaw.analyse.multipole import Multipole
9from gpaw.cluster import Cluster
10from gpaw.test import equal
11
12
13def test_multipoleH2O(in_tmp_dir):
14    h = 0.3
15
16    s = Cluster(molecule('H2O'))
17    s.minimal_box(3., h)
18
19    gpwname = 'H2O_h' + str(h) + '.gpw'
20    try:
21        # XXX check why this fails in parallel
22        calc = GPAW(gpwname + 'failsinparallel', txt=None)
23        atoms = calc.get_atoms()
24        calc.density.set_positions(atoms.get_scaled_positions() % 1.0)
25        calc.density.interpolate_pseudo_density()
26        calc.density.calculate_pseudo_charge()
27    except IOError:
28        calc = GPAW(h=h, charge=0, txt=None)
29        calc.calculate(s)
30        calc.write(gpwname)
31
32    dipole_c = calc.get_dipole_moment()
33    parprint('Dipole', dipole_c)
34
35    center = np.array([1, 1, 1]) * 50.
36    mp = Multipole(center, calc, lmax=2)
37    q_L = mp.expand(-calc.density.rhot_g)
38    parprint('Multipole', q_L)
39
40    # The dipole moment is independent of the center
41    equal(dipole_c[2], q_L[2], 1e-10)
42
43    mp.to_file(calc, mode='w')
44