TABLE OF CONTENTS
- ABINIT/m_model_screening
- m_model_screening/find_peaks
- m_model_screening/im_screening
- m_model_screening/init_peaks_even_dist
- m_model_screening/init_peaks_from_grid
- m_model_screening/init_single_peak
- m_model_screening/print_peaks
- m_model_screening/re_and_im_screening
- m_model_screening/re_and_im_screening_with_phase
- m_model_screening/re_screening
- m_model_screening/remove_phase
- m_model_screening/sequential_fitting
ABINIT/m_model_screening [ Modules ]
NAME
m_model_screening
FUNCTION
Module containing functions for calculating and fitting model dielectric functions
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_model_screening 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 use m_io_tools, only : open_file 29 30 implicit none 31 32 private 33 34 public :: im_screening ! Calc. Drude-Lorentz model function from parameters. 35 public :: re_screening ! Calc. Drude-Lorentz model function from parameters. 36 public :: re_and_im_screening ! Calc. Drude-Lorentz model function from parameters. 37 public :: re_and_im_screening_with_phase ! Calc. Drude-Lorentz model function from parameters. 38 ! with the addition of a phase 39 public :: sequential_fitting ! Fit poles one by one 40 public :: init_peaks_from_grid ! find approximate expression for parameters from 41 ! chi0 or eps^-1 on a grid in the complex plane. 42 public :: init_peaks_even_dist ! Initial guess from even distributuin of peaks 43 public :: init_single_peak ! Initialise a single peak from the maximum, the 44 ! origin, and the second value along the 45 ! imaginary axis 46 ! public :: int_screening ! Find the integral along real or complex axis 47 ! from parameters. 48 public :: remove_phase 49 50 CONTAINS !==============================================================================
m_model_screening/find_peaks [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
find_peaks
FUNCTION
Find the location of the highest peaks along gridlines starting at the real axis and then moving towards higher imaginary frequencies. Stop when enough peaks to satisfy the number of poles needed has been found
INPUTS
omega = (complex) Real and imaginary part of the frequency points fvals = The function to be fitted in the complex plane nomega = number of fit points nfreqre = number or imaginary gridlines nfreqim = number of real gridlines
OUTPUT
NOTES
SOURCE
846 subroutine find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline) 847 848 implicit none 849 850 !Arguments ------------------------------------ 851 !scalars 852 integer,intent(in) :: nomega,nfreqre,nfreqim,npoles 853 integer, intent(inout) :: iline 854 !arrays 855 integer ,intent(inout) :: ploc(npoles) 856 complex(gwpc), intent(in) :: fval(nomega) 857 858 !Local variables------------------------------- 859 !scalars 860 integer :: ire,iim,ipoles 861 integer :: idx1,idx2,idx3,ipol 862 real(gwp) :: val1,val2,val3 863 864 !arrays 865 integer :: ploc_prev(npoles) 866 real :: pval(npoles),pval_prev(npoles) 867 868 ! ********************************************************************* 869 870 ploc=-1; ploc_prev=-1 871 pval=zero; pval_prev=zero; ipol=1; ipoles=0 872 873 ! First do a line along the real axis 874 idx1 = 1; idx2 = 2 875 val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2)) 876 if (ABS(val1)>ABS(val2)) then 877 ipoles = ipoles + 1 878 ploc(1)=idx1; pval(1)=val1 879 write(std_out,*) ' pval:',pval 880 write(std_out,*) ' ploc:',ploc 881 end if 882 do ire=2,nfreqre-3 883 idx1 = ire; idx2 = idx1+1; idx3 = idx1+2 884 val1 = AIMAG(fval(idx1)) 885 val2 = AIMAG(fval(idx2)) 886 val3 = AIMAG(fval(idx3)) 887 if (((val1<val2).AND.(val2>val3))) then 888 if (sign(1.0_gwp,val2)<zero) CYCLE 889 ipoles = ipoles + 1 890 if (ANY( ABS(pval(:))<ABS(val2) )) then 891 ipol = MINLOC(ABS(pval(:)),1) 892 ploc(ipol)=idx2; pval(ipol)=val2 893 write(std_out,*) ' pval:',pval 894 write(std_out,*) ' ploc:',ploc 895 end if 896 else if (((val1>val2).AND.(val2<val3))) then 897 if (sign(1.0_gwp,val2)>zero) CYCLE 898 ipoles = ipoles + 1 899 if (ANY( ABS(pval(:))<ABS(val2) )) then 900 ipol = MINLOC(ABS(pval(:)),1) 901 ploc(ipol)=idx2; pval(ipol)=val2 902 write(std_out,*) ' pval:',pval 903 write(std_out,*) ' ploc:',ploc 904 end if 905 end if 906 end do 907 ! Check last point 908 idx1 = nfreqre-2 909 idx2 = idx1 + 1 910 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 911 val1 = AIMAG(fval(idx1)) 912 val2 = AIMAG(fval(idx2)) 913 if (ABS(val1)<ABS(val2)) then 914 ipoles = ipoles + 1 915 if (ANY( ABS(pval(:))<ABS(val2) )) then 916 ipol = MINLOC(ABS(pval(:)),1) 917 ploc(ipol)=idx2; pval(ipol)=val2 918 write(std_out,*) ' pval:',pval 919 write(std_out,*) ' ploc:',ploc 920 end if 921 end if 922 write(std_out,'(a,i0)') ' Number of poles real axis:',ipoles 923 924 925 ploc_prev = ploc; pval_prev = pval 926 927 ! Do the rest of the imaginary grid until total 928 ! number of peaks found equals npoles or less 929 do iim=1,nfreqim-1 930 ploc=-1; pval=zero; ipol=1; ipoles=0 931 idx1 = nfreqre+iim 932 idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1) 933 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 934 val1 = AIMAG(fval(idx1)) 935 val2 = AIMAG(fval(idx2)) 936 if (ABS(val1)>ABS(val2)) then 937 ipoles = ipoles + 1 938 ploc(1)=idx1; pval(1)=val1 939 write(std_out,*) ' pval:',pval 940 write(std_out,*) ' ploc:',ploc 941 end if 942 ! Do all but the last 943 do ire=1,nfreqre-4 944 idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire 945 idx2 = idx1+1 946 idx3 = idx1+2 947 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3 948 val1 = AIMAG(fval(idx1)) 949 val2 = AIMAG(fval(idx2)) 950 val3 = AIMAG(fval(idx3)) 951 if (((val1<val2).AND.(val2>val3))) then 952 if (sign(1.0_gwp,val2)<zero) CYCLE 953 ipoles = ipoles + 1 954 if (ANY( ABS(pval(:))<ABS(val2) )) then 955 ipol = MINLOC(ABS(pval(:)),1) 956 ploc(ipol)=idx2; pval(ipol)=val2 957 write(std_out,*) ' pval:',pval 958 write(std_out,*) ' ploc:',ploc 959 end if 960 else if (((val1>val2).AND.(val2<val3))) then 961 if (sign(1.0_gwp,val2)>zero) CYCLE 962 ipoles = ipoles + 1 963 if (ANY( ABS(pval(:))<ABS(val2) )) then 964 ipol = MINLOC(ABS(pval(:)),1) 965 ploc(ipol)=idx2; pval(ipol)=val2 966 write(std_out,*) ' pval:',pval 967 write(std_out,*) ' ploc:',ploc 968 end if 969 end if 970 end do 971 ! Check last point 972 idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3 973 idx2 = idx1 + 1 974 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 975 val1 = AIMAG(fval(idx1)) 976 val2 = AIMAG(fval(idx2)) 977 if (ABS(val1)<ABS(val2)) then 978 ipoles = ipoles + 1 979 if (ANY( ABS(pval(:))<ABS(val2) )) then 980 ipol = MINLOC(ABS(pval(:)),1) 981 ploc(ipol)=idx2; pval(ipol)=val2 982 write(std_out,*) ' pval:',pval 983 write(std_out,*) ' ploc:',ploc 984 end if 985 end if 986 write(std_out,'(2(a,i0))') ' Line,:',iim,' ipoles:',ipoles 987 if (ipoles<=npoles) then 988 iline = iim - 1 989 ploc = ploc_prev 990 EXIT 991 end if 992 993 ploc_prev = ploc; pval_prev = pval 994 995 end do 996 997 end subroutine find_peaks
m_model_screening/im_screening [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
im_screening
FUNCTION
Return the imaginary part of model dielectric / inverse dielectric
function as given by a Drude-Lorentx model
The function is a sum in the complex plane:
f(z) = Sum_n f_n * Im[ ((w_n^2-z^2) - i*gamma*z)^-1 ], z=a-i*b
Here, f_n is the oscillator strength, w_n the location of the peak for
the imaginary function, and gamma is related to the width
INPUTS
omega = (complex) Real and imaginary part of the frequency points
coeff = The coefficients in order: f_1,w_1,gamma_1,f_2,w_2,gamma_2,
...,f_n,w_n,gamma_n
nomega = number of fit points
ncoeff = number of coefficients
OUTPUT
NOTES
SOURCE
80 subroutine im_screening(omega,fval,nomega,coeff,ncoeff) 81 82 implicit none 83 84 !Arguments ------------------------------------ 85 !scalars 86 integer,intent(in) :: nomega,ncoeff 87 !arrays 88 complex(dpc),intent(in) :: omega(nomega) 89 real(gwp) ,intent(in) :: coeff(ncoeff) 90 real(gwp) ,intent(out) :: fval(nomega) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: io,ip 95 real(gwp) :: rez,imz,realp,imagp 96 real(gwp) :: fn,wn,gamman 97 ! ********************************************************************* 98 99 ! The expression is: -f_n*(2*rez*imz-rez*gamma_n) 100 ! /( (-imz*gamma_n+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma_n)^2 ) 101 102 do io=1,nomega 103 fval(io) = 0.0 104 rez = REAL(omega(io)) 105 imz = AIMAG(omega(io)) 106 do ip=1,ncoeff,3 107 fn = coeff(ip) 108 wn = coeff(ip+1) 109 gamman = coeff(ip+2) 110 realp = -imz*gamman+wn*wn+imz*imz-rez*rez 111 imagp = rez*(two*imz-gamman) 112 113 fval(io) = fval(io)-fn*imagp/((realp*realp)+(imagp*imagp)) 114 115 end do 116 end do 117 118 end subroutine im_screening
m_model_screening/init_peaks_even_dist [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
init_peaks_even_dist
FUNCTION
Distribute the peaks evenly along a line in the complex plane and normalise.
INPUTS
omega = (complex) Real and imaginary part of the frequency points yvals = The function to be fitted in the complex plane coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients prtvol = Verbosity of diagnostics
OUTPUT
NOTES
SOURCE
634 subroutine init_peaks_even_dist(omega,fval,nomega,nfreqre,coeff,ncoeff,prtvol) 635 636 implicit none 637 638 !Arguments ------------------------------------ 639 !scalars 640 integer,intent(in) :: nomega,nfreqre,ncoeff,prtvol 641 !arrays 642 complex(dpc) ,intent(in) :: omega(nomega) 643 real(gwp) ,intent(out) :: coeff(ncoeff) 644 complex(gwpc),intent(in) :: fval(nomega) 645 646 !Local variables------------------------------- 647 !scalars 648 integer :: npoles,ip,idx,iw 649 real(gwp) :: delta,norm,div,val1,val2,osc,pol,gam 650 651 ! ********************************************************************* 652 653 npoles = ncoeff/3 654 div = real(npoles,gwp) 655 656 delta = (omega(nfreqre)-omega(1))/(div+1.0_gwp) 657 ! Integrate function along real axis (trapezoid rule) and have normalised 658 ! oscillator strengths 659 norm = fval(1)*half 660 do iw=2,nfreqre-1 661 norm = norm + fval(iw) 662 end do 663 norm = norm + fval(nfreqre)*half 664 norm = norm*(omega(nfreqre)-omega(1))/real(nfreqre,gwp) 665 norm = norm/div 666 667 do ip=1,npoles 668 idx = 3*(ip-1)+1 669 pol = delta*ip ! Position of maximum 670 gam = 0.1_gwp ! Spread of function 671 val2 = gam*half 672 val1 = SQRT(pol*pol-val2*val2) 673 osc = norm*val1*two 674 coeff(idx ) = osc!*(-1.0_gwp)**(ip-1) 675 coeff(idx+1) = pol ! Position of maximum 676 coeff(idx+2) = -gam ! Spread of function 677 if (prtvol>9) then 678 write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip 679 write(std_out,'(a,ES16.8)') ' Osc. strength:',osc 680 write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV' 681 write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV' 682 val2 = gam*half 683 val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2)) 684 write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV' 685 write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV' 686 write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1) 687 end if 688 end do 689 690 end subroutine init_peaks_even_dist
m_model_screening/init_peaks_from_grid [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
init_peaks_from_grid
FUNCTION
Find an initial guess of coefficents from the "valleys" and "hills" in the complex plane.
INPUTS
omega = (complex) Real and imaginary part of the frequency points yvals = The function to be fitted in the complex plane coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients prtvol = Verbosity of diagnostics
OUTPUT
NOTES
SOURCE
446 subroutine init_peaks_from_grid(omega,fval,nomega,nfreqre,nfreqim,coeff,ncoeff,prtvol) 447 448 implicit none 449 450 !Arguments ------------------------------------ 451 !scalars 452 integer,intent(in) :: nomega,nfreqre,nfreqim,ncoeff,prtvol 453 !arrays 454 complex(dpc) ,intent(in) :: omega(nomega) 455 real(gwp) ,intent(out) :: coeff(ncoeff) 456 complex(gwpc),intent(in) :: fval(nomega) 457 458 !Local variables------------------------------- 459 !scalars 460 integer :: npoles,iline,idx,ip 461 real(gwp) :: pol,gam,maxv,df,dk,val2,b2,osc,val1 462 real(gwp) :: temp1,temp2,temp3 463 464 !arrays 465 integer :: ploc(ncoeff/3) 466 467 ! ********************************************************************* 468 469 npoles = ncoeff/3 470 471 ! Map the evolution of peaks if prtvol>10 472 if (prtvol>10) then 473 call print_peaks(omega,fval,nomega,nfreqre,nfreqim) 474 end if ! prtvol>10 475 476 ! Count the number of peaks per line and find the location of the 477 ! constant-imaginary frequency line wich has at least a number of 478 ! peaks commensurate with the requested number of poles 479 call find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline) 480 write(std_out,*) ' Optimum peak locations:',ploc 481 write(std_out,*) ' on iline:',iline 482 ! Now fit the peaks. A linear interpolation along the imaginary 483 ! direction is used to get a rough estimate of the width of 484 ! the peak. 485 do ip=1,npoles 486 pol = REAL(omega(ploc(ip))) 487 maxv = AIMAG(fval(ploc(ip))) 488 write(std_out,*) ' maxv:',maxv 489 if (ploc(ip)<nfreqre+1) then ! We are right on the real axis 490 if (ploc(ip)==1) then ! Peak is at origin (Drude peak, i.e. metal) 491 b2 = AIMAG(omega(nfreqre+1)) 492 val2 = AIMAG(fval(nfreqre+1)) 493 write(std_out,*) '1: ploc:',ploc(ip),' b2:',b2,' val2:',val2 494 else ! Second value will be in c-plane 495 idx = nfreqre+nfreqim+ploc(ip)-1 496 b2 = AIMAG(omega(idx)) 497 val2 = AIMAG(fval(idx)) 498 write(std_out,*) '2: ploc:',ploc(ip),' b2:',b2,' val2:',val2 499 end if 500 else if (ploc(ip)<nfreqre+nfreqim+1) then ! We are right on the imaginary axis 501 if (ploc(ip)==nfreqre+nfreqim) then 502 ABI_ERROR(' Peak in upper left corner. This should never happen') 503 end if 504 b2 = AIMAG(omega(ploc(ip)+1)) 505 val2 = AIMAG(fval(ploc(ip)+1)) 506 write(std_out,*) '3: ploc:',ploc(ip),' b2:',b2,' val2:',val2 507 else ! We are in the complex plane 508 idx = ploc(ip)+nfreqre-1 509 b2 = AIMAG(omega(idx)) 510 val2 = AIMAG(fval(idx)) 511 write(std_out,*) '4: ploc:',ploc(ip),' idx:',idx,' b2:',b2,' val2:',val2 512 end if 513 df = ABS(val2 - maxv) 514 dk = df/b2 515 gam = -ABS(val2/dk) 516 !temp1 = SQRT(-b2*b2*val2*(val2+maxv)+(maxv*maxv)**2) 517 !temp2 = b2*(two*pol*pol+b2*b2)*val2-b2*maxv*pol*pol 518 !temp3 = (pol*pol+b2*b2)*val2+maxv*pol*pol 519 !gam = -((temp1-temp2)/temp3) 520 if (gam>zero) gam = ((temp1+temp2)/temp3) 521 osc = maxv*gam*pol 522 idx = 3*(ip-1)+1 523 coeff(idx ) = osc ! Oscillator strength 524 coeff(idx+1) = pol ! Position of maximum 525 coeff(idx+2) = gam ! Spread of function 526 if (prtvol>9) then 527 write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip 528 write(std_out,'(a,ES16.8)') ' Osc. strength:',osc 529 write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV' 530 write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV' 531 val2 = gam*half 532 val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2)) 533 write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV' 534 write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV' 535 write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1) 536 end if 537 end do 538 539 end subroutine init_peaks_from_grid
m_model_screening/init_single_peak [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
init_single_peak
FUNCTION
Initialise a single peak by using the behaviour along the imaginary axis and the main peak
INPUTS
omega = (complex) Real and imaginary part of the frequency points yvals = The function to be fitted in the complex plane coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients prtvol = Verbosity of diagnostics
OUTPUT
NOTES
SOURCE
564 subroutine init_single_peak(omega,refval,imfval,nomega,nfreqre,coeff,prtvol) 565 566 implicit none 567 568 !Arguments ------------------------------------ 569 !scalars 570 integer,intent(in) :: nomega,nfreqre,prtvol 571 !arrays 572 complex(dpc),intent(in) :: omega(nomega) 573 real(gwp) ,intent(out) :: coeff(3) 574 real(gwp) ,intent(in) :: refval(nomega),imfval(nomega) 575 576 !Local variables------------------------------- 577 !scalars 578 integer :: maxpos,idx 579 real(gwp) :: pol,osc,gam,val1,val2,pol_sq 580 581 ! ********************************************************************* 582 583 maxpos = MAXLOC(ABS(imfval(1:nfreqre)),1) 584 if (maxpos==1) maxpos = MAXLOC(ABS(imfval(2:nfreqre)),1) 585 pol = REAL(omega(maxpos)) 586 pol_sq = pol*pol 587 osc = -refval(1)*pol_sq 588 idx = nfreqre+1 589 val2 = refval(idx) 590 val1 = osc+val2*pol_sq+val2*AIMAG(omega(idx))*AIMAG(omega(idx)) 591 gam = -ABS(val1/(AIMAG(omega(idx))*val2)) 592 593 coeff(1) = osc 594 coeff(2) = pol 595 coeff(3) = gam 596 597 if (prtvol>9) then 598 write(std_out,'(a,ES16.8)') ' Osc. strength:',osc 599 write(std_out,'(a,ES16.8,a)') ' Peak location:',pol*Ha_eV,' eV' 600 write(std_out,'(a,ES16.8,a)') ' Peak width:',gam*Ha_eV,' eV' 601 val2 = gam*half 602 val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2)) 603 write(std_out,'(a,ES16.8,a)') ' Re[z] for pole:',val1*Ha_eV,' eV' 604 write(std_out,'(a,ES16.8,a)') ' Im[z] for pole:',val2*Ha_eV,' eV' 605 write(std_out,'(a,ES16.8)') ' Amplitude:',osc*half/ABS(val1) 606 end if 607 608 end subroutine init_single_peak
m_model_screening/print_peaks [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
print_peaks
FUNCTION
Find and output the location of peaks on the grid in a file
INPUTS
omega = (complex) Real and imaginary part of the frequency points fvals = The function to be fitted in the complex plane nomega = number of fit points nfreqre = number or imaginary gridlines nfreqim = number of real gridlines
OUTPUT
NOTES
SOURCE
714 subroutine print_peaks(omega,fval,nomega,nfreqre,nfreqim) 715 716 implicit none 717 718 !Arguments ------------------------------------ 719 !scalars 720 integer,intent(in) :: nomega,nfreqre,nfreqim 721 !arrays 722 complex(dpc) ,intent(in) :: omega(nomega) 723 complex(gwpc),intent(in) :: fval(nomega) 724 725 !Local variables------------------------------- 726 !scalars 727 integer :: ire,iim,unt_tmp 728 integer :: idx1,idx2,idx3 729 real(gwp) :: rez,imz,val1,val2,val3 730 character(len=500) :: msg 731 732 ! ********************************************************************* 733 734 if (open_file("grid_peak_tree.dat", msg, newunit=unt_tmp) /= 0) then 735 ABI_ERROR(msg) 736 end if 737 738 do iim=nfreqim,1,-1 739 ! write(std_out,*) ' iim:',iim 740 ! Check first points 741 idx1 = nfreqre+iim 742 idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1) 743 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 744 val1 = AIMAG(fval(idx1)) 745 val2 = AIMAG(fval(idx2)) 746 if (ABS(val1)>ABS(val2)) then 747 rez = REAL(omega(idx1)) 748 imz = AIMAG(omega(idx1)) 749 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1 750 end if 751 ! Do all but the last 752 do ire=1,nfreqre-4 753 idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire 754 idx2 = idx1+1 755 idx3 = idx1+2 756 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3 757 rez = REAL(omega(idx2)) 758 imz = AIMAG(omega(idx2)) 759 val1 = AIMAG(fval(idx1)) 760 val2 = AIMAG(fval(idx2)) 761 val3 = AIMAG(fval(idx3)) 762 if (((val1<val2).AND.(val2>val3))) then 763 if (sign(1.0_gwp,val2)<zero) CYCLE 764 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 765 else if (((val1>val2).AND.(val2<val3))) then 766 if (sign(1.0_gwp,val2)>zero) CYCLE 767 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 768 end if 769 end do 770 ! Check last point 771 idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3 772 idx2 = idx1 + 1 773 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 774 rez = REAL(omega(idx2)) 775 imz = AIMAG(omega(idx2)) 776 val1 = AIMAG(fval(idx1)) 777 val2 = AIMAG(fval(idx2)) 778 if (ABS(val1)<ABS(val2)) then 779 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 780 end if 781 end do 782 ! finally, do the purely real axis 783 ! Check first points 784 idx1 = 1; idx2 = 2 785 val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2)) 786 if (ABS(val1)>ABS(val2)) then 787 rez = REAL(omega(idx1)) 788 imz = AIMAG(omega(idx1)) 789 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1 790 end if 791 do ire=2,nfreqre-3 792 idx1 = ire; idx2 = idx1+1; idx3 = idx1+2 793 rez = REAL(omega(idx2)) 794 imz = AIMAG(omega(idx2)) 795 val1 = AIMAG(fval(idx1)) 796 val2 = AIMAG(fval(idx2)) 797 val3 = AIMAG(fval(idx3)) 798 if (((val1<val2).AND.(val2>val3))) then 799 if (sign(1.0_gwp,val2)<zero) CYCLE 800 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 801 else if (((val1>val2).AND.(val2<val3))) then 802 if (sign(1.0_gwp,val2)>zero) CYCLE 803 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 804 end if 805 end do 806 ! Check last point 807 idx1 = nfreqre-2 808 idx2 = idx1 + 1 809 ! write(std_out,*) ' idx1:',idx1,' idx2:',idx2 810 rez = REAL(omega(idx2)) 811 imz = AIMAG(omega(idx2)) 812 val1 = AIMAG(fval(idx1)) 813 val2 = AIMAG(fval(idx2)) 814 if (ABS(val1)<ABS(val2)) then 815 write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2 816 end if 817 818 close(unt_tmp) 819 820 end subroutine print_peaks
m_model_screening/re_and_im_screening [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
re_and_im_screening
FUNCTION
Return the real and imaginary part of model dielectric / inverse dielectric
function as evaluated from pole coefficients.
The function is a sum of poles in the complex plane:
f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
where each pole occurs twice in a time-ordered fashion.
Here, the A are the oscillator strengths, B the real component of the position
of the pole, and C the imaginary component.
INPUTS
omega = (complex) Real and imaginary part of the frequency points coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients
OUTPUT
NOTES
SOURCE
216 subroutine re_and_im_screening(omega,fval,nomega,coeff,ncoeff) 217 218 implicit none 219 220 !Arguments ------------------------------------ 221 !scalars 222 integer,intent(in) :: nomega,ncoeff 223 !arrays 224 complex(dpc) ,intent(in) :: omega(nomega) 225 real(gwp) ,intent(in) :: coeff(ncoeff) 226 complex(gwp),intent(out) :: fval(nomega) 227 228 !Local variables------------------------------- 229 !scalars 230 integer :: io,ip 231 real(gwp) :: rez,imz,realp,imagp,refval,imfval 232 real(gwp) :: fn,wn,gamman 233 ! ********************************************************************* 234 235 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2) 236 ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 ) 237 238 do io=1,nomega 239 fval(io) = 0.0 240 rez = REAL(omega(io)) 241 imz = AIMAG(omega(io)) 242 do ip=1,ncoeff,3 243 fn = coeff(ip) 244 wn = coeff(ip+1) 245 gamman = coeff(ip+2) 246 realp = -imz*gamman+wn*wn+imz*imz-rez*rez 247 imagp = rez*(two*imz-gamman) 248 249 refval = fn*realp/((realp*realp)+(imagp*imagp)) 250 imfval = fn*imagp/((realp*realp)+(imagp*imagp)) 251 252 fval(io) = fval(io)-CMPLX(refval,imfval) 253 254 end do 255 end do 256 257 end subroutine re_and_im_screening
m_model_screening/re_and_im_screening_with_phase [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
re_and_im_screening_with_phase
FUNCTION
Return the real and imaginary part of model dielectric / inverse dielectric
function as evaluated from pole coefficients.
The function is a sum of poles in the complex plane:
f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
where each pole occurs twice in a time-ordered fashion.
Here, the A are the oscillator strengths, B the real component of the position
of the pole, and C the imaginary component.
INPUTS
omega = (complex) Real and imaginary part of the frequency points coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients
OUTPUT
NOTES
SOURCE
287 subroutine re_and_im_screening_with_phase(omega,fval,nomega,coeff,ncoeff) 288 289 implicit none 290 291 !Arguments ------------------------------------ 292 !scalars 293 integer,intent(in) :: nomega,ncoeff 294 !arrays 295 complex(dpc) ,intent(in) :: omega(nomega) 296 real(gwp) ,intent(in) :: coeff(ncoeff) 297 complex(gwpc),intent(out) :: fval(nomega) 298 299 !Local variables------------------------------- 300 !scalars 301 integer :: io,ip,npoles 302 real(gwp) :: rez,imz,realp,imagp,refval,imfval,retemp,imtemp 303 real(gwp) :: fn,wn,gamman,imrot,rerot 304 ! ********************************************************************* 305 306 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2) 307 ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 ) 308 npoles = (ncoeff-1)/3 309 310 do io=1,nomega 311 fval(io) = 0.0 312 rez = REAL(omega(io)) 313 imz = AIMAG(omega(io)) 314 do ip=1,(ncoeff-1),3 315 fn = coeff(ip) 316 wn = coeff(ip+1) 317 gamman = coeff(ip+2) 318 realp = -imz*gamman+wn*wn+imz*imz-rez*rez 319 imagp = rez*(two*imz-gamman) 320 321 refval = fn*realp/((realp*realp)+(imagp*imagp)) 322 imfval = fn*imagp/((realp*realp)+(imagp*imagp)) 323 324 fval(io) = fval(io)-CMPLX(refval,imfval) 325 326 end do 327 ! Restore phase 328 rerot = COS(coeff(npoles*3+1)) 329 imrot = SIN(coeff(npoles*3+1)) 330 retemp = REAL(fval(io)) 331 imtemp = AIMAG(fval(io)) 332 fval(io) = CMPLX(rerot*retemp-imrot*imtemp,rerot*imtemp + imrot*retemp) 333 end do 334 335 end subroutine re_and_im_screening_with_phase
m_model_screening/re_screening [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
re_screening
FUNCTION
Return the real part of model dielectric / inverse dielectric
function as evaluated from pole coefficients.
The function is a sum of poles in the complex plane:
f(z) = Sum_n[ A/(z-(B-iC)) - A/(z-(-B+iC)) ],
where each pole occurs twice in a time-ordered fashion.
Here, the A are the oscillator strengths, B the real component of the position
of the pole, and C the imaginary component.
INPUTS
omega = (complex) Real and imaginary part of the frequency points coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n nomega = number of fit points ncoeff = number of coefficients
OUTPUT
NOTES
SOURCE
148 subroutine re_screening(omega,fval,nomega,coeff,ncoeff) 149 150 implicit none 151 152 !Arguments ------------------------------------ 153 !scalars 154 integer,intent(in) :: nomega,ncoeff 155 !arrays 156 complex(dpc),intent(in) :: omega(nomega) 157 real(gwp) ,intent(in) :: coeff(ncoeff) 158 real(gwp) ,intent(out) :: fval(nomega) 159 160 !Local variables------------------------------- 161 !scalars 162 integer :: io,ip 163 real(gwp) :: rez,imz,realp,imagp 164 real(gwp) :: fn,wn,gamman 165 ! ********************************************************************* 166 167 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2) 168 ! /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 ) 169 170 do io=1,nomega 171 fval(io) = 0.0 172 rez = REAL(omega(io)) 173 imz = AIMAG(omega(io)) 174 do ip=1,ncoeff,3 175 fn = coeff(ip) 176 wn = coeff(ip+1) 177 gamman = coeff(ip+2) 178 realp = -imz*gamman+wn*wn+imz*imz-rez*rez 179 imagp = rez*(two*imz-gamman) 180 181 fval(io) = fval(io)-fn*realp/((realp*realp)+(imagp*imagp)) 182 183 end do 184 end do 185 186 end subroutine re_screening
m_model_screening/remove_phase [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
remove_phase
FUNCTION
Find out what the complex phase factor is for off-diagonal elements and unmix the components.
INPUTS
fvals = The function to be fitted in the complex plane. nomega = Number of fit points. nfreqre = number or imaginary gridlines phase = The phase angle.
OUTPUT
NOTES
SOURCE
1020 subroutine remove_phase(fval,nomega,phase) 1021 1022 implicit none 1023 1024 !Arguments ------------------------------------ 1025 !scalars 1026 integer, intent(in) :: nomega 1027 real(gwp), intent(out) :: phase 1028 !arrays 1029 complex(gwpc), intent(inout) :: fval(nomega) 1030 1031 !Local variables------------------------------- 1032 !scalars 1033 integer :: io 1034 real(gwp) :: a,b,retemp,imtemp 1035 1036 ! ********************************************************************* 1037 1038 ! The phase can be found by checking when the function is 1039 ! identically zero along the imaginary axis 1040 if (ABS(AIMAG(fval(1)))<tol14) then ! Phase is zero 1041 phase = zero 1042 RETURN 1043 else if (ABS(REAL(fval(1)))<tol14) then ! Phase is exactly pi/2 1044 phase = pi*half 1045 a = zero 1046 b = -1.0_gwp 1047 else 1048 phase = ATAN(AIMAG(fval(1))/REAL(fval(1))) 1049 a = COS(phase) 1050 b = -SIN(phase) 1051 end if 1052 ! Rotate values 1053 do io=1,nomega 1054 retemp = REAL(fval(io)) 1055 imtemp = AIMAG(fval(io)) 1056 fval(io) = CMPLX(a*retemp-b*imtemp,a*imtemp+b*retemp) 1057 end do 1058 1059 end subroutine remove_phase
m_model_screening/sequential_fitting [ Functions ]
[ Top ] [ m_model_screening ] [ Functions ]
NAME
sequential_fitting
FUNCTION
Fit a function in the complex plane pole-by-pole in such a way as to increasingly minimise the error
INPUTS
omega = (complex) Real and imaginary part of the frequency points refval = Real part of function to be fitted imfval = imaginary part of function to be fitted nomega = Total number of points in the complex plane nfreqre = Number of points along real axis nfreqim = Number of points along imaginary axis coeff = The coefficients in order: A_1,B_1,C_1,A_2,B_2,C_2,...,A_n,B_n,C_n ncoeff = number of coefficients prtvol = Diagnostics verbose level
OUTPUT
NOTES
SOURCE
364 subroutine sequential_fitting(omega,refval,imfval,nomega,nfreqre,coeff,& 365 & ncoeff,prtvol,startcoeff) 366 367 implicit none 368 369 !Arguments ------------------------------------ 370 !scalars 371 integer,intent(in) :: nomega,nfreqre,ncoeff,prtvol 372 !arrays 373 complex(dpc),intent(in) :: omega(nomega) 374 real(gwp) ,intent(out) :: coeff(ncoeff) 375 real(gwp) ,intent(inout) :: refval(nomega),imfval(nomega) 376 real(gwp),optional,intent(out) :: startcoeff(ncoeff) 377 378 !Local variables------------------------------- 379 !scalars 380 integer :: ip,npoles,idx 381 real(gwp) :: thiscoeff(3),norm,invnorm 382 real(dpc) :: re_zvals(nomega),im_zvals(nomega) 383 ! real(dpc) :: orig_refval(nomega),orig_imfval(nomega) 384 complex(gwp) :: pole_func(nomega) 385 ! ********************************************************************* 386 387 npoles = ncoeff/3 388 re_zvals(:) = REAL(omega(:)) 389 im_zvals(:) = AIMAG(omega(:)) 390 391 ! Normalise 392 norm = MAXVAL(ABS(imfval)) 393 invnorm = 1.0_gwp/norm 394 refval = invnorm*refval 395 imfval = invnorm*imfval 396 397 ! Loop over poles to fit 398 do ip=1,npoles 399 idx = 3*(ip-1)+1 400 ! Initialise pole 401 call init_single_peak(omega,refval,imfval,nomega,nfreqre,thiscoeff,prtvol) 402 if (present(startcoeff)) then 403 startcoeff(idx:idx+2) = thiscoeff(1:3) 404 end if 405 ! Make fit 406 #ifdef HAVE_LEVMAR 407 call dfit_re_and_im_screening(re_zvals,im_zvals,imfval,refval,& 408 & nomega,3,thiscoeff,prtvol) 409 #else 410 ABI_ERROR(' ABINIT was not compiled with the levmar library!') 411 #endif 412 ! Remove current fit 413 call re_and_im_screening(omega,pole_func,nomega,thiscoeff,3) 414 refval(:) = refval(:) - REAL(pole_func(:)) 415 imfval(:) = imfval(:) - AIMAG(pole_func(:)) 416 coeff(idx:idx+2) = thiscoeff(1:3) 417 coeff(idx) = norm*coeff(idx) 418 end do 419 420 end subroutine sequential_fitting