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
21program dielectric_function
22  use batch_oct_m
23  use command_line_oct_m
24  use geometry_oct_m
25  use global_oct_m
26  use io_oct_m
27  use lalg_adv_oct_m
28  use messages_oct_m
29  use namespace_oct_m
30  use parser_oct_m
31  use profiling_oct_m
32  use space_oct_m
33  use spectrum_oct_m
34  use simul_box_oct_m
35  use unit_oct_m
36  use unit_system_oct_m
37
38  implicit none
39
40  integer :: in_file, out_file, ref_file, ii, jj, kk, idir, jdir, ierr
41  integer :: time_steps, time_steps_ref, energy_steps, istart, iend, ntiter
42  FLOAT   :: dt, dt_ref, tt, ww, n0
43  FLOAT, allocatable :: vecpot(:, :), vecpot0(:), ftreal(:, :), ftimag(:, :)
44  CMPLX, allocatable :: dielectric(:, :), chi(:, :), invdielectric(:, :), fullmat(:, :)
45  FLOAT, allocatable :: vecpot_ref(:, :)
46  type(spectrum_t)  :: spectrum
47  type(block_t)     :: blk
48  type(space_t)     :: space
49  type(geometry_t)  :: geo
50  type(simul_box_t) :: sb
51  type(batch_t)     :: vecpotb, ftrealb, ftimagb
52  character(len=120) :: header
53  FLOAT :: start_time
54  character(len=MAX_PATH_LEN) :: ref_filename
55
56  ! Initialize stuff
57  call global_init(is_serial = .true.)
58
59  call getopt_init(ierr)
60  if(ierr == 0) call getopt_dielectric_function()
61  call getopt_end()
62
63  call parser_init()
64
65  call messages_init()
66
67  call io_init()
68
69  call unit_system_init(global_namespace)
70
71  call spectrum_init(spectrum, global_namespace)
72
73  call space_init(space, global_namespace)
74  call geometry_init(geo, global_namespace, space)
75  call simul_box_init(sb, global_namespace, geo, space)
76
77  SAFE_ALLOCATE(vecpot0(1:space%dim))
78
79  if(parse_block(global_namespace, 'GaugeVectorField', blk) == 0) then
80
81    do ii = 1, space%dim
82      call parse_block_float(blk, 0, ii - 1, vecpot0(ii))
83    end do
84
85    call parse_block_end(blk)
86
87  else
88
89    message(1) = "Cannot find the GaugeVectorField in the input file"
90    call messages_fatal(1)
91
92  end if
93
94  message(1) = "This program assumes that the gauge field is in the 'x'"
95  message(2) = "direction, and that the 'y' and 'z' directions are equivalent."
96  message(3) = "If this is not the case the dielectric function and the"
97  message(4) = "susceptibility will be wrong."
98  call messages_warning(4)
99
100  start_time = spectrum%start_time
101  call parse_variable(global_namespace, 'GaugeFieldDelay', start_time, spectrum%start_time )
102
103  in_file = io_open('td.general/gauge_field', global_namespace, action='read', status='old', die=.false.)
104  if(in_file < 0) then
105    message(1) = "Cannot open file '"//trim(io_workpath('td.general/gauge_field', global_namespace))//"'"
106    call messages_fatal(1)
107  end if
108  call io_skip_header(in_file)
109  call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
110
111  if(parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
112    !%Variable TransientAbsorptionReference
113    !%Type string
114    !%Default "."
115    !%Section Utilities::oct-propagation_spectrum
116    !%Description
117    !% In case of delayed kick, the calculation of the transient absorption requires
118    !% to substract a reference calculation, containing the gauge-field without the kick
119    !% This reference must be computed using GaugeFieldPropagate=yes and to have
120    !% TDOutput = gauge_field.
121    !% This variables defined the directory in which the reference gauge_field field is,
122    !% relative to the current folder
123    !%End
124
125    call parse_variable(global_namespace, 'TransientAbsorptionReference', '.', ref_filename)
126    ref_file = io_open(trim(ref_filename)//'/gauge_field', global_namespace, action='read', status='old', die=.false.)
127    if(ref_file < 0) then
128      message(1) = "Cannot open reference file '"// &
129        trim(io_workpath(trim(ref_filename)//'/gauge_field', global_namespace))//"'"
130      call messages_fatal(1)
131    end if
132    call io_skip_header(ref_file)
133    call spectrum_count_time_steps(global_namespace, ref_file, time_steps_ref, dt_ref)
134    if(time_steps_ref < time_steps) then
135      message(1) = "The reference calculation does not contain enought time steps"
136      call messages_fatal(1)
137    end if
138
139    if(dt_ref /= dt) then
140      message(1) = "The time step of the reference calculation is different from the current calculation"
141      call messages_fatal(1)
142    end if
143
144  end if
145
146  time_steps = time_steps + 1
147
148  SAFE_ALLOCATE(vecpot(1:time_steps, space%dim*3))
149
150  call io_skip_header(in_file)
151
152  do ii = 1, time_steps
153    read(in_file, *) jj, tt, (vecpot(ii, kk), kk = 1, space%dim*3)
154  end do
155
156  call io_close(in_file)
157
158  !We remove the reference
159  if(parse_is_defined(global_namespace, 'TransientAbsorptionReference')) then
160    time_steps_ref = time_steps_ref + 1
161    SAFE_ALLOCATE(vecpot_ref(1:time_steps_ref, space%dim*3))
162    call io_skip_header(ref_file)
163    do ii = 1, time_steps_ref
164      read(ref_file, *) jj, tt, (vecpot_ref(ii, kk), kk = 1, space%dim*3)
165    end do
166    call io_close(ref_file)
167    do ii = 1, time_steps
168      do kk = 1, space%dim*3
169        vecpot(ii, kk) = vecpot(ii, kk) - vecpot_ref(ii, kk)
170      end do
171    end do
172  end if
173
174  write(message(1), '(a, i7, a)') "Info: Read ", time_steps, " steps from file '"// &
175    trim(io_workpath('td.general/gauge_field', global_namespace))//"'"
176  call messages_info(1)
177
178
179  ! Find out the iteration numbers corresponding to the time limits.
180  call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
181
182  istart = max(1, istart)
183
184  energy_steps = spectrum_nenergy_steps(spectrum)
185
186  n0 = sqrt(sum(vecpot0(1:space%dim))**2)
187
188  SAFE_ALLOCATE(ftreal(1:energy_steps, 1:space%dim))
189  SAFE_ALLOCATE(ftimag(1:energy_steps, 1:space%dim))
190
191  call batch_init(vecpotb, 1, 1, space%dim, vecpot(:, space%dim+1:space%dim*2))
192  call batch_init(ftrealb, 1, 1, space%dim, ftreal)
193  call batch_init(ftimagb, 1, 1, space%dim, ftimag)
194
195  call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, dt, vecpotb)
196
197  call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_COS, spectrum%noise, &
198    istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
199    spectrum%max_energy, spectrum%energy_step, ftrealb)
200
201  call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_SIN, spectrum%noise, &
202    istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
203    spectrum%max_energy, spectrum%energy_step, ftimagb)
204
205
206  call vecpotb%end()
207  call ftrealb%end()
208  call ftimagb%end()
209
210  SAFE_ALLOCATE(invdielectric(1:space%dim, 1:energy_steps))
211  SAFE_ALLOCATE(dielectric(1:space%dim, 1:energy_steps))
212  SAFE_ALLOCATE(chi(1:space%dim, 1:energy_steps))
213  SAFE_ALLOCATE(fullmat(1:space%dim, 1:space%dim))
214
215  do kk = 1, energy_steps
216    ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
217
218    invdielectric(1:space%dim, kk) = (vecpot0(1:space%dim) + TOCMPLX(ftreal(kk, 1:space%dim), ftimag(kk, 1:space%dim)))/n0
219
220    ! calculate the full inverse dielectric matrix
221    do idir = 1, space%dim
222      do jdir = 1, space%dim
223        fullmat(idir, mod(jdir + idir - 2, space%dim) + 1) = invdielectric(jdir, kk)
224      end do
225    end do
226
227    ! and invert it
228    call lalg_sym_inverter('u', space%dim, fullmat)
229
230    ! symmetrize the rest dielectric matrix
231    do idir = 2, space%dim
232      do jdir = 1, idir - 1
233        fullmat(jdir, idir) = fullmat(idir, jdir)
234      end do
235    end do
236
237    dielectric(1:space%dim, kk) = fullmat(1:space%dim, 1)
238
239    chi(1:space%dim, kk) = (dielectric(1:space%dim, kk) - vecpot0(1:space%dim)/n0)*sb%rcell_volume/(M_FOUR*M_PI)
240  end do
241
242  SAFE_DEALLOCATE_A(fullmat)
243
244  out_file = io_open('td.general/inverse_dielectric_function', global_namespace, action='write')
245  select case(space%dim)
246  case(1)
247    write(header, '(7a15)') '#        energy', 'Re x', 'Im x'
248  case(2)
249    write(header, '(7a15)') '#        energy', 'Re x', 'Im x', 'Re y', 'Im y'
250  case(3)
251    write(header, '(7a15)') '#        energy', 'Re x', 'Im x', 'Re y', 'Im y', 'Re z', 'Im z'
252  end select
253
254
255  write(out_file,'(a)') trim(header)
256  do kk = 1, energy_steps
257    ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
258    write(out_file, '(e15.6)', advance='no') ww
259    do idir = 1, space%dim
260      write(out_file, '(2e15.6)', advance='no') TOFLOAT(invdielectric(idir, kk)), aimag(invdielectric(idir, kk))
261    end do
262    write(out_file, '()')
263  end do
264  call io_close(out_file)
265
266  out_file = io_open('td.general/dielectric_function', global_namespace, action='write')
267  write(out_file,'(a)') trim(header)
268  do kk = 1, energy_steps
269    ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
270    write(out_file, '(e15.6)', advance='no') ww
271    do idir = 1, space%dim
272      write(out_file, '(2e15.6)', advance='no') TOFLOAT(dielectric(idir, kk)), aimag(dielectric(idir, kk))
273    end do
274    write(out_file, '()')
275  end do
276  call io_close(out_file)
277
278  out_file = io_open('td.general/chi', global_namespace, action='write')
279  write(out_file,'(a)') trim(header)
280  do kk = 1, energy_steps
281    dielectric(1:3, kk) = (dielectric(1:3, kk) - M_ONE)/(CNST(4.0)*M_PI)
282    ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
283    write(out_file, '(e15.6)', advance='no') ww
284    do idir = 1, space%dim
285      write(out_file, '(2e15.6)', advance='no') TOFLOAT(chi(idir, kk)), aimag(chi(idir, kk))
286    end do
287    write(out_file, '()')
288  end do
289  call io_close(out_file)
290
291  SAFE_DEALLOCATE_A(dielectric)
292  SAFE_DEALLOCATE_A(invdielectric)
293  SAFE_DEALLOCATE_A(chi)
294  SAFE_DEALLOCATE_A(vecpot)
295  SAFE_DEALLOCATE_A(vecpot_ref)
296  SAFE_DEALLOCATE_A(vecpot0)
297  SAFE_DEALLOCATE_A(ftreal)
298  SAFE_DEALLOCATE_A(ftimag)
299
300  call simul_box_end(sb)
301  call geometry_end(geo)
302  call space_end(space)
303  call io_end()
304  call messages_end()
305
306  call parser_end()
307  call global_end()
308
309end program dielectric_function
310
311!! Local Variables:
312!! mode: f90
313!! coding: utf-8
314!! End:
315