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