Octopus
oscillator_strength.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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
22 use debug_oct_m
23 use global_oct_m
24 use io_oct_m
25 use kick_oct_m
26 use, intrinsic :: iso_fortran_env
27 use iso_c_binding
32 use parser_oct_m
35 use string_oct_m
36 use unit_oct_m
38
39 implicit none
40
41 integer, parameter :: &
42 SINE_TRANSFORM = 1, &
44
46 integer :: n_multipoles
47 integer, allocatable :: l(:), m(:)
48 real(real64), allocatable :: weight(:)
49 end type local_operator_t
50
51 integer :: observable(2)
52 type(unit_system_t) :: units
53 real(real64), allocatable :: ot(:)
54 type(kick_t) :: kick
55 integer :: time_steps
56 real(real64) :: total_time
57 integer :: mode
58 real(real64) :: dt
59
60contains
61
62 ! ---------------------------------------------------------
63 subroutine ft(omega, power)
64 real(real64), intent(in) :: omega
65 real(real64), intent(out) :: power
66
67 integer :: j
68 real(real64) :: x
69 power = m_zero
70
71 push_sub(ft)
72
73 select case (mode)
74
75 case (sine_transform)
76
77 do j = 0, time_steps
78 x = sin(omega*j*dt)
79 power = power + x*ot(j)
80 end do
81 power = power*dt / (dt*time_steps)
82
83 case (cosine_transform)
84
85 do j = 0, time_steps
86 x = cos(omega*j*dt)
87 power = power + x*ot(j)
88 end do
89 power = power*dt / (dt*time_steps)
90
91 end select
92
93 pop_sub(ft)
94 end subroutine ft
95
96
97 ! ---------------------------------------------------------
98 subroutine ft2(omega, power)
99 real(real64), intent(in) :: omega
100 real(real64), intent(out) :: power
101
102 integer :: j
103 real(real64) :: x
104 power = m_zero
105
106 push_sub(ft2)
107
108 select case (mode)
109
110 case (sine_transform)
111
112 do j = 0, time_steps
113 x = sin(omega*j*dt)
114 power = power + x*ot(j)
115 end do
116 ! The function should be *minus* the sine Fourier transform, since this
117 ! is the function to be minimized.
118 power = - (power*dt / (dt*time_steps))**2
119
120 case (cosine_transform)
121
122 do j = 0, time_steps
123 x = cos(omega*j*dt)
124 power = power + x*ot(j)
125 end do
126 power = - (power*dt / (dt*time_steps))**2
127
128 end select
129
130 pop_sub(ft2)
131 end subroutine ft2
132
133
134 ! ---------------------------------------------------------
135 subroutine local_operator_copy(o, i)
136 type(local_operator_t), intent(inout) :: o
137 type(local_operator_t), intent(inout) :: i
138
139 integer :: j
143 o%n_multipoles = i%n_multipoles
144 safe_allocate( o%l(1:o%n_multipoles))
145 safe_allocate( o%m(1:o%n_multipoles))
146 safe_allocate(o%weight(1:o%n_multipoles))
148 do j = 1, o%n_multipoles
149 o%l(j) = i%l(j)
150 o%m(j) = i%m(j)
151 o%weight(j) = i%weight(j)
152 end do
154 pop_sub(local_operator_copy)
155 end subroutine local_operator_copy
156
157 ! ---------------------------------------------------------
158 subroutine read_resonances_file(order, ffile, namespace, search_interval, final_time, nfrequencies)
159 integer, intent(inout) :: order
160 character(len=*), intent(in) :: ffile
161 type(namespace_t), intent(in):: namespace
162 real(real64), intent(inout) :: search_interval
163 real(real64), intent(in) :: final_time
164 integer, intent(in) :: nfrequencies
165
166 real(real64) :: dummy, leftbound, rightbound, w, power, dw
167 integer :: iunit, nresonances, ios, i, j, k, npairs, nspin, order_in_file, nw_subtracted, ierr
168 logical :: file_exists
169 real(real64), allocatable :: wij(:), omega(:), c0i(:)
170
171 push_sub(read_resonances_file)
172
173 if (order /= 2) then
174 write(message(1),'(a)') 'The run mode #3 is only compatible with the analysis of the'
175 write(message(2),'(a)') 'second-order response.'
176 call messages_fatal(2)
177 end if
178
179 ! First, let us check that the file "ot" exists.
180 inquire(file="ot", exist = file_exists)
181 if (.not. file_exists) then
182 write(message(1),'(a)') "Could not find 'ot' file."
183 call messages_fatal(1)
184 end if
185
186 ! Now, we should find out which units the file "ot" has.
187 call unit_system_from_file(units, "ot", namespace, ierr)
188 if (ierr /= 0) then
189 write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
190 call messages_fatal(1)
191 end if
192
194
195 iunit = io_open(trim(ffile), namespace, action='read', status='old')
196
197 call io_skip_header(iunit)
198 ! Count the number of resonances
199 nresonances = 0
200 do
201 read(iunit, *, iostat = ios) dummy, dummy
202 if (ios /= 0) exit
203 nresonances = nresonances + 1
204 end do
205
206 npairs = (nresonances*(nresonances-1))/2
207
208 safe_allocate(omega(1:nresonances))
209 safe_allocate( c0i(1:nresonances))
210 safe_allocate(wij(1:npairs))
211
212 call io_skip_header(iunit)
213 do i = 1, nresonances
214 read(iunit, *) omega(i), c0i(i)
215 end do
216
217 k = 1
218 do i = 1, nresonances
219 do j = i + 1, nresonances
220 wij(k) = omega(j) - omega(i)
221 k = k + 1
222 end do
223 end do
224
225 if (search_interval > m_zero) then
226 search_interval = units_to_atomic(units%energy, search_interval)
227 else
228 search_interval = m_half
229 end if
231 call read_ot(namespace, nspin, order_in_file, nw_subtracted)
232
233 if (order_in_file /= order) then
234 write(message(1), '(a)') 'The ot file should contain the second-order response in this run mode.'
235 call messages_fatal(1)
236 end if
237
238 if (final_time > m_zero) then
239 total_time = units_to_atomic(units%time, final_time)
240 if (total_time > dt*time_steps) then
241 total_time = dt*time_steps
242 write(message(1), '(a)') 'The requested total time to process is larger than the time available in the input file.'
243 write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
244 units_from_atomic(units%time, total_time), units_abbrev(units%time)
245 call messages_warning(2, namespace=namespace)
246 end if
247 time_steps = int(total_time / dt)
248 total_time = time_steps * dt
249 else
250 total_time = dt*time_steps
251 end if
252
253 dw = (rightbound-leftbound) / (nfrequencies - 1)
254
255 ! First, subtract zero resonance...
256 w = m_zero
257 call resonance_second_order(w, power, nw_subtracted, leftbound, rightbound, m_zero, m_zero)
258 call modify_ot(time_steps, dt, order, ot, w, power)
259 nw_subtracted = nw_subtracted + 1
260
261 ! Then, get all the others...
262 k = 1
263 do i = 1, nresonances
264 do j = i + 1, nresonances
265 leftbound = wij(k) - search_interval
266 rightbound = wij(k) + search_interval
267 call find_resonance(wij(k), leftbound, rightbound, nfrequencies)
268 call resonance_second_order(wij(k), power, nw_subtracted, leftbound, rightbound, c0i(i), c0i(j))
269 call modify_ot(time_steps, dt, order, ot, wij(k), power)
270 nw_subtracted = nw_subtracted + 1
271 k = k + 1
272 end do
273 end do
274
275 safe_deallocate_a(wij)
276 safe_deallocate_a(c0i)
277 safe_deallocate_a(omega)
278 call io_close(iunit)
279
280 pop_sub(read_resonances_file)
281 end subroutine read_resonances_file
282 ! ---------------------------------------------------------
283
284
285 ! ---------------------------------------------------------
286 subroutine analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, &
287 damping, namespace)
288 integer, intent(inout) :: order
289 real(real64), intent(inout) :: omega
290 real(real64), intent(inout) :: search_interval
291 real(real64), intent(inout) :: final_time
292 integer, intent(inout) :: nresonances
293 integer, intent(inout) :: nfrequencies
294 real(real64), intent(in) :: damping
295 type(namespace_t), intent(in) :: namespace
296
297 real(real64) :: leftbound, rightbound, dw, power
298 real(real64), allocatable :: w(:), c0I2(:)
299 integer :: nspin, i, ierr, order_in_file, nw_subtracted
300 logical :: file_exists
301
302 push_sub(analyze_signal)
303
304 ! First, let us check that the file "ot" exists.
305 inquire(file="ot", exist = file_exists)
306 if (.not. file_exists) then
307 write(message(1),'(a)') "Could not find 'ot' file."
308 call messages_fatal(1)
309 end if
310
311 ! Now, we should find out which units the file "ot" has.
312 call unit_system_from_file(units, "ot", namespace, ierr)
313 if (ierr /= 0) then
314 write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
315 call messages_fatal(1)
316 end if
317
318 if (omega > m_zero) then
319 omega = units_to_atomic(units%energy, omega)
320 else
321 omega = m_half
322 end if
323
324 if (search_interval > m_zero) then
325 search_interval = units_to_atomic(units%energy, search_interval)
326 else
327 search_interval = m_half
328 end if
329
330 leftbound = omega - search_interval
331 rightbound = omega + search_interval
332
333 call read_ot(namespace, nspin, order_in_file, nw_subtracted)
334
335 if (order_in_file /= order) then
336 write(message(1), '(a)') 'Internal error in analyze_signal'
337 call messages_fatal(1, namespace=namespace)
338 end if
339
340 if (mod(order, 2) == 1) then
341 mode = sine_transform
342 else
343 mode = cosine_transform
344 end if
345
346 if (final_time > m_zero) then
347 total_time = units_to_atomic(units%time, final_time)
348 if (total_time > dt*time_steps) then
349 total_time = dt*time_steps
350 write(message(1), '(a)') 'The requested total time to process is larger than the time available in the input file.'
351 write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
352 units_from_atomic(units%time, total_time), units_abbrev(units%time)
353 call messages_warning(2, namespace=namespace)
354 end if
355 time_steps = int(total_time / dt)
356 total_time = time_steps * dt
357 else
358 total_time = dt*time_steps
359 end if
360
361 dw = (rightbound-leftbound) / (nfrequencies - 1)
362
363 safe_allocate( w(1:nresonances))
364 safe_allocate(c0i2(1:nresonances))
365 w = m_zero
366 c0i2 = m_zero
367
368 i = 1
369 do
370 if (nw_subtracted >= nresonances) exit
371
372 if (mode == cosine_transform .and. nw_subtracted == 0) then
373 omega = m_zero
374 else
375 call find_resonance(omega, leftbound, rightbound, nfrequencies)
376 end if
377
378 select case (order)
379 case (1)
380 call resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound)
381 case (2)
382 call resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, m_zero, m_zero)
383 end select
384
385 w(i) = omega
386 c0i2(i) = power
387
388 call modify_ot(time_steps, dt, order, ot, omega, power)
389
390 nw_subtracted = nw_subtracted + 1
391 i = i + 1
392 end do
393
394 select case (order)
395 case (1)
396 call write_polarizability(namespace, nfrequencies, nresonances, dw, w, c0i2, damping)
397 end select
398
399 safe_deallocate_a(ot)
400 safe_deallocate_a(w)
401 safe_deallocate_a(c0i2)
402
403 pop_sub(analyze_signal)
404 end subroutine analyze_signal
405 ! ---------------------------------------------------------
406
407
408 ! ---------------------------------------------------------
412 subroutine write_polarizability(namespace, nfrequencies, nresonances, dw, w, c0I2, gamma)
413 type(namespace_t), intent(in) :: namespace
414 integer, intent(in) :: nfrequencies, nresonances
415 real(real64), intent(in) :: dw
416 real(real64), intent(in) :: w(nresonances), c0I2(nresonances)
417 real(real64), intent(in) :: gamma
418
419 integer :: iunit, i, j
420 real(real64) :: e
421 complex(real64) :: pol
422
423 push_sub(write_polarizability)
424
425 iunit = io_open('polarizability', namespace, action = 'write', status='replace', die=.false.)
426 write(iunit, '(a)') '# Polarizability file. Generated using the SOS formula with the following data:'
427 write(iunit, '(a)') '#'
428
429 do i = 1, nresonances
430 write(iunit, '(a1,3e20.12)') '#', w(i), sqrt(abs(c0i2(i))), c0i2(i)
431 end do
432
433 write(iunit, '(a10,f12.6)') '# Gamma = ', gamma
434 write(iunit, '(a)') '#'
435
436 do i = 1, nfrequencies
437 e = (i - 1) * dw
438 pol = m_z0
439 do j = 1, nresonances
440 pol = pol + c0i2(j) * (m_one / (w(j) - e - m_zi * gamma) + m_one/(w(j) + e + m_zi * gamma))
441 end do
442 write(iunit, '(3e20.12)') e, pol
443 end do
444
445 call io_close(iunit)
446
447 pop_sub(write_polarizability)
448 end subroutine write_polarizability
449 ! ---------------------------------------------------------
450
451
452 ! ---------------------------------------------------------
453 ! \todo This subroutine should be simplified.
454 subroutine find_resonance(omega, leftbound, rightbound, nfrequencies)
455 real(real64), intent(inout) :: omega
456 real(real64), intent(in) :: leftbound, rightbound
457 integer, intent(in) :: nfrequencies
458
459 integer :: i, ierr
460 real(real64) :: dw, w, aw, min_aw, min_w, omega_orig, min_w1, min_w2
461 real(real64), allocatable :: warray(:), tarray(:)
462
463 push_sub(find_resonance)
464
465 safe_allocate(warray(1:nfrequencies))
466 safe_allocate(tarray(1:nfrequencies))
467
468 warray = m_zero
469 tarray = m_zero
470 dw = (rightbound-leftbound) / (nfrequencies - 1)
471 do i = 1, nfrequencies
472 w = leftbound + (i-1)*dw
473 warray(i) = w
474 call ft2(w, aw)
475 tarray(i) = aw
476 end do
477
478 min_w = omega
479 min_aw = m_zero
480 do i = 1, nfrequencies
481 w = leftbound + (i-1)*dw
482 if (tarray(i) < min_aw) then
483 min_aw = tarray(i)
484 min_w = w
485 end if
486 end do
487
488 omega_orig = omega
489 omega = min_w
490 min_w1 = min_w - 2*dw
491 min_w2 = min_w + 2*dw
492 call loct_1dminimize(min_w1, min_w2, omega, ft2, ierr)
493
494 if (ierr /= 0) then
495 write(message(1),'(a)') 'Could not find a maximum.'
496 write(message(2),'(a)')
497 write(message(3), '(a,f12.8,a,f12.8,a)') ' Search interval = [', &
498 units_from_atomic(units%energy, leftbound), ',', units_from_atomic(units%energy, rightbound), ']'
499 write(message(4), '(a,f12.4,a)') ' Search discretization = ', &
500 units_from_atomic(units%energy, dw), ' '//trim(units_abbrev(units%energy))
501 call messages_fatal(4)
502 end if
503
504 safe_deallocate_a(warray)
505 safe_deallocate_a(tarray)
506
508 end subroutine find_resonance
509 ! ---------------------------------------------------------
510
511
512 ! ---------------------------------------------------------
513 subroutine resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound)
514 real(real64), intent(in) :: omega
515 real(real64), intent(out) :: power
516 integer, intent(in) :: nw_subtracted
517 real(real64), intent(in) :: dw, leftbound, rightbound
518
519 push_sub(resonance_first_order)
520
521 call ft(omega, power)
522
523 select case (mode)
524 case (sine_transform)
525 power = power / (m_one - sin(m_two*omega*total_time)/(m_two*omega*total_time))
526 case (cosine_transform)
527 call messages_not_implemented("resonance first order cosine transform")
528 end select
529
530 write(message(1), '(a)') '******************************************************************'
531 write(message(2), '(a,i3)') 'Resonance #', nw_subtracted + 1
532 write(message(3), '(a,f12.8,a,f12.8,a)') 'omega = ', units_from_atomic(units_out%energy, omega), &
533 ' '//trim(units_abbrev(units_out%energy))//' = ', omega, ' Ha'
534 write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**2, power), &
535 ' '//trim(units_abbrev(units_out%length**2))//' =', power, ' b^2'
536 write(message(5), '(a,f12.8,a,f12.8,a)') '<0|P|I> = ', units_from_atomic(units_out%length, sqrt(abs(power))), &
537 ' '//trim(units_abbrev(units_out%length))//' = ', sqrt(abs(power)),' b'
538 write(message(6), '(a,f12.8)') 'f[O->I] = ', m_two*omega*power
539 write(message(7), '(a)')
540 write(message(8), '(a,f12.8,a,f12.8,a)') ' Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', &
541 units_from_atomic(units_out%energy, rightbound), ']'
542 write(message(9), '(a,f12.4,a)') ' Search discretization = ', units_from_atomic(units_out%energy, dw), &
543 ' '//trim(units_abbrev(units_out%energy))
544 write(message(10), '(a)') '******************************************************************'
545 write(message(11), '(a)')
546 call messages_info(11)
547
548 pop_sub(resonance_first_order)
550 end subroutine resonance_first_order
551 ! ---------------------------------------------------------
552
553
554 ! ---------------------------------------------------------
555 subroutine resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, c01, c02)
556 real(real64), intent(in) :: omega
557 real(real64), intent(out) :: power
558 integer, intent(in) :: nw_subtracted
559 real(real64), intent(in) :: leftbound, rightbound
560 real(real64), intent(in) :: c01, c02
561
562 push_sub(resonance_second_order)
563
564 call ft(omega, power)
565 select case (mode)
566 case (sine_transform)
567 power = power / (m_one - sin(m_two*omega*total_time)/(m_two*omega*total_time))
568 case (cosine_transform)
569 ! WARNING: there is some difference between the omega=0 case and the rest.
570 if (abs(omega) > m_epsilon) then
571 power = power / (m_one + sin(m_two*omega*total_time)/(m_two*omega*total_time))
572 else
573 power = power / m_two
574 end if
575 end select
576
577 write(message(1), '(a)') '******************************************************************'
578 write(message(2), '(a,i3)') 'Resonance #', nw_subtracted + 1
579 write(message(3), '(a,f12.8,a,f12.8,a)') 'omega = ', units_from_atomic(units_out%energy, omega), &
580 ' '//trim(units_abbrev(units_out%energy))//' = ', omega, ' Ha'
581 write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**3, power), &
582 ' '//trim(units_abbrev(units_out%length**3))//' = ', power, ' b^3'
583 call messages_info(4)
584
585 if (abs(c01*c02) > m_epsilon) then
586 write(message(1), '(a,f12.8)') ' C(omega)/(C0i*C0j) = ', power / (c01 * c02)
587 call messages_info(1)
588 end if
589
590 write(message(1), '(a)')
591 write(message(2), '(a,f12.8,a,f12.8,a)') ' Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', &
592 units_from_atomic(units_out%energy, rightbound), ']'
593 write(message(3), '(a)') '******************************************************************'
594 write(message(4), '(a)')
595 call messages_info(4)
596
598 end subroutine resonance_second_order
599 ! ---------------------------------------------------------
600
601
602 ! ---------------------------------------------------------
603 subroutine generate_signal(namespace, order, observable)
604 type(namespace_t), intent(in) :: namespace
605 integer, intent(in) :: order
606 integer, intent(in) :: observable(2)
607
608 logical :: file_exists
609 integer :: i, j, nspin, time_steps, lmax, nfiles, k, add_lm, l, m, max_add_lm
610 integer, allocatable :: iunit(:)
611 real(real64) :: dt, lambda, dump, o0
612 type(unit_t) :: mp_unit
613 real(real64), allocatable :: q(:), mu(:), qq(:, :), c(:)
614 character(len=20) :: filename
615 type(kick_t) :: kick
616 type(unit_system_t) :: units
617 real(real64), allocatable :: multipole(:, :, :), ot(:), dipole(:, :)
618 type(local_operator_t) :: kick_operator
619 type(local_operator_t) :: obs
620
621 push_sub(generate_signal)
622
623 ! Find out how many files do we have
624 nfiles = 0
625 do
626 write(filename,'(a11,i1)') 'multipoles.', nfiles+1
627 inquire(file=trim(filename), exist = file_exists)
628 if (.not. file_exists) exit
629 nfiles = nfiles + 1
630 end do
631
632 ! WARNING: Check that order is smaller or equal to nfiles
633 if (nfiles == 0) then
634 write(message(1),'(a)') 'No multipoles.x file was found'
635 call messages_fatal(1, namespace=namespace)
636 end if
637 if (order > nfiles) then
638 write(message(1),'(a)') 'The order that you ask for is higher than the number'
639 write(message(2),'(a)') 'of multipoles.x file that you supply.'
640 call messages_fatal(2, namespace=namespace)
641 end if
642
643 ! Open the files.
644 safe_allocate(iunit(1:nfiles))
645 do j = 1, nfiles
646 write(filename,'(a11,i1)') 'multipoles.', j
647 iunit(j) = io_open(trim(filename), namespace, action='read', status='old')
648 end do
649
650 safe_allocate( q(1:nfiles))
651 safe_allocate(mu(1:nfiles))
652 safe_allocate(qq(1:nfiles, 1:nfiles))
653 safe_allocate( c(1:nfiles))
654
655 c = m_zero
656 c(order) = m_one
657
658 ! Get the basic info from the first file
659 call spectrum_mult_info(namespace, iunit(1), nspin, kick, time_steps, dt, units, lmax=lmax)
660
661 ! Sets the kick operator...
662 if (kick%n_multipoles > 0) then
663 kick_operator%n_multipoles = kick%n_multipoles
664 safe_allocate( kick_operator%l(1:kick_operator%n_multipoles))
665 safe_allocate( kick_operator%m(1:kick_operator%n_multipoles))
666 safe_allocate(kick_operator%weight(1:kick_operator%n_multipoles))
667 do i = 1, kick_operator%n_multipoles
668 kick_operator%l(i) = kick%l(i)
669 kick_operator%m(i) = kick%m(i)
670 kick_operator%weight(i) = kick%weight(i)
671 end do
672 else
673 kick_operator%n_multipoles = 3
674 safe_allocate( kick_operator%l(1:kick_operator%n_multipoles))
675 safe_allocate( kick_operator%m(1:kick_operator%n_multipoles))
676 safe_allocate(kick_operator%weight(1:kick_operator%n_multipoles))
677 kick_operator%l(1:3) = 1
678 kick_operator%m(1) = -1
679 kick_operator%m(2) = 0
680 kick_operator%m(3) = 1
681 ! WARNING: not sure if m = -1 => y, and m = 1 => x. What is the convention?
682 kick_operator%weight(1) = -sqrt((m_four*m_pi)/m_three) * kick%pol(2, kick%pol_dir)
683 kick_operator%weight(2) = sqrt((m_four*m_pi)/m_three) * kick%pol(3, kick%pol_dir)
684 kick_operator%weight(3) = -sqrt((m_four*m_pi)/m_three) * kick%pol(1, kick%pol_dir)
685 end if
686
687 ! Sets the observation operator
688 select case (observable(1))
689 case (-1)
690 ! This means that the "observation operator" should be equal
691 ! to the "perturbation operator", i.e., the kick.
692 call local_operator_copy(obs, kick_operator)
693 case (0)
694 ! This means that the observable is the dipole operator; observable(2) determines
695 ! if it is x, y or z.
696 obs%n_multipoles = 1
697 safe_allocate(obs%l(1:1))
698 safe_allocate(obs%m(1:1))
699 safe_allocate(obs%weight(1:1))
700 obs%l(1) = 1
701 select case (observable(2))
702 case (1)
703 obs%m(1) = -1
704 obs%weight(1) = -sqrt((m_four*m_pi)/m_three)
705 case (2)
706 obs%m(1) = 1
707 obs%weight(1) = sqrt((m_four*m_pi)/m_three)
708 case (3)
709 obs%m(1) = 0
710 obs%weight(1) = -sqrt((m_four*m_pi)/m_three)
711 end select
712 case default
713 ! This means that the observation operator is (l,m) = (observable(1), observable(2))
714 obs%n_multipoles = 1
715 safe_allocate(obs%l(1:1))
716 safe_allocate(obs%m(1:1))
717 safe_allocate(obs%weight(1:1))
718 obs%weight(1) = m_one
719 obs%l(1) = observable(1)
720 obs%m(1) = observable(2)
721 end select
722
723 lambda = abs(kick%delta_strength)
724 q(1) = kick%delta_strength / lambda
725
726 do j = 2, nfiles
727 call spectrum_mult_info(namespace, iunit(j), nspin, kick, time_steps, dt, units, lmax=lmax)
728 q(j) = kick%delta_strength / lambda
729 end do
730
731 do i = 1, nfiles
732 do j = 1, nfiles
733 qq(i, j) = q(j)**i
734 end do
735 end do
736
737 call lalg_inverse(nfiles, qq, 'dir')
738
739 mu = matmul(qq, c)
740
741 if (kick%n_multipoles > 0) then
742 lmax = maxval(kick%l(1:obs%n_multipoles))
743 max_add_lm = (lmax+1)**2-1
744 safe_allocate(multipole(1:max_add_lm, 0:time_steps, 1:nspin))
745 ! The units have nothing to do with the perturbing kick??
746 mp_unit = units%length**kick%l(1)
747 else
748 max_add_lm = 3
749 safe_allocate(multipole(1:3, 0:time_steps, 1:nspin))
750 mp_unit = units%length
751 end if
752 safe_allocate(ot(0:time_steps))
753 multipole = m_zero
754 ot = m_zero
755 o0 = m_zero
756
757 safe_allocate(dipole(1:3, 1:nspin))
758
759 do j = 1, nfiles
760 call io_skip_header(iunit(j))
761
762 do i = 0, time_steps
763 select case (nspin)
764 case (1)
765 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm)
766 case (2)
767 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), &
768 dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm)
769 case (4)
770 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), &
771 dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm), &
772 dump, (multipole(add_lm, i, 3), add_lm = 1, max_add_lm), &
773 dump, (multipole(add_lm, i, 4), add_lm = 1, max_add_lm)
774 end select
775 multipole(1:max_add_lm, i, :) = units_to_atomic(mp_unit, multipole(1:max_add_lm, i, :))
776
777 ! The dipole is treated differently in the multipoles file: first of all,
778 ! the program should have printed *minus* the dipole operator.
779 dipole(1:3, 1:nspin) = - multipole(1:3, i, 1:nspin)
780 ! And then it contains the "cartesian" dipole, opposed to the spherical dipole:
781 multipole(1, i, 1:nspin) = -sqrt(m_three/(m_four*m_pi)) * dipole(2, 1:nspin)
782 multipole(2, i, 1:nspin) = sqrt(m_three/(m_four*m_pi)) * dipole(3, 1:nspin)
783 multipole(3, i, 1:nspin) = -sqrt(m_three/(m_four*m_pi)) * dipole(1, 1:nspin)
784
785 dump = m_zero
786 do k = 1, obs%n_multipoles
787 add_lm = 1
788 l = 1
789 lcycle: do
790 do m = -l, l
791 if (l == obs%l(k) .and. m == obs%m(k)) exit lcycle
792 add_lm = add_lm + 1
793 end do
794 l = l + 1
795 end do lcycle
796 ! Warning: it should not add the nspin components?
797 dump = dump + obs%weight(k) * sum(multipole(add_lm, i, 1:nspin))
798 end do
799
800 if (i == 0) o0 = dump
801 ot(i) = ot(i) + mu(j)*(dump - o0)
802 end do
803
804 end do
805
806 ot = ot / lambda**order
807
808 call write_ot(namespace, nspin, time_steps, dt, kick, order, ot(0:time_steps), observable)
809
810 ! Close files and exit.
811 do j = 1, nfiles
812 call io_close(iunit(j))
813 end do
814 safe_deallocate_a(iunit)
815 safe_deallocate_a(q)
816 safe_deallocate_a(mu)
817 safe_deallocate_a(qq)
818 safe_deallocate_a(c)
819 safe_deallocate_a(ot)
820 safe_deallocate_a(multipole)
821 safe_deallocate_a(dipole)
822
823 pop_sub(generate_signal)
824 end subroutine generate_signal
825 ! ---------------------------------------------------------
826
827
828 ! ---------------------------------------------------------
829 subroutine modify_ot(time_steps, dt, order, ot, omega, power)
830 integer, intent(in) :: time_steps
831 real(real64), intent(in) :: dt
832 integer, intent(in) :: order
833 real(real64), intent(inout) :: ot(0:time_steps)
834 real(real64), intent(in) :: omega
835 real(real64), intent(in) :: power
836
837 integer :: i
838
839 push_sub(modify_ot)
840
841 select case (mod(order, 2))
842 case (1)
843 do i = 0, time_steps
844 ot(i) = ot(i) - m_two * power * sin(omega*i*dt)
845 end do
846 case (0)
847 if (abs(omega) <= m_epsilon) then
848 do i = 0, time_steps
849 ot(i) = ot(i) - power * cos(omega*i*dt)
850 end do
851 else
852 do i = 0, time_steps
853 ot(i) = ot(i) - m_two * power * cos(omega*i*dt)
854 end do
855 end if
856 end select
857
858 pop_sub(modify_ot)
859 end subroutine modify_ot
860 ! ---------------------------------------------------------
861
862
863 ! ---------------------------------------------------------
864 subroutine write_ot(namespace, nspin, time_steps, dt, kick, order, ot, observable)
865 type(namespace_t), intent(in) :: namespace
866 integer, intent(in) :: nspin, time_steps
867 real(real64), intent(in) :: dt
868 type(kick_t), intent(in) :: kick
869 integer, intent(in) :: order
870 real(real64), intent(in) :: ot(0:time_steps)
871 integer, intent(in) :: observable(2)
872
873 integer :: iunit, i
874 character(len=20) :: header_string
875 type(unit_t) :: ot_unit
876
877 push_sub(write_ot)
878
879 iunit = io_open('ot', namespace, action='write', status='replace')
880
881 write(iunit, '(a15,i2)') '# nspin ', nspin
882 write(iunit, '(a15,i2)') '# Order ', order
883 write(iunit, '(a28,i2)') '# Frequencies subtracted = ', 0
884 select case (observable(1))
885 case (-1)
886 write(iunit,'(a)') '# Observable operator = kick operator'
887 if (kick%n_multipoles > 0) then
888 ot_unit = units_out%length**kick%l(1)
889 else
890 ot_unit = units_out%length
891 end if
892 case (0)
893 select case (observable(2))
894 case (1)
895 write(iunit,'(a)') '# O = x'
896 case (2)
897 write(iunit,'(a)') '# O = y'
898 case (3)
899 write(iunit,'(a)') '# O = z'
900 end select
901 ot_unit = units_out%length
902 case default
903 ot_unit = units_out%length**observable(1)
904 write(iunit, '(a12,i1,a1,i2,a1)') '# (l, m) = (', observable(1),',',observable(2),')'
905 end select
906 call kick_write(kick, iunit)
907
908 ! Units
909 write(iunit,'(a1)', advance = 'no') '#'
910 write(iunit,'(a20)', advance = 'no') str_center('t', 19)
911 write(iunit,'(a20)', advance = 'yes') str_center('<O>(t)', 20)
912 write(iunit,'(a1)', advance = 'no') '#'
913 write(header_string, '(a)') '['//trim(units_abbrev(units_out%time))//']'
914 write(iunit,'(a20)', advance = 'no') str_center(trim(header_string), 19)
915 write(header_string, '(a)') '['//trim(units_abbrev(units_out%length))//']'
916 write(iunit,'(a20)', advance = 'yes') str_center(trim(header_string), 20)
917
918 do i = 0, time_steps
919 write(iunit, '(2e20.8)') units_from_atomic(units_out%time, i*dt), units_from_atomic(ot_unit, ot(i))
920 end do
921
922 call io_close(iunit)
923 pop_sub(write_ot)
924 end subroutine write_ot
925
926
927 ! ---------------------------------------------------------
928 subroutine read_ot(namespace, nspin, order, nw_subtracted)
929 type(namespace_t), intent(in) :: namespace
930 integer, intent(out) :: nspin
931 integer, intent(out) :: order
932 integer, intent(out) :: nw_subtracted
933
934 integer :: iunit, i, ierr
935 character(len=100) :: line
936 character(len=12) :: dummychar
937 real(real64) :: dummy, t1, t2
938 type(unit_t) :: ot_unit
939
940 push_sub(read_ot)
941
942 iunit = io_open('ot', namespace, action='read', status='old')
943
944 read(iunit, '(15x,i2)') nspin
945 read(iunit, '(15x,i2)') order
946 read(iunit, '(28x,i2)') nw_subtracted
947 read(iunit, '(a)') line
948
949 i = index(line, 'Observable')
950 if (index(line, 'Observable') /= 0) then
951 observable(1) = -1
952 elseif (index(line, '# O =') /= 0) then
953 observable(1) = 0
954 if (index(line,'x') /= 0) then
955 observable(2) = 1
956 elseif (index(line,'y') /= 0) then
957 observable(2) = 2
958 elseif (index(line,'z') /= 0) then
959 observable(2) = 3
960 end if
961 elseif (index(line, '# (l, m) = ') /= 0) then
962 read(line,'(a12,i1,a1,i2,a1)') dummychar(1:12), observable(1), dummychar(1:1), observable(2), dummychar(1:1)
963 else
964 write(message(1),'(a)') 'Problem reading "ot" file: could not figure out the shape'
965 write(message(2),'(a)') 'of the observation operator.'
966 call messages_fatal(2, namespace=namespace)
967 end if
968
969 call kick_read(kick, iunit, namespace)
970 read(iunit, '(a)') line
971 read(iunit, '(a)') line
972 call io_skip_header(iunit)
973
974 ! Figure out about the units of the file
975 call unit_system_from_file(units, "ot", namespace, ierr)
976 if (ierr /= 0) then
977 write(message(1), '(a)') 'Could not figure out the units in file "ot".'
978 call messages_fatal(1, namespace=namespace)
979 end if
980
981 select case (observable(1))
982 case (-1)
983 if (kick%n_multipoles > 0) then
984 ot_unit = units_out%length**kick%l(1)
985 else
986 ot_unit = units_out%length
987 end if
988 case (0)
989 ot_unit = units_out%length
990 case default
991 ot_unit = units_out%length**observable(1)
992 end select
993
994 ! count number of time_steps
995 time_steps = 0
996 do
997 read(iunit, *, end=100) dummy
998 time_steps = time_steps + 1
999 if (time_steps == 1) t1 = dummy
1000 if (time_steps == 2) t2 = dummy
1001 end do
1002100 continue
1003 dt = units_to_atomic(units%time, (t2 - t1)) ! units_out is OK
1004
1005 call io_skip_header(iunit)
1006
1007 safe_allocate(ot(0:time_steps))
1008
1009 do i = 0, time_steps-1
1010 read(iunit, *) dummy, ot(i)
1011 ot(i) = units_to_atomic(ot_unit, ot(i))
1012 end do
1013
1014 pop_sub(read_ot)
1015 end subroutine read_ot
1016
1017
1018 ! ---------------------------------------------------------
1019 subroutine print_omega_file(namespace, omega, search_interval, final_time, nfrequencies)
1020 type(namespace_t), intent(in) :: namespace
1021 real(real64), intent(inout) :: omega
1022 real(real64), intent(inout) :: search_interval
1023 real(real64), intent(inout) :: final_time
1024 integer, intent(inout) :: nfrequencies
1025
1026 integer :: iunit, i, ierr, nspin, order, nw_subtracted
1027 logical :: file_exists
1028 character(len=20) :: header_string
1029 real(real64), allocatable :: warray(:), tarray(:)
1030 real(real64) :: leftbound, rightbound, dw, w, aw
1031
1032 push_sub(print_omega_file)
1033
1034 ! First, let us check that the file "ot" exists.
1035 inquire(file="ot", exist = file_exists)
1036 if (.not. file_exists) then
1037 write(message(1),'(a)') "Could not find 'ot' file."
1038 call messages_fatal(1, namespace=namespace)
1039 end if
1040
1041 ! Now, we should find out which units the file "ot" has.
1042 call unit_system_from_file(units, "ot", namespace, ierr)
1043 if (ierr /= 0) then
1044 write(message(1),'(a)') "Could not retrieve units in the 'ot' file."
1045 call messages_fatal(1, namespace=namespace)
1046 end if
1047
1048 if (omega > m_zero) then
1049 omega = units_to_atomic(units%energy, omega)
1050 else
1051 omega = m_half
1052 end if
1053
1054 if (search_interval > m_zero) then
1055 search_interval = units_to_atomic(units%energy, search_interval)
1056 else
1057 search_interval = m_half
1058 end if
1059
1060 leftbound = omega - search_interval
1061 rightbound = omega + search_interval
1062
1063 safe_allocate(warray(1:nfrequencies))
1064 safe_allocate(tarray(1:nfrequencies))
1065
1066 call read_ot(namespace, nspin, order, nw_subtracted)
1067
1068 if (mod(order, 2) == 1) then
1070 else
1072 end if
1073
1074 if (final_time > m_zero) then
1075 total_time = units_to_atomic(units%time, final_time)
1076 if (total_time > dt*time_steps) then
1078 write(message(1), '(a)') 'The requested total time to process is larger than the time available in the input file.'
1079 write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', &
1080 units_from_atomic(units_out%time, total_time), trim(units_abbrev(units_out%time))
1081 call messages_warning(2, namespace=namespace)
1082 end if
1083 time_steps = int(total_time / dt)
1085 else
1087 end if
1088
1089 warray = m_zero
1090 tarray = m_zero
1091 dw = (rightbound-leftbound) / (nfrequencies - 1)
1092 do i = 1, nfrequencies
1093 w = leftbound + (i-1)*dw
1094 warray(i) = w
1095 call ft(w, aw)
1096 tarray(i) = aw
1097 end do
1098
1099 iunit = io_open('omega', namespace, action='write', status='replace')
1100 write(iunit, '(a15,i2)') '# nspin ', nspin
1101 call kick_write(kick, iunit)
1102 write(iunit, '(a)') '#%'
1103 write(iunit, '(a1,a20)', advance = 'no') '#', str_center("omega", 20)
1104 write(header_string,'(a)') 'F(omega)'
1105 write(iunit, '(a20)', advance = 'yes') str_center(trim(header_string), 20)
1106 write(iunit, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20)
1107 ! Here we should print the units of the transform.
1108 write(iunit, '(a)', advance = 'yes')
1109
1110 do i = 1, nfrequencies
1111 write(iunit,'(2e20.8)') units_from_atomic(units_out%energy, warray(i)), &
1112 tarray(i)
1113 end do
1115 safe_deallocate_a(warray)
1116 safe_deallocate_a(tarray)
1117
1118 pop_sub(print_omega_file)
1119 end subroutine print_omega_file
1120 ! ---------------------------------------------------------
1121
1123
1124! ---------------------------------------------------------
1125! ---------------------------------------------------------
1126! ---------------------------------------------------------
1127program oscillator_strength
1129 use global_oct_m
1130 use io_oct_m
1131 use messages_oct_m
1133
1134 implicit none
1135
1136 integer :: run_mode, order, nfrequencies, ierr, nresonances
1137 real(real64) :: omega, search_interval, final_time, damping
1138 integer, parameter :: &
1139 ANALYZE_NTHORDER_SIGNAL = 1, &
1140 generate_nthorder_signal = 2, &
1141 read_resonances_from_file = 3, &
1142 generate_omega_file = 4
1143 character(len=100) :: ffile
1144 character(kind=c_char) :: cfile(101)
1145
1146 ! Reads the information passed through the command line options (if available).
1147 call getopt_init(ierr)
1148 if (ierr /= 0) then
1149 message(1) = "Your Fortran compiler doesn't support command-line arguments;"
1150 message(2) = "the oct-oscillator-strength command is not available."
1151 call messages_fatal(2)
1152 end if
1153
1154 ! Set the default values
1155 run_mode = analyze_nthorder_signal
1156 omega = - m_one
1157 search_interval = - m_one
1158 order = 1
1159 nfrequencies = 1000
1160 final_time = - m_one
1161 nresonances = 1
1162 observable(1) = -1
1163 observable(2) = 0
1164 ffile = ""
1165 damping = 0.1_real64/27.2114_real64 ! This is the usual damping factor used in the literature.
1166
1167 ! Get the parameters from the command line.
1168 cfile = c_null_char
1169 call getopt_oscillator_strength(run_mode, omega, search_interval, &
1170 order, nresonances, nfrequencies, final_time, &
1171 observable(1), observable(2), damping, cfile)
1172 call string_c_to_f(cfile, ffile)
1173 call getopt_end()
1174
1175 ! Initialize stuff
1177 call parser_init()
1178 call io_init(defaults = .true.)
1179
1180 select case (run_mode)
1181 case (generate_nthorder_signal)
1182 call generate_signal(global_namespace, order, observable)
1183 case (analyze_nthorder_signal)
1184 call analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, damping, &
1185 global_namespace)
1186 case (read_resonances_from_file)
1187 call read_resonances_file(order, ffile, global_namespace, search_interval, final_time, nfrequencies)
1188 case (generate_omega_file)
1189 call print_omega_file(global_namespace, omega, search_interval, final_time, nfrequencies)
1190 case default
1191 end select
1192
1193 call io_end()
1194 call parser_end()
1195 call global_end()
1196
1197end program oscillator_strength
1198! ---------------------------------------------------------
1199
1200!! Local Variables:
1201!! mode: f90
1202!! coding: utf-8
1203!! End:
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
subroutine, public getopt_init(ierr)
Initializes the getopt machinery. Must be called before attempting to parse the options....
subroutine, public getopt_end
real(real64), parameter, public m_two
Definition: global.F90:193
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
Definition: global.F90:421
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
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:280
subroutine, public init_octopus_globals(comm)
Initialise Octopus-specific global constants and files. This routine performs no initialisation calls...
Definition: global.F90:390
complex(real64), parameter, public m_z0
Definition: global.F90:201
complex(real64), parameter, public m_zi
Definition: global.F90:205
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_three
Definition: global.F90:194
Definition: io.F90:116
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:165
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
subroutine, public io_skip_header(iunit)
Definition: io.F90:645
subroutine, public io_end()
Definition: io.F90:270
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
subroutine, public kick_read(kick, iunit, namespace)
Definition: kick.F90:820
subroutine, public kick_write(kick, iunit, out)
Definition: kick.F90:892
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
subroutine ft(omega, power)
integer, dimension(2) observable
subroutine write_ot(namespace, nspin, time_steps, dt, kick, order, ot, observable)
subroutine generate_signal(namespace, order, observable)
subroutine ft2(omega, power)
subroutine read_ot(namespace, nspin, order, nw_subtracted)
subroutine resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound)
integer, parameter cosine_transform
subroutine resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, c01, c02)
subroutine analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, damping, namespace)
subroutine modify_ot(time_steps, dt, order, ot, omega, power)
subroutine print_omega_file(namespace, omega, search_interval, final_time, nfrequencies)
subroutine find_resonance(omega, leftbound, rightbound, nfrequencies)
subroutine read_resonances_file(order, ffile, namespace, search_interval, final_time, nfrequencies)
subroutine write_polarizability(namespace, nfrequencies, nresonances, dw, w, c0I2, gamma)
Implements the SOS formula of the polarizability, and writes down to the "polarizability" file the re...
subroutine, public spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax)
Definition: spectrum.F90:2341
character(len=80) function, public str_center(s_in, l_in)
puts space around string, so that it is centered
Definition: string.F90:199
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
subroutine, public unit_system_from_file(uu, fname, namespace, ierr)
This is a very primitive procedure that attempts to find out which units were used to write an octopu...
program oscillator_strength
int true(void)