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 21module oscillator_strength_oct_m 22 use global_oct_m 23 use io_oct_m 24 use kick_oct_m 25 use lalg_adv_oct_m 26 use messages_oct_m 27 use minimizer_oct_m 28 use namespace_oct_m 29 use parser_oct_m 30 use profiling_oct_m 31 use spectrum_oct_m 32 use string_oct_m 33 use unit_oct_m 34 use unit_system_oct_m 35 36 implicit none 37 38 integer, parameter :: SINE_TRANSFORM = 1, & 39 COSINE_TRANSFORM = 2 40 41 type local_operator_t 42 integer :: n_multipoles 43 integer, pointer :: l(:), m(:) 44 FLOAT, pointer :: weight(:) 45 end type local_operator_t 46 47 integer :: observable(2) 48 type(unit_system_t) :: units 49 FLOAT, allocatable :: ot(:) 50 type(kick_t) :: kick 51 integer :: time_steps 52 FLOAT :: total_time 53 integer :: mode 54 FLOAT :: dt 55 56contains 57 58 ! --------------------------------------------------------- 59 subroutine ft(omega, power) 60 FLOAT, intent(in) :: omega 61 FLOAT, intent(out) :: power 62 63 integer :: j 64 FLOAT :: x 65 power = M_ZERO 66 67 PUSH_SUB(ft) 68 69 select case(mode) 70 71 case(SINE_TRANSFORM) 72 73 do j = 0, time_steps 74 x = sin(omega*j*dt) 75 power = power + x*ot(j) 76 end do 77 power = power*dt / (dt*time_steps) 78 79 case(COSINE_TRANSFORM) 80 81 do j = 0, time_steps 82 x = cos(omega*j*dt) 83 power = power + x*ot(j) 84 end do 85 power = power*dt / (dt*time_steps) 86 87 end select 88 89 POP_SUB(ft) 90 end subroutine ft 91 92 93 ! --------------------------------------------------------- 94 subroutine ft2(omega, power) 95 FLOAT, intent(in) :: omega 96 FLOAT, intent(out) :: power 97 98 integer :: j 99 FLOAT :: x 100 power = M_ZERO 101 102 PUSH_SUB(ft2) 103 104 select case(mode) 105 106 case(SINE_TRANSFORM) 107 108 do j = 0, time_steps 109 x = sin(omega*j*dt) 110 power = power + x*ot(j) 111 end do 112 ! The function should be *minus* the sine Fourier transform, since this 113 ! is the function to be minimized. 114 power = - (power*dt / (dt*time_steps))**2 115 116 case(COSINE_TRANSFORM) 117 118 do j = 0, time_steps 119 x = cos(omega*j*dt) 120 power = power + x*ot(j) 121 end do 122 power = - (power*dt / (dt*time_steps))**2 123 124 end select 125 126 POP_SUB(ft2) 127 end subroutine ft2 128 129 130 ! --------------------------------------------------------- 131 subroutine local_operator_copy(o, i) 132 type(local_operator_t), intent(inout) :: o 133 type(local_operator_t), intent(inout) :: i 134 135 integer :: j 136 137 PUSH_SUB(local_operator_copy) 138 139 o%n_multipoles = i%n_multipoles 140 SAFE_ALLOCATE( o%l(1:o%n_multipoles)) 141 SAFE_ALLOCATE( o%m(1:o%n_multipoles)) 142 SAFE_ALLOCATE(o%weight(1:o%n_multipoles)) 143 144 do j = 1, o%n_multipoles 145 o%l(j) = i%l(j) 146 o%m(j) = i%m(j) 147 o%weight(j) = i%weight(j) 148 end do 149 150 POP_SUB(local_operator_copy) 151 end subroutine local_operator_copy 152 153 ! --------------------------------------------------------- 154 subroutine read_resonances_file(order, ffile, namespace, search_interval, final_time, nfrequencies) 155 integer, intent(inout) :: order 156 character(len=*), intent(in) :: ffile 157 type(namespace_t), intent(in):: namespace 158 FLOAT, intent(inout) :: search_interval 159 FLOAT, intent(in) :: final_time 160 integer, intent(in) :: nfrequencies 161 162 FLOAT :: dummy, leftbound, rightbound, w, power, dw 163 integer :: iunit, nresonances, ios, i, j, k, npairs, nspin, order_in_file, nw_subtracted, ierr 164 logical :: file_exists 165 FLOAT, allocatable :: wij(:), omega(:), c0i(:) 166 167 PUSH_SUB(read_resonances_file) 168 169 if(order /= 2) then 170 write(message(1),'(a)') 'The run mode #3 is only compatible with the analysis of the' 171 write(message(2),'(a)') 'second-order response.' 172 call messages_fatal(2) 173 end if 174 175 ! First, let us check that the file "ot" exists. 176 inquire(file="ot", exist = file_exists) 177 if(.not.file_exists) then 178 write(message(1),'(a)') "Could not find 'ot' file." 179 call messages_fatal(1) 180 end if 181 182 ! Now, we should find out which units the file "ot" has. 183 call unit_system_from_file(units, "ot", namespace, ierr) 184 if(ierr /= 0) then 185 write(message(1),'(a)') "Could not retrieve units in the 'ot' file." 186 call messages_fatal(1) 187 end if 188 189 mode = COSINE_TRANSFORM 190 191 iunit = io_open(trim(ffile), namespace, action='read', status='old', die=.false.) 192 if(iunit == 0) then 193 write(message(1),'(a)') 'Could not open '//trim(ffile)//' file.' 194 call messages_fatal(1) 195 end if 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 230 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) 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 FLOAT, intent(inout) :: omega 290 FLOAT, intent(inout) :: search_interval 291 FLOAT, intent(inout) :: final_time 292 integer, intent(inout) :: nresonances 293 integer, intent(inout) :: nfrequencies 294 FLOAT, intent(in) :: damping 295 type(namespace_t), intent(in) :: namespace 296 297 FLOAT :: leftbound, rightbound, dw, power 298 FLOAT, 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) 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) 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 ! --------------------------------------------------------- 409 !> Implements the SOS formula of the polarizability, and writes 410 !! down to the "polarizability" file the real and imaginary part 411 !! of the dynamical polarizability. 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 FLOAT, intent(in) :: dw 416 FLOAT, intent(in) :: w(nresonances), c0I2(nresonances) 417 FLOAT, intent(in) :: gamma 418 419 integer :: iunit, i, j 420 FLOAT :: e 421 CMPLX :: 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 FLOAT, intent(inout) :: omega 456 FLOAT, intent(in) :: leftbound, rightbound 457 integer, intent(in) :: nfrequencies 458 459 integer :: i, ierr 460 FLOAT :: dw, w, aw, min_aw, min_w, omega_orig, min_w1, min_w2 461 FLOAT, 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; tarray = M_ZERO 469 dw = (rightbound-leftbound) / (nfrequencies - 1) 470 do i = 1, nfrequencies 471 w = leftbound + (i-1)*dw 472 warray(i) = w 473 call ft2(w, aw) 474 tarray(i) = aw 475 end do 476 477 min_w = omega 478 min_aw = M_ZERO 479 do i = 1, nfrequencies 480 w = leftbound + (i-1)*dw 481 if(tarray(i) < min_aw) then 482 min_aw = tarray(i) 483 min_w = w 484 end if 485 end do 486 487 omega_orig = omega 488 omega = min_w 489 min_w1 = min_w - 2*dw 490 min_w2 = min_w + 2*dw 491 call loct_1dminimize(min_w1, min_w2, omega, ft2, ierr) 492 493 if(ierr /= 0) then 494 write(message(1),'(a)') 'Could not find a maximum.' 495 write(message(2),'(a)') 496 write(message(3), '(a,f12.8,a,f12.8,a)') ' Search interval = [', & 497 units_from_atomic(units%energy, leftbound), ',', units_from_atomic(units%energy, rightbound), ']' 498 write(message(4), '(a,f12.4,a)') ' Search discretization = ', & 499 units_from_atomic(units%energy, dw), ' '//trim(units_abbrev(units%energy)) 500 call messages_fatal(4) 501 end if 502 503 SAFE_DEALLOCATE_A(warray) 504 SAFE_DEALLOCATE_A(tarray) 505 506 POP_SUB(find_resonance) 507 end subroutine find_resonance 508 ! --------------------------------------------------------- 509 510 511 ! --------------------------------------------------------- 512 subroutine resonance_first_order(omega, power, nw_subtracted, dw, leftbound, rightbound) 513 FLOAT, intent(in) :: omega 514 FLOAT, intent(out) :: power 515 integer, intent(in) :: nw_subtracted 516 FLOAT, intent(in) :: dw, leftbound, rightbound 517 518 PUSH_SUB(resonance_first_order) 519 520 call ft(omega, power) 521 522 select case(mode) 523 case(SINE_TRANSFORM) 524 power = power / (M_ONE - sin(M_TWO*omega*total_time)/(M_TWO*omega*total_time)) 525 case(COSINE_TRANSFORM) 526 call messages_not_implemented("resonance first order cosine transform") 527 end select 528 529 write(message(1), '(a)') '******************************************************************' 530 write(message(2), '(a,i3)') 'Resonance #', nw_subtracted + 1 531 write(message(3), '(a,f12.8,a,f12.8,a)') 'omega = ', units_from_atomic(units_out%energy, omega), & 532 ' '//trim(units_abbrev(units_out%energy))//' = ', omega, ' Ha' 533 write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**2, power), & 534 ' '//trim(units_abbrev(units_out%length**2))//' =', power, ' b^2' 535 write(message(5), '(a,f12.8,a,f12.8,a)') '<0|P|I> = ', units_from_atomic(units_out%length, sqrt(abs(power))), & 536 ' '//trim(units_abbrev(units_out%length))//' = ', sqrt(abs(power)),' b' 537 write(message(6), '(a,f12.8)') 'f[O->I] = ', M_TWO*omega*power 538 write(message(7), '(a)') 539 write(message(8), '(a,f12.8,a,f12.8,a)') ' Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', & 540 units_from_atomic(units_out%energy, rightbound), ']' 541 write(message(9), '(a,f12.4,a)') ' Search discretization = ', units_from_atomic(units_out%energy, dw), & 542 ' '//trim(units_abbrev(units_out%energy)) 543 write(message(10), '(a)') '******************************************************************' 544 write(message(11), '(a)') 545 call messages_info(11) 546 547 POP_SUB(resonance_first_order) 548 549 end subroutine resonance_first_order 550 ! --------------------------------------------------------- 551 552 553 ! --------------------------------------------------------- 554 subroutine resonance_second_order(omega, power, nw_subtracted, leftbound, rightbound, c01, c02) 555 FLOAT, intent(in) :: omega 556 FLOAT, intent(out) :: power 557 integer, intent(in) :: nw_subtracted 558 FLOAT, intent(in) :: leftbound, rightbound 559 FLOAT, intent(in) :: c01, c02 560 561 PUSH_SUB(resonance_second_order) 562 563 call ft(omega, power) 564 select case(mode) 565 case(SINE_TRANSFORM) 566 power = power / (M_ONE - sin(M_TWO*omega*total_time)/(M_TWO*omega*total_time)) 567 case(COSINE_TRANSFORM) 568 ! WARNING: there is some difference between the omega=0 case and the rest. 569 if(omega /= M_ZERO) then 570 power = power / (M_ONE + sin(M_TWO*omega*total_time)/(M_TWO*omega*total_time)) 571 else 572 power = power / M_TWO 573 end if 574 end select 575 576 write(message(1), '(a)') '******************************************************************' 577 write(message(2), '(a,i3)') 'Resonance #', nw_subtracted + 1 578 write(message(3), '(a,f12.8,a,f12.8,a)') 'omega = ', units_from_atomic(units_out%energy, omega), & 579 ' '//trim(units_abbrev(units_out%energy))//' = ', omega, ' Ha' 580 write(message(4), '(a,f12.8,a,f12.8,a)') 'C(omega) = ', units_from_atomic(units_out%length**3, power), & 581 ' '//trim(units_abbrev(units_out%length**3))//' = ', power, ' b^3' 582 call messages_info(4) 583 584 if(c01*c02 /= M_ZERO) then 585 write(message(1), '(a,f12.8)') ' C(omega)/(C0i*C0j) = ', power / (c01 * c02) 586 call messages_info(1) 587 end if 588 589 write(message(1), '(a)') 590 write(message(2), '(a,f12.8,a,f12.8,a)') ' Search interval = [', units_from_atomic(units_out%energy, leftbound), ',', & 591 units_from_atomic(units_out%energy, rightbound), ']' 592 write(message(3), '(a)') '******************************************************************' 593 write(message(4), '(a)') 594 call messages_info(4) 595 596 POP_SUB(resonance_second_order) 597 end subroutine resonance_second_order 598 ! --------------------------------------------------------- 599 600 601 ! --------------------------------------------------------- 602 subroutine generate_signal(namespace, order, observable) 603 type(namespace_t), intent(in) :: namespace 604 integer, intent(in) :: order 605 integer, intent(in) :: observable(2) 606 607 logical :: file_exists 608 integer :: i, j, nspin, time_steps, lmax, nfiles, k, add_lm, l, m, max_add_lm 609 integer, allocatable :: iunit(:) 610 FLOAT :: dt, lambda, dump, o0 611 type(unit_t) :: mp_unit 612 FLOAT, allocatable :: q(:), mu(:), qq(:, :), c(:) 613 character(len=20) :: filename 614 type(kick_t) :: kick 615 type(unit_system_t) :: units 616 FLOAT, allocatable :: multipole(:, :, :), ot(:), dipole(:, :) 617 type(local_operator_t) :: kick_operator 618 type(local_operator_t) :: obs 619 620 PUSH_SUB(generate_signal) 621 622 ! Find out how many files do we have 623 nfiles = 0 624 do 625 write(filename,'(a11,i1)') 'multipoles.', nfiles+1 626 inquire(file=trim(filename), exist = file_exists) 627 if(.not.file_exists) exit 628 nfiles = nfiles + 1 629 end do 630 631 ! WARNING: Check that order is smaller or equal to nfiles 632 if(nfiles == 0) then 633 write(message(1),'(a)') 'No multipoles.x file was found' 634 call messages_fatal(1) 635 end if 636 if(order > nfiles) then 637 write(message(1),'(a)') 'The order that you ask for is higher than the number' 638 write(message(2),'(a)') 'of multipoles.x file that you supply.' 639 call messages_fatal(2) 640 end if 641 642 ! Open the files. 643 SAFE_ALLOCATE(iunit(1:nfiles)) 644 do j = 1, nfiles 645 write(filename,'(a11,i1)') 'multipoles.', j 646 iunit(j) = io_open(trim(filename), namespace, action='read', status='old', die=.false.) 647 end do 648 649 SAFE_ALLOCATE( q(1:nfiles)) 650 SAFE_ALLOCATE(mu(1:nfiles)) 651 SAFE_ALLOCATE(qq(1:nfiles, 1:nfiles)) 652 SAFE_ALLOCATE( c(1:nfiles)) 653 654 c = M_ZERO 655 c(order) = M_ONE 656 657 ! Get the basic info from the first file 658 call spectrum_mult_info(namespace, iunit(1), nspin, kick, time_steps, dt, units, lmax=lmax) 659 660 ! Sets the kick operator... 661 if(kick%n_multipoles > 0) then 662 kick_operator%n_multipoles = kick%n_multipoles 663 SAFE_ALLOCATE( kick_operator%l(1:kick_operator%n_multipoles)) 664 SAFE_ALLOCATE( kick_operator%m(1:kick_operator%n_multipoles)) 665 SAFE_ALLOCATE(kick_operator%weight(1:kick_operator%n_multipoles)) 666 do i = 1, kick_operator%n_multipoles 667 kick_operator%l(i) = kick%l(i) 668 kick_operator%m(i) = kick%m(i) 669 kick_operator%weight(i) = kick%weight(i) 670 end do 671 else 672 kick_operator%n_multipoles = 3 673 SAFE_ALLOCATE( kick_operator%l(1:kick_operator%n_multipoles)) 674 SAFE_ALLOCATE( kick_operator%m(1:kick_operator%n_multipoles)) 675 SAFE_ALLOCATE(kick_operator%weight(1:kick_operator%n_multipoles)) 676 kick_operator%l(1:3) = 1 677 kick_operator%m(1) = -1 678 kick_operator%m(2) = 0 679 kick_operator%m(3) = 1 680 ! WARNING: not sure if m = -1 => y, and m = 1 => x. What is the convention? 681 kick_operator%weight(1) = -sqrt((M_FOUR*M_PI)/M_THREE) * kick%pol(2, kick%pol_dir) 682 kick_operator%weight(2) = sqrt((M_FOUR*M_PI)/M_THREE) * kick%pol(3, kick%pol_dir) 683 kick_operator%weight(3) = -sqrt((M_FOUR*M_PI)/M_THREE) * kick%pol(1, kick%pol_dir) 684 end if 685 686 ! Sets the observation operator 687 select case(observable(1)) 688 case(-1) 689 ! This means that the "observation operator" should be equal 690 ! to the "perturbation operator", i.e., the kick. 691 call local_operator_copy(obs, kick_operator) 692 case(0) 693 ! This means that the observable is the dipole operator; observable(2) determines 694 ! if it is x, y or z. 695 obs%n_multipoles = 1 696 SAFE_ALLOCATE(obs%l(1:1)) 697 SAFE_ALLOCATE(obs%m(1:1)) 698 SAFE_ALLOCATE(obs%weight(1:1)) 699 obs%l(1) = 1 700 select case(observable(2)) 701 case(1) 702 obs%m(1) = -1 703 obs%weight(1) = -sqrt((M_FOUR*M_PI)/M_THREE) 704 case(2) 705 obs%m(1) = 1 706 obs%weight(1) = sqrt((M_FOUR*M_PI)/M_THREE) 707 case(3) 708 obs%m(1) = 0 709 obs%weight(1) = -sqrt((M_FOUR*M_PI)/M_THREE) 710 end select 711 case default 712 ! This means that the observation operator is (l,m) = (observable(1), observable(2)) 713 obs%n_multipoles = 1 714 SAFE_ALLOCATE(obs%l(1:1)) 715 SAFE_ALLOCATE(obs%m(1:1)) 716 SAFE_ALLOCATE(obs%weight(1:1)) 717 obs%weight(1) = M_ONE 718 obs%l(1) = observable(1) 719 obs%m(1) = observable(2) 720 end select 721 722 lambda = abs(kick%delta_strength) 723 q(1) = kick%delta_strength / lambda 724 725 do j = 2, nfiles 726 call spectrum_mult_info(namespace, iunit(j), nspin, kick, time_steps, dt, units, lmax=lmax) 727 q(j) = kick%delta_strength / lambda 728 end do 729 730 do i = 1, nfiles 731 do j = 1, nfiles 732 qq(i, j) = q(j)**i 733 end do 734 end do 735 736 call lalg_inverter(nfiles, qq) 737 738 mu = matmul(qq, c) 739 740 if(kick%n_multipoles > 0) then 741 lmax = maxval(kick%l(1:obs%n_multipoles)) 742 max_add_lm = (lmax+1)**2-1 743 SAFE_ALLOCATE(multipole(1:max_add_lm, 0:time_steps, 1:nspin)) 744 ! The units have nothing to do with the perturbing kick?? 745 mp_unit = units%length**kick%l(1) 746 else 747 max_add_lm = 3 748 SAFE_ALLOCATE(multipole(1:3, 0:time_steps, 1:nspin)) 749 mp_unit = units%length 750 end if 751 SAFE_ALLOCATE(ot(0:time_steps)) 752 multipole = M_ZERO 753 ot = M_ZERO 754 755 SAFE_ALLOCATE(dipole(1:3, 1:nspin)) 756 757 do j = 1, nfiles 758 call io_skip_header(iunit(j)) 759 760 do i = 0, time_steps 761 select case(nspin) 762 case(1) 763 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm) 764 case(2) 765 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), & 766 dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm) 767 case(4) 768 read(iunit(j), *) k, dump, dump, (multipole(add_lm, i, 1), add_lm = 1, max_add_lm), & 769 dump, (multipole(add_lm, i, 2), add_lm = 1, max_add_lm), & 770 dump, (multipole(add_lm, i, 3), add_lm = 1, max_add_lm), & 771 dump, (multipole(add_lm, i, 4), add_lm = 1, max_add_lm) 772 end select 773 multipole(1:max_add_lm, i, :) = units_to_atomic(mp_unit, multipole(1:max_add_lm, i, :)) 774 775 ! The dipole is treated differently in the multipoles file: first of all, 776 ! the program should have printed *minus* the dipole operator. 777 dipole(1:3, 1:nspin) = - multipole(1:3, i, 1:nspin) 778 ! And then it contains the "cartesian" dipole, opposed to the spherical dipole: 779 multipole(1, i, 1:nspin) = -sqrt(M_THREE/(M_FOUR*M_PI)) * dipole(2, 1:nspin) 780 multipole(2, i, 1:nspin) = sqrt(M_THREE/(M_FOUR*M_PI)) * dipole(3, 1:nspin) 781 multipole(3, i, 1:nspin) = -sqrt(M_THREE/(M_FOUR*M_PI)) * dipole(1, 1:nspin) 782 783 dump = M_ZERO 784 do k = 1, obs%n_multipoles 785 add_lm = 1; l = 1 786 lcycle: do 787 do m = -l, l 788 if(l == obs%l(k) .and. m == obs%m(k)) exit lcycle 789 add_lm = add_lm + 1 790 end do 791 l = l + 1 792 end do lcycle 793 ! Warning: it should not add the nspin components? 794 dump = dump + obs%weight(k) * sum(multipole(add_lm, i, 1:nspin)) 795 end do 796 797 if(i == 0) o0 = dump 798 ot(i) = ot(i) + mu(j)*(dump - o0) 799 end do 800 801 end do 802 803 ot = ot / lambda**order 804 805 call write_ot(namespace, nspin, time_steps, dt, kick, order, ot(0:time_steps), observable) 806 807 ! Close files and exit. 808 do j = 1, nfiles 809 call io_close(iunit(j)) 810 end do 811 SAFE_DEALLOCATE_A(iunit) 812 SAFE_DEALLOCATE_A(q) 813 SAFE_DEALLOCATE_A(mu) 814 SAFE_DEALLOCATE_A(qq) 815 SAFE_DEALLOCATE_A(c) 816 SAFE_DEALLOCATE_A(ot) 817 SAFE_DEALLOCATE_A(multipole) 818 SAFE_DEALLOCATE_A(dipole) 819 820 POP_SUB(generate_signal) 821 end subroutine generate_signal 822 ! --------------------------------------------------------- 823 824 825 ! --------------------------------------------------------- 826 subroutine modify_ot(time_steps, dt, order, ot, omega, power) 827 integer, intent(in) :: time_steps 828 FLOAT, intent(in) :: dt 829 integer, intent(in) :: order 830 FLOAT, intent(inout) :: ot(0:time_steps) 831 FLOAT, intent(in) :: omega 832 FLOAT, intent(in) :: power 833 834 integer :: i 835 836 PUSH_SUB(modify_ot) 837 838 select case(mod(order, 2)) 839 case(1) 840 do i = 0, time_steps 841 ot(i) = ot(i) - M_TWO * power * sin(omega*i*dt) 842 end do 843 case(0) 844 if(abs(omega) <= M_EPSILON) then 845 do i = 0, time_steps 846 ot(i) = ot(i) - power * cos(omega*i*dt) 847 end do 848 else 849 do i = 0, time_steps 850 ot(i) = ot(i) - M_TWO * power * cos(omega*i*dt) 851 end do 852 end if 853 end select 854 855 POP_SUB(modify_ot) 856 end subroutine modify_ot 857 ! --------------------------------------------------------- 858 859 860 ! --------------------------------------------------------- 861 subroutine write_ot(namespace, nspin, time_steps, dt, kick, order, ot, observable) 862 type(namespace_t), intent(in) :: namespace 863 integer, intent(in) :: nspin, time_steps 864 FLOAT, intent(in) :: dt 865 type(kick_t), intent(in) :: kick 866 integer, intent(in) :: order 867 FLOAT, intent(in) :: ot(0:time_steps) 868 integer, intent(in) :: observable(2) 869 870 integer :: iunit, i 871 character(len=20) :: header_string 872 type(unit_t) :: ot_unit 873 874 PUSH_SUB(write_ot) 875 876 iunit = io_open('ot', namespace, action='write', status='replace') 877 878 write(iunit, '(a15,i2)') '# nspin ', nspin 879 write(iunit, '(a15,i2)') '# Order ', order 880 write(iunit, '(a28,i2)') '# Frequencies subtracted = ', 0 881 select case(observable(1)) 882 case(-1) 883 write(iunit,'(a)') '# Observable operator = kick operator' 884 if(kick%n_multipoles > 0 ) then 885 ot_unit = units_out%length**kick%l(1) 886 else 887 ot_unit = units_out%length 888 end if 889 case(0) 890 select case(observable(2)) 891 case(1); write(iunit,'(a)') '# O = x' 892 case(2); write(iunit,'(a)') '# O = y' 893 case(3); write(iunit,'(a)') '# O = z' 894 end select 895 ot_unit = units_out%length 896 case default 897 ot_unit = units_out%length**observable(1) 898 write(iunit, '(a12,i1,a1,i2,a1)') '# (l, m) = (', observable(1),',',observable(2),')' 899 end select 900 call kick_write(kick, iunit) 901 902 ! Units 903 write(iunit,'(a1)', advance = 'no') '#' 904 write(iunit,'(a20)', advance = 'no') str_center('t', 19) 905 write(iunit,'(a20)', advance = 'yes') str_center('<O>(t)', 20) 906 write(iunit,'(a1)', advance = 'no') '#' 907 write(header_string, '(a)') '['//trim(units_abbrev(units_out%time))//']' 908 write(iunit,'(a20)', advance = 'no') str_center(trim(header_string), 19) 909 write(header_string, '(a)') '['//trim(units_abbrev(units_out%length))//']' 910 write(iunit,'(a20)', advance = 'yes') str_center(trim(header_string), 20) 911 912 do i = 0, time_steps 913 write(iunit, '(2e20.8)') units_from_atomic(units_out%time, i*dt), units_from_atomic(ot_unit, ot(i)) 914 end do 915 916 call io_close(iunit) 917 POP_SUB(write_ot) 918 end subroutine write_ot 919 920 921 ! --------------------------------------------------------- 922 subroutine read_ot(namespace, nspin, order, nw_subtracted) 923 type(namespace_t), intent(in) :: namespace 924 integer, intent(out) :: nspin 925 integer, intent(out) :: order 926 integer, intent(out) :: nw_subtracted 927 928 integer :: iunit, i, ierr 929 character(len=100) :: line 930 character(len=12) :: dummychar 931 FLOAT :: dummy, t1, t2 932 type(unit_t) :: ot_unit 933 934 PUSH_SUB(read_ot) 935 936 iunit = io_open('ot', namespace, action='read', status='old') 937 if(iunit == 0) then 938 write(message(1),'(a)') 'A file called ot should be present and was not found.' 939 call messages_fatal(1) 940 end if 941 942 read(iunit, '(15x,i2)') nspin 943 read(iunit, '(15x,i2)') order 944 read(iunit, '(28x,i2)') nw_subtracted 945 read(iunit, '(a)') line 946 947 i = index(line, 'Observable') 948 if(index(line, 'Observable') /= 0) then 949 observable(1) = -1 950 elseif(index(line, '# O =') /= 0) then 951 observable(1) = 0 952 if(index(line,'x') /= 0) then 953 observable(2) = 1 954 elseif(index(line,'y') /= 0) then 955 observable(2) = 2 956 elseif(index(line,'z') /= 0) then 957 observable(2) = 3 958 end if 959 elseif(index(line, '# (l, m) = ') /= 0) then 960 read(line,'(a12,i1,a1,i2,a1)') dummychar(1:12), observable(1), dummychar(1:1), observable(2), dummychar(1:1) 961 else 962 write(message(1),'(a)') 'Problem reading "ot" file: could not figure out the shape' 963 write(message(2),'(a)') 'of the observation operator.' 964 call messages_fatal(2) 965 end if 966 967 call kick_read(kick, iunit, namespace) 968 read(iunit, '(a)') line 969 read(iunit, '(a)') line 970 call io_skip_header(iunit) 971 972 ! Figure out about the units of the file 973 call unit_system_from_file(units, "ot", namespace, ierr) 974 if(ierr /= 0) then 975 write(message(1), '(a)') 'Could not figure out the units in file "ot".' 976 call messages_fatal(1) 977 end if 978 979 select case(observable(1)) 980 case(-1) 981 if(kick%n_multipoles > 0) then 982 ot_unit = units_out%length**kick%l(1) 983 else 984 ot_unit = units_out%length 985 end if 986 case(0) 987 ot_unit = units_out%length 988 case default 989 ot_unit = units_out%length**observable(1) 990 end select 991 992 ! count number of time_steps 993 time_steps = 0 994 do 995 read(iunit, *, end=100) dummy 996 time_steps = time_steps + 1 997 if(time_steps == 1) t1 = dummy 998 if(time_steps == 2) t2 = dummy 999 end do 1000 100 continue 1001 dt = units_to_atomic(units%time, (t2 - t1)) ! units_out is OK 1002 1003 call io_skip_header(iunit) 1004 1005 SAFE_ALLOCATE(ot(0:time_steps)) 1006 1007 do i = 0, time_steps-1 1008 read(iunit, *) dummy, ot(i) 1009 ot(i) = units_to_atomic(ot_unit, ot(i)) 1010 end do 1011 1012 POP_SUB(read_ot) 1013 end subroutine read_ot 1014 1015 1016 ! --------------------------------------------------------- 1017 subroutine print_omega_file(namespace, omega, search_interval, final_time, nfrequencies) 1018 type(namespace_t), intent(in) :: namespace 1019 FLOAT, intent(inout) :: omega 1020 FLOAT, intent(inout) :: search_interval 1021 FLOAT, intent(inout) :: final_time 1022 integer, intent(inout) :: nfrequencies 1023 1024 integer :: iunit, i, ierr, nspin, order, nw_subtracted 1025 logical :: file_exists 1026 character(len=20) :: header_string 1027 FLOAT, allocatable :: warray(:), tarray(:) 1028 FLOAT :: leftbound, rightbound, dw, w, aw 1029 1030 PUSH_SUB(print_omega_file) 1031 1032 ! First, let us check that the file "ot" exists. 1033 inquire(file="ot", exist = file_exists) 1034 if(.not.file_exists) then 1035 write(message(1),'(a)') "Could not find 'ot' file." 1036 call messages_fatal(1) 1037 end if 1038 1039 ! Now, we should find out which units the file "ot" has. 1040 call unit_system_from_file(units, "ot", namespace, ierr) 1041 if(ierr /= 0) then 1042 write(message(1),'(a)') "Could not retrieve units in the 'ot' file." 1043 call messages_fatal(1) 1044 end if 1045 1046 if(omega > M_ZERO) then 1047 omega = units_to_atomic(units%energy, omega) 1048 else 1049 omega = M_HALF 1050 end if 1051 1052 if(search_interval > M_ZERO) then 1053 search_interval = units_to_atomic(units%energy, search_interval) 1054 else 1055 search_interval = M_HALF 1056 end if 1057 1058 leftbound = omega - search_interval 1059 rightbound = omega + search_interval 1060 1061 SAFE_ALLOCATE(warray(1:nfrequencies)) 1062 SAFE_ALLOCATE(tarray(1:nfrequencies)) 1063 1064 call read_ot(namespace, nspin, order, nw_subtracted) 1065 1066 if(mod(order, 2) == 1) then 1067 mode = SINE_TRANSFORM 1068 else 1069 mode = COSINE_TRANSFORM 1070 end if 1071 1072 if(final_time > M_ZERO) then 1073 total_time = units_to_atomic(units%time, final_time) 1074 if(total_time > dt*time_steps) then 1075 total_time = dt*time_steps 1076 write(message(1), '(a)') 'The requested total time to process is larger than the time available in the input file.' 1077 write(message(2), '(a,f8.4,a)') 'The time has been adjusted to ', & 1078 units_from_atomic(units_out%time, total_time), trim(units_abbrev(units_out%time)) 1079 call messages_warning(2) 1080 end if 1081 time_steps = int(total_time / dt) 1082 total_time = time_steps * dt 1083 else 1084 total_time = dt*time_steps 1085 end if 1086 1087 warray = M_ZERO; tarray = M_ZERO 1088 dw = (rightbound-leftbound) / (nfrequencies - 1) 1089 do i = 1, nfrequencies 1090 w = leftbound + (i-1)*dw 1091 warray(i) = w 1092 call ft(w, aw) 1093 tarray(i) = aw 1094 end do 1095 1096 iunit = io_open('omega', namespace, action='write', status='replace') 1097 write(iunit, '(a15,i2)') '# nspin ', nspin 1098 call kick_write(kick, iunit) 1099 write(iunit, '(a)') '#%' 1100 write(iunit, '(a1,a20)', advance = 'no') '#', str_center("omega", 20) 1101 write(header_string,'(a)') 'F(omega)' 1102 write(iunit, '(a20)', advance = 'yes') str_center(trim(header_string), 20) 1103 write(iunit, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20) 1104 ! Here we should print the units of the transform. 1105 write(iunit, '(a)', advance = 'yes') 1106 1107 do i = 1, nfrequencies 1108 write(iunit,'(2e20.8)') units_from_atomic(units_out%energy, warray(i)), & 1109 tarray(i) 1110 end do 1111 1112 SAFE_DEALLOCATE_A(warray) 1113 SAFE_DEALLOCATE_A(tarray) 1114 1115 POP_SUB(print_omega_file) 1116 end subroutine print_omega_file 1117 ! --------------------------------------------------------- 1118 1119end module oscillator_strength_oct_m 1120 1121! --------------------------------------------------------- 1122! --------------------------------------------------------- 1123! --------------------------------------------------------- 1124program oscillator_strength 1125 use command_line_oct_m 1126 use global_oct_m 1127 use io_oct_m 1128 use messages_oct_m 1129 use oscillator_strength_oct_m 1130 1131 implicit none 1132 1133 integer :: run_mode, order, nfrequencies, ierr, nresonances 1134 FLOAT :: omega, search_interval, final_time, damping 1135 integer, parameter :: ANALYZE_NTHORDER_SIGNAL = 1, & 1136 GENERATE_NTHORDER_SIGNAL = 2, & 1137 READ_RESONANCES_FROM_FILE = 3, & 1138 GENERATE_OMEGA_FILE = 4 1139 character(len=100) :: ffile 1140 1141 ! Reads the information passed through the command line options (if available). 1142 call getopt_init(ierr) 1143 if(ierr /= 0) then 1144 message(1) = "Your Fortran compiler doesn't support command-line arguments;" 1145 message(2) = "the oct-oscillator-strength command is not available." 1146 call messages_fatal(2) 1147 end if 1148 1149 ! Set the default values 1150 run_mode = ANALYZE_NTHORDER_SIGNAL 1151 omega = - M_ONE 1152 search_interval = - M_ONE 1153 order = 1 1154 nfrequencies = 1000 1155 final_time = - M_ONE 1156 nresonances = 1 1157 observable(1) = -1 1158 observable(2) = 0 1159 ffile = "" 1160 damping = CNST(0.1)/CNST(27.2114) ! This is the usual damping factor used in the literature. 1161 1162 ! Get the parameters from the command line. 1163 call getopt_oscillator_strength(run_mode, omega, search_interval, & 1164 order, nresonances, nfrequencies, final_time, & 1165 observable(1), observable(2), damping, ffile) 1166 call getopt_end() 1167 1168 ! Initialize stuff 1169 call global_init(is_serial = .true.) 1170 call parser_init() 1171 call io_init(defaults = .true.) 1172 1173 select case(run_mode) 1174 case(GENERATE_NTHORDER_SIGNAL) 1175 call generate_signal(global_namespace, order, observable) 1176 case(ANALYZE_NTHORDER_SIGNAL) 1177 call analyze_signal(order, omega, search_interval, final_time, nresonances, nfrequencies, damping, & 1178 global_namespace) 1179 case(READ_RESONANCES_FROM_FILE) 1180 call read_resonances_file(order, ffile, global_namespace, search_interval, final_time, nfrequencies) 1181 case(GENERATE_OMEGA_FILE) 1182 call print_omega_file(global_namespace, omega, search_interval, final_time, nfrequencies) 1183 case default 1184 end select 1185 1186 call io_end() 1187 call parser_end() 1188 call global_end() 1189 1190end program oscillator_strength 1191! --------------------------------------------------------- 1192 1193!! Local Variables: 1194!! mode: f90 1195!! coding: utf-8 1196!! End: 1197