Octopus
dielectric_function.F90
Go to the documentation of this file.
1!! Copyright (C) 2008 X. Andrade
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
30 use batch_oct_m
32 use global_oct_m
33 use io_oct_m
36 use math_oct_m
39 use parser_oct_m
41 use space_oct_m
43 use unit_oct_m
45
46 implicit none
47
48 integer :: ierr
49
50 ! Initialize stuff
52
53 call getopt_init(ierr)
54 if (ierr == 0) call getopt_dielectric_function()
55 call getopt_end()
56
57 call parser_init()
58
59 call messages_init()
60
61 call io_init()
63
65
67
69 call io_end()
70 call messages_end()
71
72 call parser_end()
73 call global_end()
74
75contains
76
78
79 integer :: in_file, out_file, ref_file, ii, jj, kk
80 integer :: time_steps, time_steps_ref, energy_steps, istart, iend, ntiter
81 real(real64) :: dt, dt_ref, tt, ww, norm_Aext
82 real(real64), allocatable :: vecpot(:, :), Aext(:), Eind_real(:, :), Eind_imag(:, :)
83 complex(real64), allocatable :: dielectric(:), chi(:), invdielectric(:)
84 real(real64), allocatable :: vecpot_ref(:, :)
85 type(spectrum_t) :: spectrum
86 type(block_t) :: blk
87 type(space_t) :: space
88 type(lattice_vectors_t) :: latt
89 type(batch_t) :: vecpotb, Eind_realb, Eind_imagb
90 character(len=120) :: header
91 real(real64) :: start_time
92 character(len=MAX_PATH_LEN) :: ref_filename
93 complex(real64),allocatable :: Eind(:)
94
95 call spectrum_init(spectrum, global_namespace)
96
99
100 safe_allocate(aext(1:space%dim))
101 safe_allocate(eind(1:space%dim))
102
103 if (parse_block(global_namespace, 'GaugeVectorField', blk) == 0) then
104
105 do ii = 1, space%dim
106 call parse_block_float(blk, 0, ii - 1, aext(ii))
107 end do
108
109 call parse_block_end(blk)
110
111 else
112
113 message(1) = "Cannot find the GaugeVectorField in the input file"
114 call messages_fatal(1)
115
116 end if
117
118 if(space%dim > 1) then
119 message(1) = "This program assumes that the gauge field is along a "
120 message(2) = "direction for which the dielectric tensor has no off-diagonal terms."
121 message(3) = "If this is not the case the dielectric function and the"
122 message(4) = "susceptibility will be wrong."
123 call messages_warning(4)
124 end if
125
126 start_time = spectrum%start_time
127 call parse_variable(global_namespace, 'GaugeFieldDelay', start_time, spectrum%start_time)
128
129 in_file = io_open('td.general/gauge_field', global_namespace, action='read', status='old')
130
131 ! Get the number of iterations and the time-step
132 call io_skip_header(in_file)
133 call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
134
135 ! Fix the correct time at which the kick was applied.
136 ! This guaranties that the first time has exactly a zero vector potential
137 ! If we do not do this, we can get some artifacts for transient absorption.
138 spectrum%start_time = ceiling(spectrum%start_time/dt)*dt
139
140 if (parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
141 !%Variable TransientAbsorptionReference
142 !%Type string
143 !%Default "."
144 !%Section Utilities::oct-propagation_spectrum
145 !%Description
146 !% In case of delayed kick, the calculation of the transient absorption requires
147 !% to substract a reference calculation, containing the gauge-field without the kick
148 !% This reference must be computed using GaugeFieldPropagate=yes and to have
149 !% TDOutput = gauge_field.
150 !% This variables defined the directory in which the reference gauge_field field is,
151 !% relative to the current folder
152 !%End
153
154 call parse_variable(global_namespace, 'TransientAbsorptionReference', '.', ref_filename)
155 ref_file = io_open(trim(ref_filename)//'/gauge_field', global_namespace, action='read', status='old')
156 call io_skip_header(ref_file)
157 call spectrum_count_time_steps(global_namespace, ref_file, time_steps_ref, dt_ref)
158 if (time_steps_ref < time_steps) then
159 message(1) = "The reference calculation does not contain enought time steps"
160 call messages_fatal(1)
161 end if
162
163 if (.not. is_close(dt_ref, dt)) then
164 message(1) = "The time step of the reference calculation is different from the current calculation"
165 call messages_fatal(1)
166 end if
167
168 end if
169
170 ! Add one as zero is included in the time signal
171 time_steps = time_steps + 1
172
173 safe_allocate(vecpot(1:time_steps, 1:space%dim*3))
174
175 call io_skip_header(in_file)
176
177 do ii = 1, time_steps
178 read(in_file, *) jj, tt, (vecpot(ii, kk), kk = 1, space%dim*3)
179 end do
180
181 call io_close(in_file)
182
183 !We remove the reference
184 if (parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
185 time_steps_ref = time_steps_ref + 1
186 safe_allocate(vecpot_ref(1:time_steps_ref, 1:space%dim*3))
187 call io_skip_header(ref_file)
188 do ii = 1, time_steps_ref
189 read(ref_file, *) jj, tt, (vecpot_ref(ii, kk), kk = 1, space%dim*3)
190 end do
191 call io_close(ref_file)
192 do ii = 1, time_steps
193 do kk = 1, space%dim*3
194 vecpot(ii, kk) = vecpot(ii, kk) - vecpot_ref(ii, kk)
195 end do
196 end do
197 end if
198
199 write(message(1), '(a, i7, a)') "Info: Read ", time_steps, " steps from file '"// &
200 trim(io_workpath('td.general/gauge_field', global_namespace))//"'"
201 call messages_info(1)
202
203
204 ! Find out the iteration numbers corresponding to the time limits.
205 ! Max time correspond to time_steps-1, as we start from 0.
206 call spectrum_fix_time_limits(spectrum, time_steps-1, dt, istart, iend, ntiter)
207
208 ! We need to fix istart because the damping starts are (istart-1)*dt, which needs to be spectrum%start_time
209 istart = istart + 1
210
211 energy_steps = spectrum_nenergy_steps(spectrum)
212
213 norm_aext = norm2(aext(1:space%dim))
214
215 safe_allocate(eind_real(1:energy_steps, 1:space%dim))
216 safe_allocate(eind_imag(1:energy_steps, 1:space%dim))
217
218 ! We select dA/dt = E
219 call batch_init(vecpotb, 1, 1, space%dim, vecpot(:, space%dim+1:space%dim*2))
220 call batch_init(eind_realb, 1, 1, space%dim, eind_real)
221 call batch_init(eind_imagb, 1, 1, space%dim, eind_imag)
222
223 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, dt, vecpotb)
224
225 call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
226 istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
227 spectrum%max_energy, spectrum%energy_step, eind_realb)
228
229 call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
230 istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
231 spectrum%max_energy, spectrum%energy_step, eind_imagb)
232
233
234 call vecpotb%end()
235 call eind_realb%end()
236 call eind_imagb%end()
237
238 safe_allocate(invdielectric(1:energy_steps))
239 safe_allocate(dielectric(1:energy_steps))
240 safe_allocate(chi(1:energy_steps))
241
242 do kk = 1, energy_steps
243 ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
244
245 ! We compute 1/\epsilon(\omega), see Eq. (4) in Bertsch et al. PRB 62, 7998
246 ! More precisely, we have \epsilon^{-1}_{\alpha\beta} = \frac{E_{tot,\alpha}}{E_{ext,\beta}}
247
248 ! Here we can only assume that the Gauge field is along an axis without off-diagonal terms
249 ! Else, we need several calulations
250 eind = cmplx(eind_real(kk, 1:space%dim), eind_imag(kk, 1:space%dim), real64)
251 invdielectric(kk) = dot_product(aext, aext + eind) /norm_aext**2
252
253 dielectric(kk) = m_one / invdielectric(kk)
254
255 chi(kk) = (dielectric(kk) - m_one)*latt%rcell_volume/(m_four*m_pi)
256 end do
257
258 out_file = io_open('td.general/inverse_dielectric_function', global_namespace, action='write')
259 write(header, '(7a15)') '# energy', 'Re', 'Im'
260
261 write(out_file,'(a)') trim(header)
262 do kk = 1, energy_steps
263 ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
264 write(out_file, '(e15.6)', advance='no') ww
265 write(out_file, '(2e15.6)', advance='no') real(invdielectric(kk), real64), aimag(invdielectric(kk))
266 write(out_file, '()')
267 end do
268 call io_close(out_file)
269
270 out_file = io_open('td.general/dielectric_function', global_namespace, action='write')
271 write(out_file,'(a)') trim(header)
272 do kk = 1, energy_steps
273 ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
274 write(out_file, '(e15.6)', advance='no') ww
275 write(out_file, '(2e15.6)', advance='no') real(dielectric(kk), real64), aimag(dielectric(kk))
276 write(out_file, '()')
277 end do
278 call io_close(out_file)
279
280 out_file = io_open('td.general/chi', global_namespace, action='write')
281 write(out_file,'(a)') trim(header)
282 do kk = 1, energy_steps
283 ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
284 write(out_file, '(e15.6)', advance='no') ww
285 write(out_file, '(2e15.6)', advance='no') real(chi(kk), real64), aimag(chi(kk))
286 write(out_file, '()')
287 end do
288 call io_close(out_file)
289
290 safe_deallocate_a(dielectric)
291 safe_deallocate_a(invdielectric)
292 safe_deallocate_a(chi)
293 safe_deallocate_a(vecpot)
294 safe_deallocate_a(vecpot_ref)
295 safe_deallocate_a(aext)
296 safe_deallocate_a(eind)
297 safe_deallocate_a(eind_real)
298 safe_deallocate_a(eind_imag)
299
300 end subroutine dielectric_function_compute
301
302end program dielectric_function
303
304!! Local Variables:
305!! mode: f90
306!! coding: utf-8
307!! End:
subroutine dielectric_function_compute()
program dielectric_function
A utility used to obtain the dielectric function from a kick calculation using the gauge field approa...
initialize a batch with existing memory
Definition: batch.F90:267
This module implements batches of mesh functions.
Definition: batch.F90:133
subroutine, public getopt_init(ierr)
Initializes the getopt machinery. Must be called before attempting to parse the options....
subroutine, public getopt_end
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
Definition: global.F90:382
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
type(mpi_comm), parameter, public serial_dummy_comm
Alias MPI_COMM_UNDEFINED for the specific use case of initialising Octopus utilities with no MPI supp...
Definition: global.F90:265
subroutine, public init_octopus_globals(comm)
Initialise Octopus-specific global constants and files. This routine performs no initialisation calls...
Definition: global.F90:354
real(real64), parameter, public m_one
Definition: global.F90:189
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
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public messages_end()
Definition: messages.F90:277
subroutine, public messages_init(output_dir)
Definition: messages.F90:224
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
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_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
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 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
subroutine, public spectrum_signal_damp(damp_type, damp_factor, time_start, time_end, t0, time_step, time_function)
Definition: spectrum.F90:2551
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
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.
subroutine, public unit_system_init(namespace)