Octopus
tdtdm.F90
Go to the documentation of this file.
1!! Copyright (C) 2019-2021 N. Tancogne-Dejean
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21program tdtdm
22 use batch_oct_m
24 use comm_oct_m
25 use debug_oct_m
28 use fft_oct_m
29 use global_oct_m
30 use grid_oct_m
32 use io_oct_m
37 use mesh_oct_m
39 use mpi_oct_m
42 use parser_oct_m
51 use types_oct_m
52 use unit_oct_m
54 use xc_oct_m
55
56 implicit none
57
58 integer :: in_file, ii, jj, kk, ierr, ip_h, irow, ifreq, nrow, it
59 integer :: ik, ist, uist, istep, ikpoint, irep, out_file, iop, idim
60 integer :: time_steps, energy_steps, istart, iend, ntiter, Nreplica, Ntrans
61 real(real64) :: dt, tt, weight, kpoint(3), kpoint_sym(3), kred(3), kred_sym(3)
62 real(real64) :: xx_h_sym(3)
63 integer :: irep_h, ip_h_sym, rankmin
64 real(real64) :: start_time, dmin
65 real(real64), allocatable :: Et(:), ftreal(:, :, :), ftimag(:, :, :), tmp(:), omega(:)
66 complex(real64), allocatable :: Xiak(:,:,:), Yiak(:,:,:)
67 real(real64), allocatable :: proj_r(:,:,:,:), proj_i(:,:,:,:)
68 real(real64), allocatable :: proj_r_corr(:,:), proj_i_corr(:,:), centers(:,:)
69 complex(real64), allocatable :: tdm(:,:), tdm_1D(:,:,:,:)
70 complex(real64), allocatable, target :: psi(:,:), upsi(:,:)
71 complex(real64), allocatable :: phase(:,:,:), ftcmplx(:,:)
72 complex(real64), pointer :: psi_sym(:,:), upsi_sym(:,:)
73 type(spectrum_t) :: spectrum
74 type(electrons_t), pointer :: sys
75 type(batch_t) :: projb_r, projb_i, ftrealb, ftimagb
76 character(len=MAX_PATH_LEN) :: fname
77 type(states_elec_t), pointer :: st
78 type(states_elec_t) :: gs_st
79 type(restart_t) :: restart
80 type(unit_t) :: fn_unit
81 integer :: kpt_start, kpt_end, supercell(3), nomega, ncols
82 type(block_t) :: blk
83 real(real64) :: pos_h(3), norm
84
85 ! Initializion
86 call global_init()
87 call parser_init()
88
89 call messages_init()
90 call io_init()
91
93
94 call messages_experimental("oct-tdtdm utility")
98
99 call calc_mode_par%set_parallelization(p_strategy_states, default = .false.)
101 call sys%init_parallelization(mpi_world)
102
103 call spectrum_init(spectrum, global_namespace)
104
105 st => sys%st
106
107 if(sys%st%d%ispin == spinors) then
108 call messages_not_implemented('oct-tdtdm with spinors')
109 end if
110
111 if(st%parallel_in_states) then
112 call messages_not_implemented("oct-tdtdm with states parallelization")
113 end if
115 if(sys%gr%parallel_in_domains) then
116 call messages_not_implemented("oct-tdtdm with domain parallelization")
117 end if
118
119 !%Variable TDTDMFrequencies
120 !%Type block
121 !%Section Utilities::oct-tdtdm
122 !%Description
123 !% This block defines for which frequencies the analysis is performed.
124 !%
125 !% Each row of the block indicates a frequency.
126 !%End
127 if (parse_block(global_namespace, 'TDTDMFrequencies', blk) == 0) then
128
129 nrow = parse_block_n(blk)
130 nomega = nrow
131
132 safe_allocate(omega(1:nrow))
133 !read frequencies
134 do irow = 0, nrow-1
135 call parse_block_float(blk, irow, 0, omega(irow+1))
136 end do
137
138 call parse_block_end(blk)
139 else
140 message(1) = "oct-tdtdm: TDTDMFrequencies must be defined."
141 call messages_fatal(1)
142 end if
143
144 ! We check that the resonant and antiresonant transitions are contained in the
145 ! energy range of the Fourier transforms
146 if(any(omega > spectrum%max_energy)) then
147 message(1) = "One requested frequecy is larger than PropagationSpectrumMaxEnergy."
148 message(2) = "Please increase the value of PropagationSpectrumMaxEnergy."
149 call messages_fatal(2)
150 end if
151 if(any(omega > -spectrum%min_energy)) then
152 message(1) = "One requested frequency is larger than -PropagationSpectrumMinEnergy."
153 message(2) = "Please decrease the value of PropagationSpectrumMinEnergy."
154 call messages_fatal(2)
155 end if
156
157
158 call states_elec_copy(gs_st, st, exclude_wfns = .true., exclude_eigenval = .true.)
159
160 safe_deallocate_a(gs_st%node)
161
162 call restart_init(restart, global_namespace, restart_proj, restart_type_load, sys%mc, ierr, mesh=sys%gr)
163 if(ierr == 0) call states_elec_look(restart, ii, jj, gs_st%nst, ierr)
164 if(ierr /= 0) then
165 message(1) = "oct-tdtdm: Unable to read states information."
166 call messages_fatal(1)
167 end if
168
169 ! allocate memory
170 safe_allocate(gs_st%occ(1:gs_st%nst, 1:gs_st%nik))
171 safe_allocate(gs_st%eigenval(1:gs_st%nst, 1:gs_st%nik))
172
173 ! We want all the task to have all the states
174 ! States can be distibuted for the states we propagate.
175 safe_allocate(gs_st%node(1:gs_st%nst))
176 gs_st%node(:) = 0
177 call kpoints_distribute(gs_st, sys%mc)
179
180 kpt_start = gs_st%d%kpt%start
181 kpt_end = gs_st%d%kpt%end
182
183 gs_st%eigenval = huge(gs_st%eigenval)
184 gs_st%occ = m_zero
185 if(gs_st%d%ispin == spinors) then
186 safe_deallocate_a(gs_st%spin)
187 safe_allocate(gs_st%spin(1:3, 1:gs_st%nst, 1:gs_st%nik))
188 end if
189
190 call states_elec_allocate_wfns(gs_st, sys%gr, type_cmplx)
191 call states_elec_load(restart, global_namespace, sys%space, gs_st, sys%gr, sys%kpoints, ierr)
192 if(ierr /= 0 .and. ierr /= (gs_st%st_end-gs_st%st_start+1)*(kpt_end-kpt_start+1)*gs_st%d%dim) then
193 message(1) = "oct-tdtdm: Unable to read wavefunctions for TDOutput."
194 call messages_fatal(1)
195 end if
196 call restart_end(restart)
197
198
199 in_file = io_open('td.general/projections', action='read', status='old')
200 call io_skip_header(in_file)
201 call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
202 dt = units_to_atomic(units_out%time, dt)
203
204
205 safe_allocate(tmp(1:st%nst*gs_st%nst*st%nik*2))
206 safe_allocate(proj_r(1:time_steps, 1:gs_st%nst, 1:st%nst, 1:st%nik))
207 safe_allocate(proj_i(1:time_steps, 1:gs_st%nst, 1:st%nst, 1:st%nik))
208
209
210 call io_skip_header(in_file)
211
212 do ii = 1, time_steps
213 read(in_file, *) jj, tt, (tmp(kk), kk = 1, st%nst*gs_st%nst*st%nik*2)
214 do ik = 1, st%nik
215 do ist = 1, st%nst
216 do uist = 1, gs_st%nst
217 jj = (ik-1)*st%nst*gs_st%nst + (ist-1)*gs_st%nst + uist
218 proj_r(ii, uist, ist, ik) = tmp((jj-1)*2+1)
219 ! Here we add a minus sign, as we want to get <\phi_0 | \psi(t)>
220 ! and td_occup computes the complex conjugaute of this
221 proj_i(ii, uist, ist, ik) = -tmp((jj-1)*2+2)
222 end do
223 end do
224 end do
225 end do
226 safe_deallocate_a(tmp)
227
228 call io_close(in_file)
229
230 write(message(1), '(a, i7, a)') "oct-tdtdm: Read ", time_steps, " steps from file '"// &
231 trim(io_workpath('td.general/projections'))//"'"
232 call messages_info(1)
233
234 start_time = spectrum%start_time
235
236 ! Phase correction of the projections before doing the Fourier transforms
237 ! See Eq. (5) of Williams et al., JCTC 17, 1795 (2021)
238 ! We need to multiply C_ik(t)e^{-ie_kt} (the projection of \phi_i(t) on \phi_k^GS)
239 ! by e^{ie_it}, which is obtained by the cc of the projection of \phi_i(t) on \phi_i^GS
240 ! Here we only care about optical transitions (so TD occupied to GS unocc)
241 safe_allocate(proj_r_corr(1:time_steps, 1:gs_st%nst*st%nst*(kpt_end-kpt_start+1)))
242 safe_allocate(proj_i_corr(1:time_steps, 1:gs_st%nst*st%nst*(kpt_end-kpt_start+1)))
243 proj_r_corr = m_zero
244 proj_i_corr = m_zero
245 do ik = kpt_start, kpt_end
246 do ist = 1, st%nst
247 do uist = ist+1, gs_st%nst
248 jj = (ik-kpt_start)*st%nst*gs_st%nst+(ist-1)*gs_st%nst+uist
249 do ii = 1, time_steps
250 norm = hypot(proj_r(ii, ist, ist, ik),proj_i(ii, ist, ist, ik))
251 proj_r_corr(ii, jj) = (proj_r(ii, uist, ist, ik) * proj_r(ii, ist, ist, ik) &
252 + proj_i(ii, uist, ist, ik) * proj_i(ii, ist, ist, ik))/norm
253 proj_i_corr(ii, jj) =(-proj_r(ii, uist, ist, ik) * proj_i(ii, ist, ist, ik) &
254 + proj_i(ii, uist, ist, ik) * proj_r(ii, ist, ist, ik))/norm
255 end do
256 end do
257 end do
258 end do
259
260 safe_deallocate_a(proj_r)
261 safe_deallocate_a(proj_i)
262
263 ! Find out the iteration numbers corresponding to the time limits.
264 call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
265 istart = max(1, istart)
266 energy_steps = spectrum_nenergy_steps(spectrum)
267
268 safe_allocate(ftreal(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1), 1:2))
269 safe_allocate(ftimag(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1), 1:2))
270
271 call batch_init(projb_r, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), proj_r_corr)
272 call batch_init(projb_i, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), proj_i_corr)
273 call batch_init(ftrealb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftreal(:,:,1))
274 call batch_init(ftimagb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftimag(:,:,1))
275
276 write(message(1), '(a)') "oct-tdtdm: Fourier transforming real part of the projections"
277 call messages_info(1)
278
279 call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
280 istart, iend, spectrum%start_time, dt, projb_r, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftrealb)
281
282 call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
283 istart, iend, spectrum%start_time, dt, projb_r, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftimagb)
284
285 call ftrealb%end()
286 call ftimagb%end()
287
288 safe_allocate(ftcmplx(1:energy_steps, 1:st%nst*gs_st%nst*(kpt_end-kpt_start+1)))
289 do ii = 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1)
290 ftcmplx(1:energy_steps,ii) = cmplx(ftreal(1:energy_steps,ii,1), ftimag(1:energy_steps,ii,1), real64)
291 end do
292
293 write(message(1), '(a)') "oct-tdtdm: Fourier transforming imaginary part of the projections"
294 call messages_info(1)
295
296 call batch_init(ftrealb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftreal(:,:,2))
297 call batch_init(ftimagb, 1, 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1), ftimag(:,:,2))
298
299 call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
300 istart, iend, spectrum%start_time, dt, projb_i, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftrealb)
301
302 call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
303 istart, iend, spectrum%start_time, dt, projb_i, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftimagb)
304
305 call projb_i%end()
306 call projb_r%end()
307 call ftrealb%end()
308 call ftimagb%end()
309 safe_deallocate_a(proj_r_corr)
310 safe_deallocate_a(proj_i_corr)
311
312 do ii = 1, st%nst*gs_st%nst*(kpt_end-kpt_start+1)
313 ftcmplx(1:energy_steps,ii) = ftcmplx(1:energy_steps,ii) + m_zi*ftreal(1:energy_steps,ii,2) - ftimag(1:energy_steps,ii,2)
314 end do
315
316 safe_deallocate_a(ftreal)
317 safe_deallocate_a(ftimag)
318
319 write(message(1), '(a)') "oct-tdtdm: Constructing the two-particle wavefunctions."
320 call messages_info(1)
321
322 !%Variable SupercellDimensions
323 !%Type block
324 !%Default KPointsGrid
325 !%Section Utilities::oct-tdtdm
326 !%Description
327 !% This block allows to specify the size of the supercell used to plot excitonic wavefunctions.
328 !% If not specified, the code uses the number of k-points for defining the size of the supercell.
329 !%End
330 if (parse_is_defined(sys%namespace, 'SupercellDimensions')) then
331 if (parse_block(sys%namespace, 'SupercellDimensions', blk) == 0) then
332 ncols = parse_block_cols(blk, 0)
333 if (ncols /= sys%space%dim) then
334 write(message(1),'(a,i3,a,i3)') 'SupercellDimensions has ', ncols, ' columns but must have ', sys%space%dim
335 call messages_fatal(1, namespace=sys%namespace)
336 end if
337 do ii = 1, sys%space%dim
338 call parse_block_integer(blk, 0, ii - 1, supercell(ii))
339 end do
340
341 call parse_block_end(blk)
342 end if
343 else
344 supercell(1:sys%space%dim) = sys%kpoints%nik_axis(1:sys%space%dim)
345 end if
346
347 nreplica = product(supercell(1:sys%space%dim))
348
349 ! The center of each replica of the unit cell
350 safe_allocate(centers(1:sys%space%dim, 1:nreplica))
351 irep = 1
352 do ii = 0, supercell(1)-1
353 do jj = 0, supercell(2)-1
354 do kk = 0, supercell(3)-1
355 centers(1, irep) = -floor((supercell(1)-1)/m_two)+ii
356 centers(2, irep) = -floor((supercell(2)-1)/m_two)+jj
357 centers(3, irep) = -floor((supercell(3)-1)/m_two)+kk
358 centers(:, irep) = matmul(sys%ions%latt%rlattice, centers(:, irep))
359 irep = irep + 1
360 end do
361 end do
362 end do
363
364 ! The phase for each center
365 irep = 0
366 do ik = kpt_start, kpt_end
367 ikpoint = gs_st%d%get_kpoint_index(ik)
368 irep = max(irep, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint))
369 end do
370 safe_allocate(phase(kpt_start:kpt_end, 1:irep, 1:nreplica))
371 do irep = 1, nreplica
372 do ik = kpt_start, kpt_end
373 ikpoint = gs_st%d%get_kpoint_index(ik)
374 kpoint(1:sys%space%dim) = sys%kpoints%get_point(ikpoint)
375 do ii = 1, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
376 iop = kpoints_get_symmetry_ops(sys%kpoints, ikpoint, ii)
377
378 if (sys%kpoints%use_symmetries) then !We apply the symmetry
379 call kpoints_to_reduced(sys%kpoints%latt, kpoint, kred)
380 call symmetries_apply_kpoint_red(sys%kpoints%symm, iop, kred, kred_sym)
381 call kpoints_to_absolute(sys%kpoints%latt, kred_sym, kpoint_sym)
382 else
383 kpoint_sym = kpoint
384 end if
385 phase(ik, ii, irep) = exp(-m_zi*sum(kpoint_sym(1:sys%space%dim)*centers(:, irep)))
386 end do
387 end do
388 end do
389
390 ! Position of the hole, here assumed to be on top of the first atom
391 ! To be obtained from the input file
392 if(sys%space%dim > 1) then
393 call tdtdm_get_hole_position(pos_h, ip_h)
394 end if
395
396 ntrans = 0
397 ! Here we assume that there is a clear gap, so the information at Gamma is enough
398 do ist = 1, gs_st%nst
399 if(abs(gs_st%occ(ist, 1)) < m_epsilon) cycle
400
401 do uist = 1, gs_st%nst
402 if(abs(gs_st%occ(uist, 1)) > m_epsilon) cycle
403 weight = gs_st%kweights(1) * (gs_st%occ(ist, 1)-gs_st%occ(uist, 1))
404 if(abs(weight) < m_epsilon) cycle
405 ntrans = ntrans + 1
406 end do
407 end do
408 if(ntrans == 0) then
409 write(message(1), '(a)') "oct-tdtdm: No transition found."
410 write(message(2), '(a)') "Please check that unoccupied states are included in the ground state calculation."
411 call messages_fatal(2)
412 end if
413
414 safe_allocate(xiak(1:st%nst, 1:gs_st%nst, 1:st%nik))
415 safe_allocate(yiak(1:st%nst, 1:gs_st%nst, 1:st%nik))
416 safe_allocate(et(1:ntrans*st%nik))
417 safe_allocate(psi(1:sys%gr%np, 1:gs_st%d%dim))
418 safe_allocate(upsi(1:sys%gr%np, 1:gs_st%d%dim))
419
420 if(sys%kpoints%use_symmetries) then
421 safe_allocate(psi_sym(1:sys%gr%np, 1:st%d%dim))
422 safe_allocate(upsi_sym(1:sys%gr%np, 1:st%d%dim))
423 end if
424
425 select case(sys%space%dim)
426 case(2,3)
427 safe_allocate(tdm(1:sys%gr%np, 1:nreplica))
428 case(1)
429 safe_allocate(tdm_1d(1:sys%gr%np, 1:sys%gr%np, 1:nreplica, 1:nreplica))
430 end select
431
432 do ifreq = 1, nomega
433
434 write(message(1), '(a, f6.4, a)') "oct-tdtdm: Constructing the two-particle wavefunction at ", omega(ifreq), " Ha."
435 call messages_info(1)
436
437 select case(sys%space%dim)
438 case(2,3)
439 tdm = m_z0
440 case(1)
441 tdm_1d = m_z0
442 end select
443
444 et = m_zero
445 xiak = m_z0
446 yiak = m_z0
447
448 ! Local transition index
449 it = (kpt_start-1)*ntrans + 1
450
451 do ik = kpt_start, kpt_end
452 ikpoint = st%d%get_kpoint_index(ik)
453
454 do ist = 1, st%nst
455 if(abs(gs_st%occ(ist, ik)) < m_epsilon) cycle
456
457 call states_elec_get_state(gs_st, sys%gr, ist, ik, psi)
458 if (sys%hm%phase%is_allocated()) then
459 call sys%hm%phase%apply_to_single(psi, sys%gr%np, gs_st%d%dim, ik, .false.)
460 end if
461
462 do uist = 1, gs_st%nst
463 if(abs(gs_st%occ(uist, ik)) > m_epsilon) cycle
464
465 ! For a given requested frequency, we get the corresponding values of Xia and Yia
466 ! One correspond to the +\Omega frequency, the other one to the -\Omega frequency
467 ! For Xiak, we use the fact that TF[f*](\Omega) = (TF[f](-\Omega))^*
468 jj = (ik-kpt_start)*st%nst*gs_st%nst+(ist-1)*gs_st%nst+uist
469 istep = int((+omega(ifreq)-spectrum%min_energy)/spectrum%energy_step)
470 xiak(ist, uist, ik) = conjg(ftcmplx(istep, jj))
471 istep = int((+omega(ifreq)-spectrum%min_energy)/spectrum%energy_step)
472 yiak(ist, uist, ik) = ftcmplx(istep, jj)
473
474
475 weight = gs_st%kweights(ik) * (gs_st%occ(ist, ik)-gs_st%occ(uist, ik)) &
476 / kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
477 if(abs(weight) < m_epsilon) cycle
478
479 call states_elec_get_state(gs_st, sys%gr, uist, ik, upsi)
480 if(sys%hm%phase%is_allocated()) then
481 call sys%hm%phase%apply_to_single(upsi, sys%gr%np, st%d%dim, ik, .false.)
482 end if
483
484 do ii = 1, kpoints_get_num_symmetry_ops(sys%kpoints, ikpoint)
485 iop = kpoints_get_symmetry_ops(sys%kpoints, ikpoint, ii)
486
487 if(sys%kpoints%use_symmetries) then
488 do idim = 1, st%d%dim
489 call zgrid_symmetrize_single(sys%gr, iop, psi(:,idim), psi_sym(:,idim))
490 call zgrid_symmetrize_single(sys%gr, iop, upsi(:,idim), upsi_sym(:,idim))
491 end do
492
493 ! We need to get the position of the hole after applying the symmetry operation too
494 xx_h_sym = symm_op_apply_cart(sys%kpoints%symm%ops(iop), pos_h)
495 xx_h_sym = sys%ions%latt%fold_into_cell(xx_h_sym)
496 ! At the moment, we ignore rankmin
497 assert(.not.sys%gr%parallel_in_domains)
498 ip_h_sym = mesh_nearest_point(sys%gr, xx_h_sym, dmin, rankmin)
499 else
500 psi_sym => psi
501 upsi_sym => upsi
502 ip_h_sym = ip_h
503 end if
504
505 ! We now compute the single mode TDTDM
506 ! See Eq. (5) of Williams et al., JCTC 17, 1795 (2021)
507 ! We take here the complex conjugate of the 2-body wavefunction
508 select case(sys%space%dim)
509 case(2,3)
510 do irep = 1, nreplica
511 call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * weight &
512 * conjg(xiak(ist,uist,ik))*conjg(psi_sym(ip_h_sym,1)), upsi_sym(:, 1), tdm(:,irep))
513 call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * weight &
514 * yiak(ist,uist,ik)*conjg(upsi_sym(ip_h_sym,1)), psi_sym(:, 1), tdm(:,irep))
515 end do
516 case(1)
517 ! In the 1D case, we contruct the full TDTDM of r_e, r_h
518 do irep_h = 1, nreplica
519 do irep = 1, nreplica
520 do ip_h = 1, sys%gr%np
521 call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * conjg(phase(ik, ii, irep_h)) &
522 * weight * conjg(xiak(ist,uist,ik)) * conjg(psi_sym(ip_h,1)), &
523 upsi_sym(:, 1), tdm_1d(:, ip_h, irep, irep_h))
524 call lalg_axpy(sys%gr%np, phase(ik, ii, irep) * conjg(phase(ik, ii, irep_h)) &
525 * weight * conjg(yiak(ist,uist,ik)) * conjg(upsi_sym(ip_h,1)), &
526 psi_sym(:, 1), tdm_1d(:, ip_h, irep, irep_h))
527 end do
528 end do
529 end do
530 end select
531
532 end do ! ii
533
534 et(it) = gs_st%eigenval(uist, ik) - gs_st%eigenval(ist, ik)
535 it = it + 1
536 end do
537 end do
538 end do
539
540 if(gs_st%d%kpt%parallel) then
541 if(sys%space%dim > 1) then
542 call comm_allreduce(gs_st%d%kpt%mpi_grp, tdm)
543 else
544 call comm_allreduce(gs_st%d%kpt%mpi_grp, tdm_1d)
545 end if
546 call comm_allreduce(gs_st%d%kpt%mpi_grp, et)
547 call comm_allreduce(gs_st%d%kpt%mpi_grp, xiak)
548 call comm_allreduce(gs_st%d%kpt%mpi_grp, yiak)
549 end if
550
552
554
555 end do ! ifreq
556
557 safe_deallocate_a(et)
558 safe_deallocate_a(xiak)
559 safe_deallocate_a(yiak)
560 safe_deallocate_a(tdm)
561 safe_deallocate_a(tdm_1d)
562
563 safe_deallocate_a(psi)
564 safe_deallocate_a(upsi)
565 if(sys%kpoints%use_symmetries) then
566 safe_deallocate_p(psi_sym)
567 safe_deallocate_p(upsi_sym)
568 end if
569 safe_deallocate_a(ftcmplx)
570 safe_deallocate_a(centers)
571 safe_deallocate_a(phase)
572 safe_deallocate_a(omega)
573
574 safe_deallocate_p(sys)
575 call states_elec_end(gs_st)
576 call fft_all_end()
577 call io_end()
579 call messages_end()
580 call parser_end()
581 call global_end()
582
583contains
584
585 ! -----------------------------------------------------------------
586 ! Determines the position of the hole, either from the input or using the
587 ! first atom in the cell.
588 ! This returns the index of the point in the mesh closest to the position.
589 subroutine tdtdm_get_hole_position(xx_h, ip_h)
590 real(real64), intent(out) :: xx_h(1:sys%space%dim)
591 integer, intent(out) :: ip_h
592
593 real(real64) :: dmin
594 integer :: idir, rankmin
595
597
598 !%Variable TDTDMHoleCoordinates
599 !%Type float
600 !%Section Utilities::oct-tdtdm
601 !%Description
602 !% The position of the hole used to compute the TDTDM,
603 !% in Cartesian coordinates.
604 !% Note that the code will use the closest grid point.
605 !%
606 !% The coordinates of the hole are specified in the following way
607 !% <tt>%TDTDMHoleCoordinates
608 !% <br>&nbsp;&nbsp;hole_x | hole_y | hole_z
609 !% <br>%</tt>
610 !%
611 !% If TDTDMHoleCoordinates or TDTDMHoleReducedCoordinates are not specified,
612 !% the code will use the coordinate of the first atom in the cell.
613 !%End
614
615 if(parse_block(global_namespace, 'TDTDMHoleCoordinates', blk) == 0) then
616 if(parse_block_cols(blk,0) < sys%space%dim) then
617 call messages_input_error(global_namespace, 'TDTDMHoleCoordinates')
618 end if
619 do idir = 1, sys%space%dim
620 call parse_block_float(blk, 0, idir - 1, xx_h(idir), units_inp%length)
621 end do
622 call parse_block_end(blk)
623 else
624 !%Variable TDTDMHoleReducedCoordinates
625 !%Type float
626 !%Section Utilities::oct-tdtdm
627 !%Description
628 !% Same as TDTDMHoleCoordinates, except that coordinates are given in reduced coordinates
629 !%End
630
631 if(parse_block(global_namespace, 'TDTDMHoleReducedCoordinates', blk) == 0) then
632 if(parse_block_cols(blk,0) < sys%space%dim) then
633 call messages_input_error(global_namespace, 'TDTDMHoleReducedCoordinates')
634 end if
635 do idir = 1, sys%space%dim
636 call parse_block_float(blk, 0, idir - 1, xx_h(idir), units_inp%length)
637 end do
638 call parse_block_end(blk)
639 xx_h = sys%ions%latt%red_to_cart(xx_h)
640 else
641 xx_h(1:sys%space%dim) = sys%ions%pos(1:sys%space%dim, 1)
642 end if
643 end if
644
645 ! We bring back the hole into the cell
646 xx_h = sys%ions%latt%fold_into_cell(xx_h)
647
648 ! At the moment, we ignore rankmin
649 assert(.not.sys%gr%parallel_in_domains)
650 ip_h = mesh_nearest_point(sys%gr, xx_h, dmin, rankmin)
651 write(message(1), '(a, 3(1x,f7.4,a))') "oct-tdtdm: Requesting the hole at (", xx_h(1), &
652 ",", xx_h(2), ",", xx_h(3), ")."
653 call mesh_r(sys%gr, ip_h, dmin, coords=xx_h)
654 write(message(2), '(a, 3(1x,f7.4,a))') "oct-tdtdm: Setting the hole at (", xx_h(1), &
655 ",", xx_h(2), ",", xx_h(3), ")."
656
657 call messages_info(2)
658
660 end subroutine tdtdm_get_hole_position
661
662 subroutine tdtdm_output_density()
663 real(real64), allocatable :: den(:,:), den_1d(:,:,:,:)
664 real(real64) :: norm, xx(3), xx_h(3)
665 integer :: iunit
666
667 push_sub(tdtdm_output_density)
668
669 ! We compute the TDM density
670 select case(sys%space%dim)
671 case(2,3)
672 safe_allocate(den(1:sys%gr%np, 1:nreplica))
673 do irep = 1, nreplica
674 do ii = 1, sys%gr%np
675 den(ii, irep) = real(tdm(ii, irep)*conjg(tdm(ii, irep)), real64)
676 end do
677 end do
678
679 ! Here we renormalize to avoid too small numbers in the outputs
680 norm = maxval(den)
681 call lalg_scal(sys%gr%np, nreplica, m_one/norm, den)
683 case(1)
684 safe_allocate(den_1d(1:sys%gr%np, 1:sys%gr%np, 1:nreplica, 1:nreplica))
685 do irep_h = 1, nreplica
686 do irep = 1, nreplica
687 do ip_h = 1, sys%gr%np
688 do ii = 1, sys%gr%np
689 tdm_1d(ii, ip_h, irep, irep_h) = conjg(tdm_1d(ii, ip_h, irep, irep_h))
690 den_1d(ii, ip_h, irep, irep_h) = real(tdm_1d(ii, ip_h, irep, irep_h)*conjg(tdm_1d(ii,ip_h, irep, irep_h)), real64)
691 end do
692 end do
693 end do
694 end do
695 end select
696
697 fn_unit = units_out%length**(-sys%space%dim)
698
699 select case(sys%space%dim)
700 case(2,3)
701 write(fname, '(a, f0.4)') 'tdm_density-0', omega(ifreq)
702 call io_function_output_supercell(io_function_fill_how("XCrySDen"), "td.general", fname, &
703 sys%gr, sys%space, sys%ions%latt, den, centers, supercell, fn_unit, &
704 ierr, global_namespace, pos=sys%ions%pos, atoms=sys%ions%atom, grp = st%dom_st_kpt_mpi_grp, extra_atom=pos_h)
705
706 call io_function_output_supercell(io_function_fill_how("PlaneZ"), "td.general", fname, &
707 sys%gr, sys%space, sys%ions%latt, den, centers, supercell, fn_unit, &
708 ierr, global_namespace, grp = st%dom_st_kpt_mpi_grp)
709
710 safe_deallocate_a(den)
711
712 case(1)
713
714 call tdtdm_get_hole_position(pos_h, ip_h)
715 irep_h = floor(supercell(1)/m_two)
716
717 write(fname, '(a, f0.4)') 'tdm_density-0', omega(ifreq)
718 call io_function_output_supercell(io_function_fill_how("AxisX"), "td.general", fname, &
719 sys%gr, sys%space, sys%ions%latt, &
720 den_1d(:,ip_h,:,irep_h), centers, supercell, fn_unit, ierr, global_namespace, &
721 grp = st%dom_st_kpt_mpi_grp)
722
723 write(fname, '(a, f0.4)') 'tdm_wfn-0', omega(ifreq)
724 call io_function_output_supercell(io_function_fill_how("AxisX"), "td.general", fname, &
725 sys%gr, sys%space, sys%ions%latt, &
726 tdm_1d(:,ip_h,:,irep_h), centers, supercell, fn_unit, ierr, global_namespace, &
727 grp = st%dom_st_kpt_mpi_grp)
728
729 assert(.not.sys%gr%parallel_in_domains)
730 if (mpi_grp_is_root(mpi_world)) then
731 write(fname, '(a, f0.4)') 'td.general/tdm_density-0', omega(ifreq)
732 iunit = io_open(fname, action='write')
733 write(iunit, '(a)', iostat=ierr) '# r_e r_h Re(\Psi(r_e,r_h)) Im(\Psi(r_e,r_h)) |\Psi(r_e,r_h)|^2'
734
735 do irep_h = 1, nreplica
736 do ip_h = 1, sys%gr%np
737 xx_h = units_from_atomic(units_out%length, mesh_x_global(sys%gr, i4_to_i8(ip_h)) &
738 + centers(1:sys%space%dim, irep_h))
739
740 do irep = 1, nreplica
741 do ii = 1, sys%gr%np
742 xx = units_from_atomic(units_out%length, mesh_x_global(sys%gr, i4_to_i8(ii)) &
743 + centers(1:sys%space%dim, irep))
744 write(iunit, '(5es23.14E3)', iostat=ierr) xx(1), xx_h(1), &
745 real(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h)), real64) ,&
746 aimag(units_from_atomic(fn_unit, tdm_1D(ii, ip_h, irep, irep_h))), &
747 units_from_atomic(fn_unit, den_1D(ii, ip_h, irep, irep_h))
748 end do
749 end do
750 end do
751 end do
752 end if
753
754 safe_deallocate_a(den_1d)
755 end select
756
757
758 pop_sub(tdtdm_output_density)
759 end subroutine tdtdm_output_density
760
761 subroutine tdtdm_excitonic_weight()
762 real(real64), allocatable :: weight(:,:)
763
764 if (.not. mpi_grp_is_root(mpi_world)) return
765
766 push_sub(tdtdm_excitonic_weight)
767
768 safe_allocate(weight(1:st%nik, 1:gs_st%nst))
769 weight = m_zero
770
771 do ik = 1, st%nik
772 do ist = 1, st%nst
773 if(abs(gs_st%occ(ist, ik)) < m_epsilon) cycle
774
775 do uist = ist+1, gs_st%nst
776 if(abs(gs_st%occ(uist, ik)) > m_epsilon) cycle
777
778 weight(ik, ist) = weight(ik, ist) + abs(xiak(ist, uist, ik))**2
779 weight(ik, uist) = weight(ik,uist) + abs(yiak(ist, uist, ik))**2
780 end do
781 end do
782 end do
783
784 write(fname, '(a, f0.4)') 'td.general/tdm_weights-0', omega(ifreq)
785 out_file = io_open(fname, action='write')
786 write(out_file, '(a)') '# ik - kx - ky - kz - sum weights - eigenval and weights(ist,ik) '
787 do ik = 1, st%nik
788 ikpoint = st%d%get_kpoint_index(ik)
789 kpoint(1:sys%space%dim) = sys%kpoints%reduced%point1BZ(1:sys%space%dim,ikpoint)
790 write(out_file, '(i4,4e15.6)', advance='no') ik, kpoint(1:3), sum(weight(ik, 1:gs_st%nst))
791 do uist = 1, gs_st%nst-1
792 write(out_file, '(2e15.6)', advance='no') gs_st%eigenval(uist, ik), weight(ik, uist)
793 end do
794 write(out_file, '(e15.6)') weight(ik, uist)
795 end do
796 call io_close(out_file)
797
798 safe_deallocate_a(weight)
799
801 end subroutine tdtdm_excitonic_weight
802
803end program tdtdm
804
805!! Local Variables:
806!! mode: f90
807!! coding: utf-8
808!! End:
initialize a batch with existing memory
Definition: batch.F90:267
constant times a vector plus a vector
Definition: lalg_basic.F90:171
scales a vector by a constant
Definition: lalg_basic.F90:157
double hypot(double __x, double __y) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double floor(double __x) __attribute__((__nothrow__
This module implements batches of mesh functions.
Definition: batch.F90:133
This module handles the calculation mode.
type(calc_mode_par_t), public calc_mode_par
Singleton instance of parallel calculation mode.
integer, parameter, public p_strategy_states
parallelization in states
integer, parameter, public spinors
Fast Fourier Transform module. This module provides a single interface that works with different FFT ...
Definition: fft.F90:118
subroutine, public fft_all_init(namespace)
initialize the table
Definition: fft.F90:278
subroutine, public fft_all_end()
delete all plans
Definition: fft.F90:391
real(real64), parameter, public m_two
Definition: global.F90:190
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
Definition: global.F90:382
real(real64), parameter, public m_zero
Definition: global.F90:188
complex(real64), parameter, public m_z0
Definition: global.F90:198
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_epsilon
Definition: global.F90:204
subroutine, public global_init(communicator)
Initialise Octopus.
Definition: global.F90:325
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public zgrid_symmetrize_single(gr, iop, field, symm_field, suppress_warning)
Definition: grid.F90:859
integer(int64) function, public io_function_fill_how(where)
Use this function to quickly plot functions for debugging purposes: call dio_function_output(io_funct...
Definition: io.F90:114
subroutine, public io_init(defaults)
If the argument defaults is present and set to true, then the routine will not try to read anything f...
Definition: io.F90:161
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
subroutine, public io_skip_header(iunit)
Definition: io.F90:597
subroutine, public io_end()
Definition: io.F90:258
character(len=max_path_len) function, public io_workpath(path, namespace)
Definition: io.F90:270
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
integer pure function, public kpoints_get_num_symmetry_ops(this, ik)
Definition: kpoints.F90:1559
integer pure function, public kpoints_get_symmetry_ops(this, ik, index)
Definition: kpoints.F90:1572
subroutine, public kpoints_to_reduced(latt, kin, kout)
Definition: kpoints.F90:1044
subroutine, public kpoints_to_absolute(latt, kin, kout)
Definition: kpoints.F90:1031
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
integer function, public mesh_nearest_point(mesh, pos, dmin, rankmin)
Returns the index of the point which is nearest to a given vector position pos.
Definition: mesh.F90:380
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:336
real(real64) function, dimension(1:mesh%box%dim), public mesh_x_global(mesh, ipg)
Definition: mesh.F90:804
subroutine, public messages_end()
Definition: messages.F90:277
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_init(output_dir)
Definition: messages.F90:224
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
type(namespace_t), public global_namespace
Definition: namespace.F90:132
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
subroutine, public parser_init()
Initialise the Octopus parser.
Definition: parser.F90:450
subroutine, public parser_end()
End the Octopus parser.
Definition: parser.F90:481
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public profiling_end(namespace)
Definition: profiling.F90:413
subroutine, public profiling_init(namespace)
Create profiling subdirectory.
Definition: profiling.F90:255
subroutine, public restart_module_init(namespace)
Definition: restart.F90:308
subroutine, public restart_init(restart, namespace, data_type, type, mc, ierr, mesh, dir, exact)
Initializes a restart object.
Definition: restart.F90:516
integer, parameter, public restart_proj
Definition: restart.F90:200
integer, parameter, public restart_type_load
Definition: restart.F90:245
subroutine, public restart_end(restart)
Definition: restart.F90:722
subroutine, public spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
Definition: spectrum.F90:2507
subroutine, public spectrum_fourier_transform(method, transform, noise, time_start, time_end, t0, time_step, time_function, energy_start, energy_end, energy_step, energy_function)
Computes the sine, cosine, (or "exponential") Fourier transform of the real function given in the tim...
Definition: spectrum.F90:2637
subroutine, public spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
Definition: spectrum.F90:213
integer default
Definition: spectrum.F90:207
integer, parameter, public spectrum_transform_cos
Definition: spectrum.F90:171
integer, parameter, public spectrum_transform_sin
Definition: spectrum.F90:171
subroutine, public spectrum_count_time_steps(namespace, iunit, time_steps, dt)
Definition: spectrum.F90:2385
pure integer function, public spectrum_nenergy_steps(spectrum)
Definition: spectrum.F90:2946
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_distribute_nodes(st, namespace, mc)
@Brief. Distribute states over the processes for states parallelization
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public kpoints_distribute(this, mc)
distribute k-points over the nodes in the corresponding communicator
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_look(restart, nik, dim, nst, ierr)
Reads the 'states' file in the restart directory, and finds out the nik, dim, and nst contained in it...
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_load(restart, namespace, space, st, mesh, kpoints, ierr, iter, lr, lowest_missing, label, verbose, skip)
returns in ierr: <0 => Fatal error, or nothing read =0 => read all wavefunctions >0 => could only rea...
subroutine, public symmetries_apply_kpoint_red(this, iop, aa, bb)
Definition: symmetries.F90:548
type(type_t), public type_cmplx
Definition: types.F90:134
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
subroutine, public unit_system_init(namespace)
type(unit_system_t), public units_inp
the units systems for reading and writing
Class describing the electron system.
Definition: electrons.F90:218
int true(void)
subroutine tdtdm_excitonic_weight()
Definition: tdtdm.F90:855
subroutine tdtdm_output_density()
Definition: tdtdm.F90:756
subroutine tdtdm_get_hole_position(xx_h, ip_h)
Definition: tdtdm.F90:683
program tdtdm
Definition: tdtdm.F90:114