Coverage for gpaw/test/response/test_two-aluminum_chi_RPA.py: 99%
69 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import pytest
2import numpy as np
3import time
5from ase.build import bulk
6from ase.parallel import parprint
8from gpaw import GPAW, PW
9from gpaw.test import findpeak
10from gpaw.mpi import size, world
12from gpaw.response import ResponseGroundStateAdapter
13from gpaw.response.chiks import ChiKSCalculator
14from gpaw.response.susceptibility import ChiFactory
15from gpaw.response.pair_functions import read_susceptibility_array
18@pytest.mark.kspair
19@pytest.mark.response
20def test_response_two_aluminum_chi_RPA(in_tmp_dir):
21 assert size <= 4**3
23 # Ground state calculation
25 t1 = time.time()
27 a = 4.043
28 atoms1 = bulk('Al', 'fcc', a=a)
29 atoms2 = atoms1.repeat((2, 1, 1))
31 calc1 = GPAW(mode=PW(200),
32 nbands=12,
33 gpts=(12, 12, 12),
34 kpts=(8, 8, 8),
35 convergence={'density': 1e-5,
36 'eigenstates': 1e-7,
37 'bands': 8},
38 parallel={'domain': 1},
39 xc='LDA')
41 atoms1.calc = calc1
42 atoms1.get_potential_energy()
44 t2 = time.time()
46 calc2 = GPAW(mode=PW(200),
47 nbands=24,
48 gpts=(24, 12, 12),
49 kpts=(4, 8, 8),
50 convergence={'density': 1e-5,
51 'eigenstates': 1e-7,
52 'bands': 16},
53 parallel={'domain': 1},
54 xc='LDA')
56 atoms2.calc = calc2
57 atoms2.get_potential_energy()
59 t3 = time.time()
61 # Excited state calculation
62 q1_qc = [np.array([1 / 8., 0., 0.]), np.array([3 / 8., 0., 0.])]
63 q2_qc = [np.array([1 / 4., 0., 0.]), np.array([-1 / 4., 0., 0.])]
64 w = np.linspace(0, 24, 241)
66 # Calculate susceptibility using Al1
67 calculate_chi(calc1, q1_qc, w, 8, filename_prefix='Al1')
69 t4 = time.time()
71 # Calculate susceptibility using Al2
72 calculate_chi(calc2, q2_qc, w, 16, filename_prefix='Al2')
74 t5 = time.time()
76 parprint('')
77 parprint('Ground state calc 1 took', (t2 - t1), 'seconds')
78 parprint('Ground state calc 2 took', (t3 - t2), 'seconds')
79 parprint('Susceptibility calc 1 took', (t4 - t3), 'seconds')
80 parprint('Susceptibility calc 2 took', (t5 - t4), 'seconds')
82 # Check that results are consistent, when structure is simply repeated
84 # Read results
85 w11_w, G11_Gc, chi11_wGG = read_susceptibility_array('Al1_chiGG_q1.npz')
86 w21_w, G21_Gc, chi21_wGG = read_susceptibility_array('Al2_chiGG_q1.npz')
87 w12_w, G12_Gc, chi12_wGG = read_susceptibility_array('Al1_chiGG_q2.npz')
88 w22_w, G22_Gc, chi22_wGG = read_susceptibility_array('Al2_chiGG_q2.npz')
90 # Check that reciprocal lattice vectors remain as assumed in check below
91 assert np.linalg.norm(G11_Gc[0]) == pytest.approx(0., abs=1e-6)
92 assert np.linalg.norm(G21_Gc[0]) == pytest.approx(0., abs=1e-6)
93 assert np.linalg.norm(G12_Gc[0]) == pytest.approx(0., abs=1e-6)
94 assert np.linalg.norm(G22_Gc[1] - np.array([1, 0, 0])) == pytest.approx(
95 0., abs=1e-6)
97 # Check plasmon peaks remain the same
98 wpeak11, Ipeak11 = findpeak(w11_w, -chi11_wGG[:, 0, 0].imag)
99 wpeak21, Ipeak21 = findpeak(w21_w, -chi21_wGG[:, 0, 0].imag)
100 assert wpeak11 == pytest.approx(wpeak21, abs=0.02)
101 assert Ipeak11 == pytest.approx(Ipeak21, abs=1.0)
102 wpeak12, Ipeak12 = findpeak(w12_w, -chi12_wGG[:, 0, 0].imag)
103 wpeak22, Ipeak22 = findpeak(w22_w, -chi22_wGG[:, 1, 1].imag)
104 assert wpeak12 == pytest.approx(wpeak22, abs=0.05)
105 assert Ipeak12 == pytest.approx(Ipeak22, abs=1.0)
108def calculate_chi(calc, q_qc, w, nbands,
109 eta=0.2, ecut=50,
110 spincomponent='00', fxc=None,
111 filename_prefix=None, reduced_ecut=25):
112 gs = ResponseGroundStateAdapter(calc)
113 chiks_calc = ChiKSCalculator(gs, ecut=ecut, nbands=nbands)
114 chi_factory = ChiFactory(chiks_calc)
116 if filename_prefix is None:
117 filename = 'chiGG_qXXX.npz'
118 else:
119 filename = filename_prefix + '_chiGG_qXXX.npz'
121 for q, q_c in enumerate(q_qc):
122 fname = filename.replace('XXX', str(q + 1))
123 _, chi = chi_factory(spincomponent, q_c, w + 1.j * eta, fxc=fxc)
124 chi = chi.copy_with_reduced_ecut(reduced_ecut)
125 chi.write_array(fname)
126 world.barrier()