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