Coverage for gpaw/test/lcaotddft/test_circular_dichroism.py: 100%
131 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import numpy as np
2import pytest
4from ase import Atoms
6from gpaw import GPAW
7from gpaw.mpi import world, serial_comm
8from gpaw.lcaotddft import LCAOTDDFT
9from gpaw.lcaotddft.densitymatrix import DensityMatrix
10from gpaw.lcaotddft.magneticmomentwriter import MagneticMomentWriter
11from gpaw.lcaotddft.magneticmomentwriter import parse_header
12from gpaw.tddft.spectrum import rotatory_strength_spectrum
13from gpaw.tddft.units import as_to_au, eV_to_au, au_to_eV, rot_au_to_cgs
15from gpaw.test import only_on_master
16from . import parallel_options, check_txt_data, copy_and_cut_file
18pytestmark = pytest.mark.usefixtures('module_tmp_path')
20parallel_i = parallel_options()
23# Parameterize over spin polarized and unpolarized calculations
24@pytest.fixture(scope='module', params=[True, False])
25@only_on_master(world)
26def initialize_system(request):
27 comm = serial_comm
29 atoms = Atoms('LiNaNaNa',
30 positions=[[0.0, 0.0, 0.0],
31 [2.0, 1.0, 0.0],
32 [4.0, 0.0, 1.0],
33 [6.0, -1.0, 0.0]])
34 atoms.center(vacuum=4.0)
36 if request.param:
37 atoms.set_initial_magnetic_moments([0.01] * len(atoms))
39 calc = GPAW(nbands=2,
40 h=0.4,
41 setups={'Na': '1'},
42 basis='sz(dzp)',
43 mode='lcao',
44 convergence={'density': 1e-12},
45 communicator=comm,
46 symmetry={'point_group': False},
47 txt='gs.out')
48 atoms.calc = calc
49 atoms.get_potential_energy()
50 calc.write('gs.gpw', mode='all')
52 for gauge in ['length', 'velocity']:
53 if gauge == 'velocity':
54 add = '_' + gauge
55 else:
56 add = ''
57 td_calc = LCAOTDDFT('gs.gpw', communicator=comm,
58 txt='td' + add + '.out')
59 dmat = DensityMatrix(td_calc)
60 MagneticMomentWriter(td_calc, 'mm' + add + '.dat', dmat=dmat)
61 MagneticMomentWriter(td_calc, 'mm_grid' + add + '.dat',
62 calculate_on_grid=True)
63 MagneticMomentWriter(td_calc, 'mm_origin' + add + '.dat',
64 origin='zero', origin_shift=[1.0, 2.0, 3.0])
65 td_calc.absorption_kick([1e-5, 0., 0.], gauge=gauge)
66 td_calc.propagate(100, 3)
67 td_calc.write('td' + add + '.gpw', mode='all')
68 td_calc.propagate(100, 2)
71@pytest.mark.rttddft
72def test_magnetic_moment_velocity_gauge(initialize_system, module_tmp_path,
73 in_tmp_dir):
74 with open('mm_ref.dat', 'w', encoding='utf-8') as f:
75 f.write('''
76# MagneticMomentWriter[version=5](**{"origin": "COM", "origin_shift": [0.0, 0.0, 0.0], "calculate_on_grid": false, "only_pseudo": false})
77# origin_v = [7.634300, 5.000000, 4.302858] Å
78# time mmx mmy mmz
79# Start; Time = 0.00000000
80 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
81# Kick = [ 1.000000000000e-05, 0.000000000000e+00, 0.000000000000e+00]; Time = 0.00000000
82 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
83 4.13413733 -3.063974663768e-07 -3.574171353385e-07 1.245812807572e-06
84 8.26827467 -1.197667430552e-06 -1.414585723558e-06 4.779431710530e-06
85 12.40241200 -2.606870392812e-06 -3.132296132758e-06 1.010969895228e-05
86 16.53654934 -4.424269197168e-06 -5.428714567811e-06 1.652651689986e-05
87 20.67068667 -6.501209589454e-06 -8.162871785725e-06 2.323149911180e-05
88'''.strip()) # noqa: E501
89 check_txt_data(module_tmp_path / 'mm_velocity.dat',
90 'mm_ref.dat', atol=2e-14)
92 with open('mm_grid_ref.dat', 'w', encoding='utf-8') as f:
93 f.write('''
94# MagneticMomentWriter[version=5](**{"origin": "COM", "origin_shift": [0.0, 0.0, 0.0], "calculate_on_grid": true, "only_pseudo": false})
95# origin_v = [7.634300, 5.000000, 4.302858] Å
96# time mmx mmy mmz
97# Start; Time = 0.00000000
98 0.00000000 -0.000000000000e+00 -0.000000000000e+00 -0.000000000000e+00
99# Kick = [ 1.000000000000e-05, 0.000000000000e+00, 0.000000000000e+00]; Time = 0.00000000
100 0.00000000 -0.000000000000e+00 -0.000000000000e+00 -0.000000000000e+00
101 4.13413733 -3.063224380042e-07 -3.573133338861e-07 1.245412048950e-06
102 8.26827467 -1.197361494602e-06 -1.414167218866e-06 4.777891626962e-06
103 12.40241200 -2.606162974554e-06 -3.131346591800e-06 1.010643346002e-05
104 16.53654934 -4.422978770321e-06 -5.427020998550e-06 1.652116134041e-05
105 20.67068667 -6.499162446085e-06 -8.160246795614e-06 2.322393973892e-05
106'''.strip()) # noqa: E501
108 check_txt_data(module_tmp_path / 'mm_grid_velocity.dat',
109 'mm_grid_ref.dat', atol=2e-14)
111 with open('mm_origin_ref.dat', 'w', encoding='utf-8') as f:
112 f.write('''
113# MagneticMomentWriter[version=5](**{"origin": "zero", "origin_shift": [1.0, 2.0, 3.0], "calculate_on_grid": false, "only_pseudo": false})
114# origin_v = [1.000000, 2.000000, 3.000000] Å
115# time mmx mmy mmz
116# Start; Time = 0.00000000
117 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
118# Kick = [ 1.000000000000e-05, 0.000000000000e+00, 0.000000000000e+00]; Time = 0.00000000
119 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
120 4.13413733 3.055448269266e-07 -2.985624975165e-06 4.181522917320e-06
121 8.26827467 1.055413708309e-06 -1.144251258824e-05 1.639709003121e-05
122 12.40241200 1.830150711621e-06 -2.413141874059e-05 3.586908883647e-05
123 16.53654934 2.098630938582e-06 -3.928697753788e-05 6.127421238813e-05
124 20.67068667 1.409938523718e-06 -5.497359580260e-05 9.073479845925e-05
125'''.strip()) # noqa: E501
127 check_txt_data(module_tmp_path / 'mm_origin_velocity.dat',
128 'mm_origin_ref.dat', atol=1e-13)
131@pytest.mark.rttddft
132def test_magnetic_moment_values(initialize_system, module_tmp_path,
133 in_tmp_dir):
134 with open('mm_ref.dat', 'w', encoding='utf-8') as f:
135 f.write('''
136# MagneticMomentWriter[version=4](origin='COM')
137# origin_v = [7.634300, 5.000000, 4.302858] Å
138# time mmx mmy mmz
139 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
140# Kick = [ 1.000000000000e-05, 0.000000000000e+00, 0.000000000000e+00]; Time = 0.00000000
141 0.00000000 8.192189793082e-06 1.038446327373e-05 -2.730498071751e-05
142 4.13413733 7.838837723234e-06 1.000765310013e-05 -2.573300722038e-05
143 8.26827467 6.809084660174e-06 8.879683492897e-06 -2.128890950807e-05
144 12.40241200 5.175350632237e-06 7.009694921954e-06 -1.462938416394e-05
145 16.53654934 3.058296873929e-06 4.443905967036e-06 -6.697375210691e-06
146 20.67068667 6.247451722277e-07 1.298788405738e-06 1.460017881082e-06
147'''.strip()) # noqa: E501
149 check_txt_data(module_tmp_path / 'mm.dat', 'mm_ref.dat', atol=1e-13)
152@pytest.mark.rttddft
153def test_magnetic_moment_grid_evaluation(initialize_system, module_tmp_path):
154 dpath = module_tmp_path
155 check_txt_data(dpath / 'mm.dat', dpath / 'mm_grid.dat', atol=2e-8)
158@pytest.mark.rttddft
159@pytest.mark.parametrize('parallel', parallel_i)
160@pytest.mark.parametrize('gauge', ['velocity', 'length'])
161def test_magnetic_moment_parallel(initialize_system, module_tmp_path, parallel,
162 in_tmp_dir, gauge):
164 td_calc = LCAOTDDFT(module_tmp_path / 'gs.gpw',
165 parallel=parallel,
166 txt='td.out')
167 print(parallel)
168 add = '_velocity' if gauge == 'velocity' else ''
169 MagneticMomentWriter(td_calc, f'mm{add}.dat')
170 MagneticMomentWriter(td_calc, f'mm_grid{add}.dat', calculate_on_grid=True)
171 MagneticMomentWriter(td_calc, f'mm_origin{add}.dat',
172 origin='zero', origin_shift=[1.0, 2.0, 3.0])
173 td_calc.absorption_kick([1e-5, 0., 0.], gauge=gauge)
174 td_calc.propagate(100, 5)
176 for fname in [f'mm{add}.dat', f'mm_grid{add}.dat', f'mm_origin{add}.dat']:
177 check_txt_data(module_tmp_path / fname, fname, atol=7e-14)
180@pytest.mark.rttddft
181@pytest.mark.parametrize('parallel', parallel_i)
182def test_magnetic_moment_restart(initialize_system, module_tmp_path, parallel,
183 in_tmp_dir):
184 td_calc = LCAOTDDFT(module_tmp_path / 'td.gpw',
185 parallel=parallel,
186 txt='td.out')
187 for fname in ['mm.dat', 'mm_grid.dat', 'mm_origin.dat']:
188 if world.rank == 0:
189 copy_and_cut_file(module_tmp_path / fname, fname, cut_lines=3)
190 world.barrier()
191 MagneticMomentWriter(td_calc, fname)
192 td_calc.propagate(100, 2)
194 for fname in ['mm.dat', 'mm_grid.dat', 'mm_origin.dat']:
195 check_txt_data(module_tmp_path / fname, fname, atol=7e-14)
198@only_on_master(world)
199def test_spectrum(in_tmp_dir, rng):
200 from gpaw.utilities.folder import Folder
202 # Parameters for test data
203 kick_strength = 1e-5
204 frequency_v = np.array([1.0, 2.0, 3.0]) * eV_to_au
205 strength_v = np.array([1.0, 2.0, 3.0])
207 # Create dummy magnetic moment files
208 for v, kick in enumerate('xyz'):
209 kick_v = [0.0, 0.0, 0.0]
210 kick_v[v] = kick_strength
211 time_t = np.arange(0, 31e3, 10.0) * as_to_au
212 data_tv = np.zeros((len(time_t), 4))
213 data_tv[:, 0] = time_t
214 # Fill unused columns with random values
215 data_tv[:, 1:] = rng.random((len(time_t), 3))
216 # Diagonal column has the data used for spectrum
217 data_tv[:, v + 1] = (kick_strength * strength_v[v]
218 * np.cos(frequency_v[v] * time_t))
219 with open(f'mm-{kick}.dat', 'w', encoding='utf-8') as f:
220 f.write(f'''
221# MagneticMomentWriter[version=4](origin='COM')
222# time mmx mmy mmz
223 0.00000000 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
224# Kick = {kick_v}; Time = 0.00000000
225''') # noqa: E501
226 np.savetxt(f, data_tv)
228 # Calculate spectrum
229 folding_kwargs = dict(folding='Gauss', width=0.2)
230 rotatory_strength_spectrum(['mm-x.dat', 'mm-y.dat', 'mm-z.dat'],
231 'spec.dat',
232 **folding_kwargs,
233 e_min=0.0, e_max=5.0, delta_e=0.01)
235 # Reference spectrum
236 energy_e = np.arange(0, 5.0 + 1e-8, 0.01)
237 f = Folder(**folding_kwargs).fold_values
238 spec_e = (f(frequency_v * au_to_eV, strength_v, energy_e)[1]
239 + f(-frequency_v * au_to_eV, strength_v, energy_e)[1])
240 spec_e *= 0.5
241 spec_e *= rot_au_to_cgs * 1e40
243 # Compare spectrum
244 data_ei = np.loadtxt('spec.dat')
245 assert np.allclose(data_ei[:, 0], energy_e)
246 assert np.allclose(data_ei[:, 1], spec_e)
248 # Test failure
249 with pytest.raises(RuntimeError):
250 rotatory_strength_spectrum(['mm-x.dat', 'mm-x.dat', 'mm-z.dat'],
251 'spec.dat')
252 with pytest.raises(RuntimeError):
253 rotatory_strength_spectrum(['mm-y.dat', 'mm-x.dat', 'mm-z.dat'],
254 'spec.dat')
255 with pytest.raises(RuntimeError):
256 rotatory_strength_spectrum(['mm-x.dat', 'mm-y.dat', 'mm-x.dat'],
257 'spec.dat')
260def test_parse_header():
261 line = 'SomeWriter[version=42](**{"str": "value", "vector": [1.2, 3.4], "boolean": true})' # noqa: E501
262 name, version, kwargs = parse_header(line)
263 assert name == 'SomeWriter'
264 assert version == 42
265 assert kwargs['str'] == 'value'
266 assert kwargs['vector'] == [1.2, 3.4]
267 assert kwargs['boolean']
269 # Test failure
270 with pytest.raises(ValueError):
271 parse_header('wrong line')
272 with pytest.raises(ValueError):
273 parse_header('A[version=1](**{wrong json})')