1!! Copyright (C) 2015 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 21 program conductivity 22 use batch_oct_m 23 use command_line_oct_m 24 use geometry_oct_m 25 use global_oct_m 26 use grid_oct_m 27 use io_oct_m 28 use messages_oct_m 29 use namespace_oct_m 30 use parser_oct_m 31 use profiling_oct_m 32 use simul_box_oct_m 33 use space_oct_m 34 use spectrum_oct_m 35 use species_oct_m 36 use states_elec_oct_m 37 use unit_oct_m 38 use unit_system_oct_m 39 40 implicit none 41 42 integer :: iunit, ierr, ii, jj, iter, read_iter, ntime, nvel, ivel 43 integer :: istart, iend, energy_steps, out_file 44 FLOAT, allocatable :: time(:), velocities(:, :) 45 FLOAT, allocatable :: total_current(:, :), ftcurr(:, :, :), curr(:, :) 46 FLOAT, allocatable :: heat_current(:,:), ftheatcurr(:,:,:), heatcurr(:,:,:) 47 CMPLX, allocatable :: invdielectric(:, :) 48 type(geometry_t) :: geo 49 type(space_t) :: space 50 type(simul_box_t) :: sb 51 type(spectrum_t) :: spectrum 52 type(block_t) :: blk 53 type(batch_t) :: currb, ftcurrb, heatcurrb, ftheatcurrb 54 FLOAT :: ww, curtime, deltat, velcm(1:MAX_DIM), vel0(1:MAX_DIM), current(1:MAX_DIM), integral(1:2), v0 55 integer :: ifreq, max_freq 56 integer :: skip 57 logical :: from_forces 58 character(len=120) :: header 59 FLOAT :: excess_charge, val_charge, qtot 60 61 ! Initialize stuff 62 call global_init(is_serial = .true.) 63 64 call getopt_init(ierr) 65 call getopt_end() 66 67 call parser_init() 68 69 call messages_init() 70 71 call messages_experimental('oct-conductivity') 72 73 call io_init() 74 75 call unit_system_init(global_namespace) 76 77 call spectrum_init(spectrum, global_namespace, default_energy_step = CNST(0.0001), default_max_energy = CNST(1.0)) 78 79 !%Variable ConductivitySpectrumTimeStepFactor 80 !%Type integer 81 !%Default 1 82 !%Section Utilities::oct-conductivity_spectrum 83 !%Description 84 !% In the calculation of the conductivity, it is not necessary 85 !% to read the velocity at every time step. This variable controls 86 !% the integer factor between the simulation time step and the 87 !% time step used to calculate the conductivity. 88 !%End 89 90 call messages_obsolete_variable(global_namespace, 'PropagationSpectrumTimeStepFactor', 'ConductivitySpectrumTimeStepFactor') 91 call parse_variable(global_namespace, 'ConductivitySpectrumTimeStepFactor', 1, skip) 92 if(skip <= 0) call messages_input_error(global_namespace, 'ConductivitySpectrumTimeStepFactor') 93 94 !%Variable ConductivityFromForces 95 !%Type logical 96 !%Default no 97 !%Section Utilities::oct-conductivity_spectrum 98 !%Description 99 !% (Experimental) If enabled, Octopus will attempt to calculate the conductivity from the forces instead of the current. 100 !%End 101 call parse_variable(global_namespace, 'ConductivityFromForces', .false., from_forces) 102 if(from_forces) call messages_experimental('ConductivityFromForces') 103 104 max_freq = spectrum_nenergy_steps(spectrum) 105 106 if (spectrum%end_time < M_ZERO) spectrum%end_time = huge(spectrum%end_time) 107 108 call space_init(space, global_namespace) 109 call geometry_init(geo, global_namespace, space) 110 call simul_box_init(sb, global_namespace, geo, space) 111 112 !We need the total charge 113 call parse_variable(global_namespace, 'ExcessCharge', M_ZERO, excess_charge) 114 call geometry_val_charge(geo, val_charge) 115 qtot = -(val_charge + excess_charge) 116 117 if(from_forces) then 118 119 call messages_write('Info: Reading coordinates from td.general/coordinates') 120 call messages_info() 121 122 ! Opens the coordinates files. 123 iunit = io_open('td.general/coordinates', global_namespace, action='read') 124 125 call io_skip_header(iunit) 126 call spectrum_count_time_steps(global_namespace, iunit, ntime, deltat) 127 call io_close(iunit) 128 129 nvel = geo%natoms*space%dim 130 131 SAFE_ALLOCATE(time(1:ntime)) 132 SAFE_ALLOCATE(velocities(1:nvel, 1:ntime)) 133 134 ! Opens the coordinates files. 135 iunit = io_open('td.general/coordinates', global_namespace, action='read', status='old', die=.false.) 136 137 call io_skip_header(iunit) 138 139 ntime = 1 140 iter = 1 141 142 do 143 read(unit = iunit, iostat = ierr, fmt = *) read_iter, curtime, & 144 ((geo%atom(ii)%x(jj), jj = 1, 3), ii = 1, geo%natoms), & 145 ((geo%atom(ii)%v(jj), jj = 1, 3), ii = 1, geo%natoms) 146 147 curtime = units_to_atomic(units_out%time, curtime) 148 149 if(ierr /= 0 .or. curtime >= spectrum%end_time) then 150 iter = iter - 1 ! last iteration is not valid 151 ntime = ntime - 1 152 exit 153 end if 154 155 ASSERT(iter == read_iter + 1) 156 157 if (curtime >= spectrum%start_time .and. mod(iter, skip) == 0) then 158 159 time(ntime) = curtime 160 ivel = 1 161 do ii = 1, geo%natoms 162 do jj = 1, space%dim 163 velocities(ivel, ntime) = units_to_atomic(units_out%velocity, geo%atom(ii)%v(jj)) 164 ivel = ivel + 1 165 end do 166 end do 167 168 ntime = ntime + 1 169 end if 170 171 iter = iter + 1 !counts number of timesteps (with time larger than zero up to SpecEndTime) 172 end do 173 174 call io_close(iunit) 175 176 call messages_write(' done.') 177 call messages_info() 178 179 else !from_forces 180 181 vel0(:) = M_ONE 182 183 call messages_write('Info: Reading total current from td.general/total_current') 184 call messages_info() 185 186 iunit = io_open('td.general/total_current', global_namespace, action='read') 187 188 if(iunit > 0) then 189 190 call io_skip_header(iunit) 191 call spectrum_count_time_steps(global_namespace, iunit, ntime, deltat) 192 deltat = units_to_atomic(units_out%time, deltat) 193 ntime = ntime + 1 194 195 call io_skip_header(iunit) 196 197 SAFE_ALLOCATE(total_current(1:space%dim, 1:ntime)) 198 SAFE_ALLOCATE(heat_current(1:space%dim, 1:ntime)) 199 SAFE_ALLOCATE(time(1:ntime)) 200 201 do iter = 1, ntime 202 read(iunit, *) read_iter, time(iter), (total_current(ii, iter), ii = 1, space%dim) 203 time(iter) = units_to_atomic(units_out%time, time(iter)) 204 do ii = 1, space%dim 205 total_current(ii, iter) = units_to_atomic(units_out%length/units_out%time, total_current(ii, iter)) 206 end do 207 end do 208 209 call io_close(iunit) 210 211 else 212 213 call messages_write("Cannot find the 'td.general/total_current' file.") 214 call messages_write("Conductivity will only be calculated from the forces") 215 call messages_warning() 216 217 ntime = 1 218 SAFE_ALLOCATE(total_current(1:space%dim, 1:ntime)) 219 SAFE_ALLOCATE(heat_current(1:space%dim, 1:ntime)) 220 SAFE_ALLOCATE(time(1:ntime)) 221 222 total_current(1:space%dim, 1:ntime) = M_ZERO 223 224 end if 225 226 iunit = io_open('td.general/total_heat_current', global_namespace, action='read', status='old', die=.false.) 227 228 if(iunit > 0) then 229 230 if(ntime == 1) then 231 call io_skip_header(iunit) 232 call spectrum_count_time_steps(global_namespace, iunit, ntime, deltat) 233 ntime = ntime + 1 234 end if 235 236 call io_skip_header(iunit) 237 238 do iter = 1, ntime 239 read(iunit, *) read_iter, time(iter), (heat_current(ii, iter), ii = 1, space%dim) 240 end do 241 242 call io_close(iunit) 243 244 else 245 246 call messages_write("Cannot find the 'td.general/heat_current' file.") 247 call messages_write("Thermal conductivity will only be calculated from the forces") 248 call messages_warning() 249 250 heat_current(1:space%dim, 1:ntime) = M_ZERO 251 252 end if 253 end if 254 255 256 SAFE_ALLOCATE(curr(ntime, 1:space%dim)) 257 SAFE_ALLOCATE(heatcurr(ntime, 1:space%dim, 1:1)) 258 integral = M_ZERO 259 260 if(from_forces) iunit = io_open('td.general/current_from_forces', global_namespace, action='write') 261 262 do iter = 1, ntime 263 264 if(from_forces) then 265 if(iter == 1) then 266 vel0 = M_ZERO 267 ivel = 1 268 do ii = 1, geo%natoms 269 do jj = 1, space%dim 270 vel0(jj) = vel0(jj) + velocities(ivel, iter)/TOFLOAT(geo%natoms) 271 ivel = ivel + 1 272 end do 273 end do 274 end if 275 276 velcm = M_ZERO 277 current = M_ZERO 278 279 ivel = 1 280 do ii = 1, geo%natoms 281 do jj = 1, space%dim 282 velcm(jj) = velcm(jj) + velocities(ivel, iter)/TOFLOAT(geo%natoms) 283 current(jj) = current(jj) + species_mass(geo%atom(ii)%species)/sb%rcell_volume*(velocities(ivel, iter) - vel0(jj)) 284 ivel = ivel + 1 285 end do 286 end do 287 288 integral(1) = integral(1) + deltat/vel0(1)*(vel0(1)*qtot/sb%rcell_volume + current(1)) 289 integral(2) = integral(2) + deltat/vel0(1)*(vel0(1)*qtot/sb%rcell_volume - total_current(1, iter)/sb%rcell_volume) 290 291 curr(iter, 1:space%dim) = vel0(1:space%dim)*qtot/sb%rcell_volume + current(1:space%dim) 292 293 else 294 curr(iter, 1:space%dim) = total_current(1:space%dim, iter)/sb%rcell_volume 295 heatcurr(iter,1:space%dim, 1) = heat_current(1:space%dim, iter)/sb%rcell_volume 296 end if 297 298 if(from_forces) write(iunit,*) iter, iter*deltat, curr(iter, 1:space%dim) 299 300 end do 301 302 if(from_forces) call io_close(iunit) 303 304 ! Find out the iteration numbers corresponding to the time limits. 305 call spectrum_fix_time_limits(spectrum, ntime, deltat, istart, iend, max_freq) 306 istart = max(1, istart) 307 energy_steps = spectrum_nenergy_steps(spectrum) 308 309 SAFE_ALLOCATE(ftcurr(1:energy_steps, 1:space%dim, 1:2)) 310 311 call batch_init(currb, 1, 1, space%dim, curr) 312 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, deltat, currb) 313 314 call batch_init(ftcurrb, 1, 1, space%dim, ftcurr(:, :, 1)) 315 call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_COS, spectrum%noise, & 316 istart, iend, spectrum%start_time, deltat, currb, spectrum%min_energy, spectrum%max_energy, & 317 spectrum%energy_step, ftcurrb) 318 call ftcurrb%end() 319 320 call batch_init(ftcurrb, 1, 1, space%dim, ftcurr(:, :, 2)) 321 call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_SIN, spectrum%noise, & 322 istart, iend, spectrum%start_time, deltat, currb, spectrum%min_energy, spectrum%max_energy, & 323 spectrum%energy_step, ftcurrb) 324 325 call ftcurrb%end() 326 call currb%end() 327 328 329 !and print the spectrum 330 iunit = io_open('td.general/conductivity', global_namespace, action='write') 331 332 333 write(unit = iunit, iostat = ierr, fmt = '(a)') & 334 '###########################################################################################################################' 335 write(unit = iunit, iostat = ierr, fmt = '(8a)') '# HEADER' 336 write(unit = iunit, iostat = ierr, fmt = '(a,a,a)') & 337 '# Energy [', trim(units_abbrev(units_out%energy)), '] Conductivity [a.u.] ReX ImX ReY ImY ReZ ImZ' 338 write(unit = iunit, iostat = ierr, fmt = '(a)') & 339 '###########################################################################################################################' 340 341 v0 = sqrt(sum(vel0(1:space%dim)**2)) 342 if(.not. from_forces .or. v0 < epsilon(v0)) v0 = CNST(1.0) 343 344 if( .not. from_forces .and. parse_is_defined(global_namespace, 'GaugeVectorField')) then 345 if(parse_block(global_namespace, 'GaugeVectorField', blk) == 0) then 346 do ii = 1, space%dim 347 call parse_block_float(blk, 0, ii - 1, vel0(ii)) 348 end do 349 call parse_block_end(blk) 350 end if 351 vel0(1:space%dim) = vel0(1:space%dim) / P_C 352 v0 = sqrt(sum(vel0(1:space%dim)**2)) 353 end if 354 355 356 do ifreq = 1, energy_steps 357 ww = spectrum%energy_step*(ifreq - 1) + spectrum%min_energy 358 write(unit = iunit, iostat = ierr, fmt = '(7e20.10)') units_from_atomic(units_out%energy, ww), & 359 transpose(ftcurr(ifreq, 1:space%dim, 1:2)/v0) 360 end do 361 362 call io_close(iunit) 363 364 365 !Compute the inverse dielectric function from the conductivity 366 ! We have \chi = -i \sigma / \omega 367 ! and \epsilon^-1 = 1 + 4 \pi \chi 368 SAFE_ALLOCATE(invdielectric(1:space%dim, 1:energy_steps)) 369 do ifreq = 1, energy_steps 370 ww = max((ifreq-1)*spectrum%energy_step + spectrum%min_energy, M_EPSILON) 371 372 invdielectric(1:space%dim, ifreq) = (vel0(1:space%dim) + M_FOUR * M_PI * & 373 TOCMPLX(ftcurr(ifreq, 1:space%dim, 2),-ftcurr(ifreq, 1:space%dim, 1)) / ww )/v0 374 375 end do 376 377 out_file = io_open('td.general/inverse_dielectric_function_from_current', global_namespace, action='write') 378 write(header, '(7a15)') '# energy', 'Re x', 'Im x', 'Re y', 'Im y', 'Re z', 'Im z' 379 write(out_file,'(a)') trim(header) 380 do ifreq = 1, energy_steps 381 ww = (ifreq-1)*spectrum%energy_step + spectrum%min_energy 382 write(out_file, '(7e15.6)') ww, & 383 TOFLOAT(invdielectric(1, ifreq)), aimag(invdielectric(1, ifreq)), & 384 TOFLOAT(invdielectric(2, ifreq)), aimag(invdielectric(2, ifreq)), & 385 TOFLOAT(invdielectric(3, ifreq)), aimag(invdielectric(3, ifreq)) 386 end do 387 call io_close(out_file) 388 389 SAFE_DEALLOCATE_A(ftcurr) 390 SAFE_DEALLOCATE_A(invdielectric) 391 392 393!!!!!!!!!!!!!!!!!!!!! 394 SAFE_ALLOCATE(ftheatcurr(1:energy_steps, 1:3, 1:2)) 395 396 ftheatcurr = M_ONE 397 398 call batch_init(heatcurrb, 3, 1, 1, heatcurr) 399 400 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, 1, ntime, M_ZERO, deltat, heatcurrb) 401 402 call batch_init(ftheatcurrb, 3, 1, 1, ftheatcurr) 403 call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_COS, spectrum%noise, & 404 1, ntime, M_ZERO, deltat, heatcurrb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftheatcurrb) 405 call ftheatcurrb%end() 406 407 call batch_init(ftheatcurrb, 3, 1, 1, ftheatcurr(:, :, 2:2)) 408 call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_SIN, spectrum%noise, & 409 1, ntime, M_ZERO, deltat, heatcurrb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftheatcurrb) 410 call ftheatcurrb%end() 411 412 call heatcurrb%end() 413 414 415 !and print the spectrum 416 iunit = io_open('td.general/heat_conductivity', global_namespace, action='write') 417 418 419 write(unit = iunit, iostat = ierr, fmt = '(a)') & 420 '###########################################################################################################################' 421 write(unit = iunit, iostat = ierr, fmt = '(8a)') '# HEADER' 422 write(unit = iunit, iostat = ierr, fmt = '(a,a,a)') & 423 '# Energy [', trim(units_abbrev(units_out%energy)), '] Conductivity [a.u.] ReX ImX ReY ImY ReZ ImZ' 424 write(unit = iunit, iostat = ierr, fmt = '(a)') & 425 '###########################################################################################################################' 426 427 v0 = sqrt(sum(vel0(1:space%dim)**2)) 428 if(.not. from_forces .or. v0 < epsilon(v0)) v0 = CNST(1.0) 429 430 do ifreq = 1, energy_steps 431 ww = spectrum%energy_step*(ifreq - 1) + spectrum%min_energy 432 write(unit = iunit, iostat = ierr, fmt = '(7e20.10)') units_from_atomic(units_out%energy, ww), & 433 transpose(ftheatcurr(ifreq, 1:3, 1:2)/v0) 434 !print *, ifreq, ftheatcurr(ifreq, 1:3, 1:2) 435 end do 436 437 call io_close(iunit) 438 439 SAFE_DEALLOCATE_A(ftheatcurr) 440!!!!!!!!!!!!!!!!!!!!!!!!!!! 441 442 443 call simul_box_end(sb) 444 call geometry_end(geo) 445 call space_end(space) 446 447 SAFE_DEALLOCATE_A(time) 448 449 call io_end() 450 call messages_end() 451 452 call parser_end() 453 call global_end() 454 455 end program conductivity 456 457!! Local Variables: 458!! mode: f90 459!! coding: utf-8 460!! End: 461