Coverage for gpaw/test/lcaotddft/test_molecule.py: 90%
214 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1import numpy as np
2import pytest
4from ase.utils import workdir
6from gpaw import GPAW
7from gpaw.mpi import world, serial_comm, broadcast
8from gpaw.lcaotddft.wfwriter import WaveFunctionReader
9from gpaw.lcaotddft.densitymatrix import DensityMatrix
10from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
11from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
13from gpaw.test import only_on_master
14from . import (parallel_options, calculate_time_propagation, calculate_error,
15 check_txt_data, check_wfs)
17pytestmark = pytest.mark.usefixtures('module_tmp_path')
19parallel_i = parallel_options(fix_sl_auto=True)
22@pytest.fixture(scope='module')
23def nacl_spin(gpw_files):
24 return gpw_files['nacl_spin']
27@pytest.fixture(scope='module')
28def nacl_nospin(gpw_files):
29 return gpw_files['nacl_nospin']
32@pytest.fixture(scope='module')
33@only_on_master(world)
34def initialize_system(nacl_nospin):
35 comm = serial_comm
36 calc = GPAW(nacl_nospin, communicator=comm)
37 fdm = calculate_time_propagation(nacl_nospin,
38 kick=np.ones(3) * 1e-5,
39 communicator=comm,
40 do_fdm=True)
42 # Calculate ground state with full unoccupied space
43 unocc_calc = calc.fixed_density(nbands='nao',
44 communicator=comm,
45 txt='unocc.out')
46 unocc_calc.write('unocc.gpw', mode='all')
47 return unocc_calc, fdm
50@pytest.mark.rttddft
51def test_propagated_wave_function(initialize_system, module_tmp_path):
52 wfr = WaveFunctionReader(module_tmp_path / 'wf.ulm')
53 coeff = wfr[-1].wave_functions.coefficients
54 # Pick a few coefficients corresponding to non-degenerate states;
55 # degenerate states should be normalized so that they can be compared
56 coeff = coeff[np.ix_([0], [0], [0, 1, 4], [0, 1, 2])]
57 # Normalize the wave function sign
58 coeff = np.sign(coeff.real[..., 0, np.newaxis]) * coeff
59 ref = [[[[1.6564776755628504e-02 + 1.2158943340143986e-01j,
60 4.7464497657284752e-03 + 3.4917799444496286e-02j,
61 8.2152048273399657e-07 - 1.6344333784831069e-06j],
62 [1.5177089239371724e-01 + 7.6502712023931621e-02j,
63 8.0497556154952932e-01 + 4.0573839188792121e-01j,
64 -5.1505952970811632e-06 - 1.1507918955641119e-05j],
65 [2.5116252101774323e+00 + 3.6776360873471503e-01j,
66 1.9024613198566329e-01 + 2.7843314959952882e-02j,
67 -1.3848736953929574e-05 - 2.6402210145403184e-05j]]]]
68 err = calculate_error(coeff, ref)
69 assert err < 1e-4
72@pytest.mark.rttddft
73@pytest.mark.parametrize('parallel', parallel_i)
74def test_propagation(initialize_system, module_tmp_path, parallel,
75 gpw_files, in_tmp_dir):
76 calculate_time_propagation(gpw_files['nacl_nospin'],
77 kick=np.ones(3) * 1e-5,
78 parallel=parallel)
79 check_wfs(module_tmp_path / 'wf.ulm', 'wf.ulm', atol=1e-12)
82@pytest.fixture(scope='module')
83@only_on_master(world, broadcast=broadcast)
84def dipole_moment_reference(initialize_system):
85 from gpaw.tddft.spectrum import \
86 read_dipole_moment_file, calculate_fourier_transform
88 unocc_calc, fdm = initialize_system
89 _, time_t, _, dm_tv = read_dipole_moment_file('dm.dat')
90 dm_tv = dm_tv - dm_tv[0]
91 dm_wv = calculate_fourier_transform(time_t, dm_tv,
92 fdm.foldedfreqs_f[0])
93 return dm_wv
96@pytest.fixture(scope='module')
97@only_on_master(world)
98def ksd_reference(initialize_system):
99 unocc_calc, fdm = initialize_system
100 ksd = KohnShamDecomposition(unocc_calc)
101 ksd.initialize(unocc_calc)
102 return ksd, fdm
105def ksd_transform_fdm(ksd, fdm):
106 rho_iwp = np.empty((2, len(fdm.freq_w), len(ksd.w_p)), dtype=complex)
107 rho_iwp[:] = np.nan + 1j * np.nan
108 for i, rho_wuMM in enumerate([fdm.FReDrho_wuMM, fdm.FImDrho_wuMM]):
109 for w in range(len(fdm.freq_w)):
110 rho_uMM = rho_wuMM[w]
111 rho_up = ksd.transform(rho_uMM)
112 rho_iwp[i, w, :] = rho_up[0]
113 return rho_iwp
116@pytest.fixture(scope='module')
117@only_on_master(world, broadcast=broadcast)
118def ksd_transform_reference(ksd_reference):
119 ksd, fdm = ksd_reference
120 ref_rho_iwp = ksd_transform_fdm(ksd, fdm)
121 return ref_rho_iwp
124@pytest.fixture(scope='module', params=parallel_i)
125def build_ksd(initialize_system, request):
126 calc = GPAW('unocc.gpw', parallel=request.param, txt=None)
127 if not calc.wfs.ksl.using_blacs and calc.wfs.bd.comm.size > 1:
128 pytest.xfail('Band parallelization without scalapack '
129 'is not supported')
130 ksd = KohnShamDecomposition(calc)
131 ksd.initialize(calc)
132 ksd.write('ksd.ulm')
135@pytest.fixture(scope='module', params=parallel_i)
136def load_ksd(build_ksd, request):
137 calc = GPAW('unocc.gpw', parallel=request.param, txt=None)
138 # Initialize positions in order to calculate density
139 calc.initialize_positions()
140 ksd = KohnShamDecomposition(calc, 'ksd.ulm')
141 dmat = DensityMatrix(calc)
142 fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
143 return ksd, fdm
146@pytest.fixture(scope='module')
147def ksd_transform(load_ksd):
148 ksd, fdm = load_ksd
149 rho_iwp = ksd_transform_fdm(ksd, fdm)
150 return rho_iwp
153@pytest.mark.skip(reason='See #933')
154@pytest.mark.rttddft
155def test_ksd_transform(ksd_transform, ksd_transform_reference):
156 ref_iwp = ksd_transform_reference
157 rho_iwp = ksd_transform
158 err = calculate_error(rho_iwp, ref_iwp)
159 atol = 1e-18
160 assert err < atol
163@pytest.mark.skip(reason='See #933')
164@pytest.mark.rttddft
165def test_ksd_transform_real_only(load_ksd, ksd_transform_reference):
166 ksd, fdm = load_ksd
167 ref_iwp = ksd_transform_reference
168 rho_iwp = np.empty((2, len(fdm.freq_w), len(ksd.w_p)), dtype=complex)
169 rho_iwp[:] = np.nan + 1j * np.nan
170 for i, rho_wuMM in enumerate([fdm.FReDrho_wuMM, fdm.FImDrho_wuMM]):
171 for w in range(len(fdm.freq_w)):
172 rho_uMM = rho_wuMM[w]
173 rho_p = ksd.transform([rho_uMM[0].real], broadcast=True)[0] \
174 + 1j * ksd.transform([rho_uMM[0].imag], broadcast=True)[0]
175 rho_iwp[i, w, :] = rho_p
176 err = calculate_error(rho_iwp, ref_iwp)
177 atol = 1e-18
178 assert err < atol
181@pytest.mark.rttddft
182def test_dipole_moment_from_ksd(ksd_transform, load_ksd,
183 dipole_moment_reference):
184 ksd, fdm = load_ksd
185 dm_wv = np.empty((len(fdm.freq_w), 3), dtype=complex)
186 dm_wv[:] = np.nan + 1j * np.nan
187 rho_wp = ksd_transform[0]
188 for w in range(len(fdm.freq_w)):
189 dm_v = ksd.get_dipole_moment([rho_wp[w]])
190 dm_wv[w, :] = dm_v
192 ref_wv = dipole_moment_reference
193 err = calculate_error(dm_wv, ref_wv)
194 atol = 1e-7
195 assert err < atol
198def get_density_fdm(ksd, fdm, kind):
199 assert kind in ['dmat', 'ksd']
200 rho_wg = fdm.dmat.density.finegd.empty(len(fdm.freq_w), dtype=complex)
201 rho_wg[:] = np.nan + 1j * np.nan
202 for w in range(len(fdm.freq_w)):
203 rho_uMM = fdm.FReDrho_wuMM[w]
204 if kind == 'dmat':
205 rho_g = fdm.dmat.get_density([rho_uMM[0].real]) \
206 + 1j * fdm.dmat.get_density([rho_uMM[0].imag])
207 elif kind == 'ksd':
208 rho_up = ksd.transform(rho_uMM, broadcast=True)
209 rho_g = ksd.get_density(fdm.dmat.wfs, [rho_up[0].real]) \
210 + 1j * ksd.get_density(fdm.dmat.wfs, [rho_up[0].imag])
211 rho_wg[w, :] = rho_g
212 return rho_wg
215@pytest.fixture(scope='module')
216@only_on_master(world, broadcast=broadcast)
217def density_reference(ksd_reference):
218 ksd, fdm = ksd_reference
219 dmat_rho_wg = get_density_fdm(ksd, fdm, 'dmat')
220 ksd_rho_wg = get_density_fdm(ksd, fdm, 'ksd')
221 return dict(dmat=dmat_rho_wg, ksd=ksd_rho_wg)
224@pytest.mark.rttddft
225def test_ksd_vs_dmat_density(density_reference):
226 ref_wg = density_reference['dmat']
227 rho_wg = density_reference['ksd']
228 err = calculate_error(rho_wg, ref_wg)
229 atol = 2e-10
230 assert err < atol
233@pytest.fixture(scope='module')
234def density(load_ksd):
235 ksd, fdm = load_ksd
236 if ksd.ksl.using_blacs:
237 pytest.xfail('Scalapack is not supported')
238 dmat_rho_wg = get_density_fdm(ksd, fdm, 'dmat')
239 ksd_rho_wg = get_density_fdm(ksd, fdm, 'ksd')
240 return dict(dmat=dmat_rho_wg, ksd=ksd_rho_wg)
243@pytest.mark.rttddft
244@pytest.mark.parametrize('kind', ['ksd', 'dmat'])
245def test_density(kind, density, load_ksd, density_reference):
246 ksd, fdm = load_ksd
247 ref_wg = density_reference[kind]
248 rho_wg = fdm.dmat.density.finegd.collect(density[kind])
249 err = calculate_error(rho_wg, ref_wg)
250 atol = 3e-19
251 assert err < atol
254@pytest.mark.rttddft
255@pytest.mark.parametrize('kind', ['ksd', 'dmat'])
256def test_dipole_moment_from_density(kind, density, load_ksd,
257 dipole_moment_reference):
258 ksd, fdm = load_ksd
259 rho_wg = density[kind]
260 dm_wv = np.empty((len(fdm.freq_w), 3), dtype=complex)
261 dm_wv[:] = np.nan + 1j * np.nan
262 for w in range(len(fdm.freq_w)):
263 dm_v = ksd.density.finegd.calculate_dipole_moment(rho_wg[w])
264 dm_wv[w, :] = dm_v
266 ref_wv = dipole_moment_reference
267 err = calculate_error(dm_wv, ref_wv)
268 atol = 5e-7
269 assert err < atol
272@pytest.mark.rttddft
273@only_on_master(world)
274def test_read_ksd(ksd_reference):
275 # Build a KohnShamDecomposition object from the calculator
276 ksd, _ = ksd_reference
278 # Now save it and read it without the calculator
279 ksd.write('ksd_save.ulm')
280 ksd_read = KohnShamDecomposition(filename='ksd_save.ulm')
282 np.testing.assert_equal(ksd.atoms, ksd_read.atoms)
284 for attr in ['S_uMM', 'C0_unM', 'eig_un', 'occ_un', 'C0S_unM']:
285 ref = getattr(ksd, attr)
286 test = getattr(ksd_read, attr)
288 np.testing.assert_almost_equal(ref, test)
291@pytest.fixture(scope='module')
292@only_on_master(world)
293@workdir('spinpol', mkdir=True)
294def initialize_system_spinpol(nacl_spin):
295 comm = serial_comm
296 calculate_time_propagation(nacl_spin,
297 kick=np.ones(3) * 1e-5,
298 communicator=comm,
299 do_fdm=True)
302@pytest.mark.rttddft
303def test_spinpol_dipole_moment(initialize_system, initialize_system_spinpol,
304 module_tmp_path):
305 # The test system has even number of electrons and is non-magnetic
306 # so spin-paired and spin-polarized calculation should give same result
307 check_txt_data(module_tmp_path / 'dm.dat',
308 module_tmp_path / 'spinpol' / 'dm.dat',
309 atol=1.0001e-12)
312@pytest.mark.rttddft
313@pytest.mark.parametrize('parallel', parallel_i)
314def test_spinpol_propagation(initialize_system_spinpol, module_tmp_path,
315 parallel, in_tmp_dir, gpw_files):
316 ref_path = module_tmp_path / 'spinpol'
317 calculate_time_propagation(gpw_files['nacl_spin'],
318 kick=np.ones(3) * 1e-5,
319 parallel=parallel)
320 check_wfs(ref_path / 'wf.ulm', 'wf.ulm', atol=1e-12)