Coverage for gpaw/test/exx/test_forces_ut.py: 100%
119 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import pytest
2import numpy as np
3from ase import Atoms
4from ase.units import Ha, Bohr
6from gpaw.kpt_descriptor import KPointDescriptor
7from gpaw.grid_descriptor import GridDescriptor
8from gpaw.kpoint import KPoint
9from gpaw.symmetry import Symmetry
10from gpaw.wavefunctions.arrays import PlaneWaveExpansionWaveFunctions
11from gpaw.pw.descriptor import PWDescriptor
12from gpaw.pw.lfc import PWLFC
13from gpaw.projections import Projections
14from gpaw.mpi import world
15from gpaw.spline import Spline
16from gpaw.hybrids.energy import calculate_energy
17from gpaw.hybrids.forces import calculate_forces
18from gpaw.hybrids.coulomb import coulomb_interaction
19from gpaw.hybrids.paw import calculate_paw_stuff
20from gpaw.hybrids.symmetry import Symmetry as Sym
21from gpaw.hybrids.kpts import get_kpt
24N = 20
25L = 2.5
26nb = 2
27r2 = np.linspace(0, 2, 101)**2
30@pytest.mark.ci
31@pytest.mark.skipif(world.size > 1, reason='Not parallelized')
32def test_force():
33 x = 0.2
34 y = 0.35
35 z = 0.1
36 calc = Calc([[x, y, z]])
37 omega = 0.2
38 wfs = calc.wfs
39 kd = wfs.kd
40 coulomb = coulomb_interaction(omega, wfs.gd, kd)
41 sym = Sym(kd)
42 paw_s = calculate_paw_stuff(wfs, calc.density)
43 F_v = calculate_forces(wfs, coulomb, sym, paw_s)[0] * Ha / Bohr
44 print(F_v)
45 dx = 0.0001
46 y -= dx / 2
47 calc = Calc([[x, y, z]])
48 em = energy(calc, sym, coulomb)
49 y += dx
50 calc = Calc([[x, y, z]])
51 ep = energy(calc, sym, coulomb)
52 error = (em - ep) / (dx * L * Bohr) - F_v[1]
53 print(error)
54 assert abs(error) < 6e-10
57class Calc:
58 def __init__(self, spos_ac):
59 self.spos_ac = np.array(spos_ac)
60 self.wfs = WFS(self.spos_ac)
61 self.density = Dens(self.wfs)
62 self.hamiltonian = Ham()
65class WFS:
66 nspins = 1
67 world = world
69 def __init__(self, spos_ac):
70 self.spos_ac = spos_ac
71 self.gd = GridDescriptor([N, N, N], np.eye(3) * L)
72 sym = Symmetry([], self.gd.cell_cv)
73 self.kd = KPointDescriptor(None)
74 self.kd.set_symmetry(Atoms(pbc=True), sym)
75 self.setups = [Setup()]
76 pd = PWDescriptor(50, self.gd, complex, self.kd)
77 data = pd.zeros(nb)
78 R = self.gd.get_grid_point_coordinates()
79 R -= L / 2
80 R2 = (R**2).sum(0)
81 data[0] = pd.fft(np.exp(-R2 * (10 / L**2)))
82 data[1] = pd.fft(np.exp(-R2 * (10 / L**2))) + 0.5
83 psit = PlaneWaveExpansionWaveFunctions(nb, pd, data=data)
84 proj = Projections(nb, [1], AP(), world, dtype=complex)
85 self.pt = PWLFC([[Spline.from_data(0, 2.0, np.exp(-r2 * 10))]], pd)
86 self.pt.set_positions(self.spos_ac)
87 psit.matrix_elements(self.pt, out=proj)
88 f_n = np.array([1.0, 0.5])
89 kpt = KPoint(1.0, 1.0, 0, 0, 0, None)
90 kpt.psit = psit
91 kpt.projections = proj
92 kpt.f_n = f_n
93 self.kpt_u = [kpt]
94 self.kpt_qs = [[kpt]]
97class Ham:
98 e_total_free = -100.0
99 e_xc = -20.0
102class Dens:
103 nspins = 1
105 def __init__(self, wfs):
106 kpt = wfs.kpt_u[0]
107 self.D_asp = D(kpt)
108 self.finegd = wfs.gd.refine()
109 self.nt_sg = self.finegd.zeros(1)
112class D(dict):
113 def __init__(self, kpt):
114 dict.__init__(self)
115 P_ni = kpt.P_ani[0]
116 self[0] = kpt.f_n.dot(P_ni * P_ni.conj()).real[np.newaxis]
117 self.partition = AP()
120class AP:
121 my_indices = [0]
122 comm = world
123 rank_a = [0]
126class Setup:
127 Delta_iiL = np.zeros((1, 1, 1)) + 0.1
128 X_p = np.zeros(1) + 0.3
129 ExxC = -10.0
130 ghat_l = [Spline.from_data(0, 2.0, np.exp(-r2 * 10))]
131 xc_correction = None
132 M_pp = np.zeros((1, 1)) + 0.3
135def energy(calc, sym, coulomb):
136 wfs = calc.wfs
137 paw_s = calculate_paw_stuff(wfs, calc.density)
138 kpts = [get_kpt(wfs, 0, 0, 0, nb)]
139 e1, e2 = calculate_energy(kpts, paw_s[0], wfs, sym, coulomb, wfs.spos_ac)
140 evc = e1 * 2 * Ha
141 evv = e2 * 2 * Ha
142 return evc + evv