TABLE OF CONTENTS


ABINIT/m_model_screening [ Modules ]

[ Top ] [ Modules ]

NAME

 m_model_screening

FUNCTION

  Module containing functions for calculating and fitting model dielectric functions

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_model_screening
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_io_tools,       only : open_file
34 
35  implicit none
36 
37  private
38 
39  public :: im_screening         ! Calc. Drude-Lorentz model function from parameters.
40  public :: re_screening         ! Calc. Drude-Lorentz model function from parameters.
41  public :: re_and_im_screening  ! Calc. Drude-Lorentz model function from parameters.
42  public :: re_and_im_screening_with_phase  ! Calc. Drude-Lorentz model function from parameters.
43                                            !  with the addition of a phase
44  public :: sequential_fitting   ! Fit poles one by one
45  public :: init_peaks_from_grid ! find approximate expression for parameters from
46                                 !  chi0 or eps^-1 on a grid in the complex plane.
47  public :: init_peaks_even_dist ! Initial guess from even distributuin of peaks
48  public :: init_single_peak     ! Initialise a single peak from the maximum, the
49                                 !  origin, and the second value along the
50                                 !  imaginary axis
51 ! public :: int_screening        ! Find the integral along real or complex axis
52                                 !  from parameters.
53  public :: remove_phase
54 
55 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

PARENTS

      m_model_screening

CHILDREN

SOURCE

 959 subroutine find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline)
 960 
 961 
 962 !This section has been created automatically by the script Abilint (TD).
 963 !Do not modify the following lines by hand.
 964 #undef ABI_FUNC
 965 #define ABI_FUNC 'find_peaks'
 966 !End of the abilint section
 967 
 968   implicit none
 969 
 970 !Arguments ------------------------------------
 971 !scalars
 972   integer,intent(in)     :: nomega,nfreqre,nfreqim,npoles
 973   integer, intent(inout) :: iline
 974 !arrays
 975   integer    ,intent(inout) :: ploc(npoles)
 976   complex(gwpc), intent(in) :: fval(nomega)
 977 
 978 !Local variables-------------------------------
 979 !scalars
 980   integer    :: ire,iim,ipoles
 981   integer    :: idx1,idx2,idx3,ipol
 982   real(gwp) :: val1,val2,val3
 983 
 984 !arrays
 985   integer :: ploc_prev(npoles)
 986   real    :: pval(npoles),pval_prev(npoles)
 987 
 988 ! *********************************************************************
 989 
 990   ploc=-1; ploc_prev=-1
 991   pval=zero; pval_prev=zero; ipol=1; ipoles=0
 992 
 993 ! First do a line along the real axis
 994   idx1 = 1; idx2 = 2
 995   val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2))
 996   if (ABS(val1)>ABS(val2)) then
 997     ipoles = ipoles + 1
 998     ploc(1)=idx1; pval(1)=val1
 999     write(std_out,*) ' pval:',pval
1000     write(std_out,*) ' ploc:',ploc
1001   end if
1002   do ire=2,nfreqre-3
1003     idx1 = ire; idx2 = idx1+1; idx3 = idx1+2
1004     val1 = AIMAG(fval(idx1))
1005     val2 = AIMAG(fval(idx2))
1006     val3 = AIMAG(fval(idx3))
1007     if (((val1<val2).AND.(val2>val3))) then
1008       if (sign(1.0_gwp,val2)<zero) CYCLE
1009       ipoles = ipoles + 1
1010       if (ANY( ABS(pval(:))<ABS(val2) )) then
1011         ipol = MINLOC(ABS(pval(:)),1)
1012         ploc(ipol)=idx2; pval(ipol)=val2
1013         write(std_out,*) ' pval:',pval
1014         write(std_out,*) ' ploc:',ploc
1015       end if
1016     else if (((val1>val2).AND.(val2<val3))) then
1017       if (sign(1.0_gwp,val2)>zero) CYCLE
1018       ipoles = ipoles + 1
1019       if (ANY( ABS(pval(:))<ABS(val2) )) then
1020         ipol = MINLOC(ABS(pval(:)),1)
1021         ploc(ipol)=idx2; pval(ipol)=val2
1022         write(std_out,*) ' pval:',pval
1023         write(std_out,*) ' ploc:',ploc
1024       end if
1025     end if
1026   end do
1027 ! Check last point
1028   idx1 = nfreqre-2
1029   idx2 = idx1 + 1
1030 !  write(std_out,*) ' idx1:',idx1,' idx2:',idx2
1031   val1 = AIMAG(fval(idx1))
1032   val2 = AIMAG(fval(idx2))
1033   if (ABS(val1)<ABS(val2)) then
1034     ipoles = ipoles + 1
1035     if (ANY( ABS(pval(:))<ABS(val2) )) then
1036       ipol = MINLOC(ABS(pval(:)),1)
1037       ploc(ipol)=idx2; pval(ipol)=val2
1038        write(std_out,*) ' pval:',pval
1039        write(std_out,*) ' ploc:',ploc
1040     end if
1041   end if
1042   write(std_out,'(a,i0)') ' Number of poles real axis:',ipoles
1043 
1044 
1045   ploc_prev = ploc; pval_prev = pval
1046 
1047 ! Do the rest of the imaginary grid until total
1048 !  number of peaks found equals npoles or less
1049   do iim=1,nfreqim-1
1050     ploc=-1; pval=zero; ipol=1; ipoles=0
1051     idx1 = nfreqre+iim
1052     idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)
1053 !    write(std_out,*) ' idx1:',idx1,' idx2:',idx2
1054     val1 = AIMAG(fval(idx1))
1055     val2 = AIMAG(fval(idx2))
1056     if (ABS(val1)>ABS(val2)) then
1057       ipoles = ipoles + 1
1058       ploc(1)=idx1; pval(1)=val1
1059       write(std_out,*) ' pval:',pval
1060       write(std_out,*) ' ploc:',ploc
1061     end if
1062 !   Do all but the last
1063     do ire=1,nfreqre-4
1064        idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire
1065        idx2 = idx1+1
1066        idx3 = idx1+2
1067 !       write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3
1068        val1 = AIMAG(fval(idx1))
1069        val2 = AIMAG(fval(idx2))
1070        val3 = AIMAG(fval(idx3))
1071        if (((val1<val2).AND.(val2>val3))) then
1072          if (sign(1.0_gwp,val2)<zero) CYCLE
1073          ipoles = ipoles + 1
1074          if (ANY( ABS(pval(:))<ABS(val2) )) then
1075            ipol = MINLOC(ABS(pval(:)),1)
1076            ploc(ipol)=idx2; pval(ipol)=val2
1077            write(std_out,*) ' pval:',pval
1078            write(std_out,*) ' ploc:',ploc
1079          end if
1080        else if (((val1>val2).AND.(val2<val3))) then
1081          if (sign(1.0_gwp,val2)>zero) CYCLE
1082          ipoles = ipoles + 1
1083          if (ANY( ABS(pval(:))<ABS(val2) )) then
1084            ipol = MINLOC(ABS(pval(:)),1)
1085            ploc(ipol)=idx2; pval(ipol)=val2
1086            write(std_out,*) ' pval:',pval
1087            write(std_out,*) ' ploc:',ploc
1088          end if
1089        end if
1090     end do
1091 !   Check last point
1092     idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3
1093     idx2 = idx1 + 1
1094 !   write(std_out,*) ' idx1:',idx1,' idx2:',idx2
1095     val1 = AIMAG(fval(idx1))
1096     val2 = AIMAG(fval(idx2))
1097     if (ABS(val1)<ABS(val2)) then
1098       ipoles = ipoles + 1
1099       if (ANY( ABS(pval(:))<ABS(val2) )) then
1100         ipol = MINLOC(ABS(pval(:)),1)
1101         ploc(ipol)=idx2; pval(ipol)=val2
1102          write(std_out,*) ' pval:',pval
1103          write(std_out,*) ' ploc:',ploc
1104       end if
1105     end if
1106     write(std_out,'(2(a,i0))') ' Line,:',iim,' ipoles:',ipoles
1107     if (ipoles<=npoles) then
1108       iline = iim - 1
1109       ploc = ploc_prev
1110       EXIT
1111     end if
1112 
1113     ploc_prev = ploc; pval_prev = pval
1114 
1115  end do
1116 
1117 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

PARENTS

CHILDREN

SOURCE

 89 subroutine im_screening(omega,fval,nomega,coeff,ncoeff)
 90 
 91 
 92 !This section has been created automatically by the script Abilint (TD).
 93 !Do not modify the following lines by hand.
 94 #undef ABI_FUNC
 95 #define ABI_FUNC 'im_screening'
 96 !End of the abilint section
 97 
 98   implicit none
 99 
100 !Arguments ------------------------------------
101 !scalars
102   integer,intent(in)   :: nomega,ncoeff
103 !arrays
104   complex(dpc),intent(in)  :: omega(nomega)
105   real(gwp)  ,intent(in)  :: coeff(ncoeff)
106   real(gwp)  ,intent(out) :: fval(nomega)
107 
108 !Local variables-------------------------------
109 !scalars
110   integer :: io,ip
111   real(gwp) :: rez,imz,realp,imagp
112   real(gwp) :: fn,wn,gamman
113 ! *********************************************************************
114 
115 ! The expression is: -f_n*(2*rez*imz-rez*gamma_n)
116 !    /( (-imz*gamma_n+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma_n)^2 )
117 
118   do io=1,nomega
119     fval(io) = 0.0
120     rez =  REAL(omega(io))
121     imz = AIMAG(omega(io))
122     do ip=1,ncoeff,3
123       fn     = coeff(ip)
124       wn     = coeff(ip+1)
125       gamman = coeff(ip+2)
126       realp  = -imz*gamman+wn*wn+imz*imz-rez*rez
127       imagp  = rez*(two*imz-gamman)
128 
129         fval(io) = fval(io)-fn*imagp/((realp*realp)+(imagp*imagp))
130 
131     end do
132   end do
133 
134 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

PARENTS

CHILDREN

SOURCE

723 subroutine init_peaks_even_dist(omega,fval,nomega,nfreqre,coeff,ncoeff,prtvol)
724 
725 
726 !This section has been created automatically by the script Abilint (TD).
727 !Do not modify the following lines by hand.
728 #undef ABI_FUNC
729 #define ABI_FUNC 'init_peaks_even_dist'
730 !End of the abilint section
731 
732   implicit none
733 
734 !Arguments ------------------------------------
735 !scalars
736   integer,intent(in)   :: nomega,nfreqre,ncoeff,prtvol
737 !arrays
738   complex(dpc) ,intent(in)  :: omega(nomega)
739   real(gwp)   ,intent(out) :: coeff(ncoeff)
740   complex(gwpc),intent(in)  :: fval(nomega)
741 
742 !Local variables-------------------------------
743 !scalars
744   integer :: npoles,ip,idx,iw
745   real(gwp) :: delta,norm,div,val1,val2,osc,pol,gam
746 
747 ! *********************************************************************
748 
749   npoles = ncoeff/3
750   div = real(npoles,gwp)
751 
752   delta = (omega(nfreqre)-omega(1))/(div+1.0_gwp)
753   ! Integrate function along real axis (trapezoid rule) and have normalised
754   ! oscillator strengths
755   norm = fval(1)*half
756   do iw=2,nfreqre-1
757     norm = norm + fval(iw)
758   end do
759   norm = norm + fval(nfreqre)*half
760   norm = norm*(omega(nfreqre)-omega(1))/real(nfreqre,gwp)
761   norm = norm/div
762 
763   do ip=1,npoles
764     idx = 3*(ip-1)+1
765     pol = delta*ip   ! Position of maximum
766     gam = 0.1_gwp ! Spread of function
767     val2 = gam*half
768     val1 = SQRT(pol*pol-val2*val2)
769     osc = norm*val1*two
770     coeff(idx  ) = osc!*(-1.0_gwp)**(ip-1)
771     coeff(idx+1) = pol   ! Position of maximum
772     coeff(idx+2) = -gam ! Spread of function
773     if (prtvol>9) then
774       write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip
775       write(std_out,'(a,ES16.8)')    '  Osc. strength:',osc
776       write(std_out,'(a,ES16.8,a)')  '  Peak location:',pol*Ha_eV,' eV'
777       write(std_out,'(a,ES16.8,a)')  '     Peak width:',gam*Ha_eV,' eV'
778       val2 = gam*half
779       val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
780       write(std_out,'(a,ES16.8,a)')  ' Re[z] for pole:',val1*Ha_eV,' eV'
781       write(std_out,'(a,ES16.8,a)')  ' Im[z] for pole:',val2*Ha_eV,' eV'
782       write(std_out,'(a,ES16.8)')    '      Amplitude:',osc*half/ABS(val1)
783     end if
784   end do
785 
786 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

PARENTS

CHILDREN

SOURCE

512 subroutine init_peaks_from_grid(omega,fval,nomega,nfreqre,nfreqim,coeff,ncoeff,prtvol)
513 
514 
515 !This section has been created automatically by the script Abilint (TD).
516 !Do not modify the following lines by hand.
517 #undef ABI_FUNC
518 #define ABI_FUNC 'init_peaks_from_grid'
519 !End of the abilint section
520 
521   implicit none
522 
523 !Arguments ------------------------------------
524 !scalars
525   integer,intent(in)   :: nomega,nfreqre,nfreqim,ncoeff,prtvol
526 !arrays
527   complex(dpc) ,intent(in)  :: omega(nomega)
528   real(gwp)   ,intent(out) :: coeff(ncoeff)
529   complex(gwpc),intent(in)  :: fval(nomega)
530 
531 !Local variables-------------------------------
532 !scalars
533   integer :: npoles,iline,idx,ip
534   real(gwp) :: pol,gam,maxv,df,dk,val2,b2,osc,val1
535   real(gwp) :: temp1,temp2,temp3
536 
537 !arrays
538   integer :: ploc(ncoeff/3)
539 
540 ! *********************************************************************
541 
542   npoles = ncoeff/3
543 
544 ! Map the evolution of peaks if prtvol>10
545   if (prtvol>10) then
546     call print_peaks(omega,fval,nomega,nfreqre,nfreqim)
547   end if ! prtvol>10
548 
549 ! Count the number of peaks per line and find the location of the
550 ! constant-imaginary frequency line wich has at least a number of
551 ! peaks commensurate with the requested number of poles
552   call find_peaks(fval,nomega,nfreqre,nfreqim,ploc,npoles,iline)
553   write(std_out,*) ' Optimum peak locations:',ploc
554   write(std_out,*) '               on iline:',iline
555 ! Now fit the peaks. A linear interpolation along the imaginary
556 !  direction is used to get a rough estimate of the width of
557 !  the peak.
558   do ip=1,npoles
559     pol  = REAL(omega(ploc(ip)))
560     maxv = AIMAG(fval(ploc(ip)))
561     write(std_out,*) ' maxv:',maxv
562     if (ploc(ip)<nfreqre+1) then ! We are right on the real axis
563       if (ploc(ip)==1) then ! Peak is at origin (Drude peak, i.e. metal)
564         b2   = AIMAG(omega(nfreqre+1))
565         val2 = AIMAG(fval(nfreqre+1))
566         write(std_out,*) '1: ploc:',ploc(ip),' b2:',b2,' val2:',val2
567       else ! Second value will be in c-plane
568         idx  = nfreqre+nfreqim+ploc(ip)-1
569         b2   = AIMAG(omega(idx))
570         val2 = AIMAG(fval(idx))
571         write(std_out,*) '2: ploc:',ploc(ip),' b2:',b2,' val2:',val2
572       end if
573     else if (ploc(ip)<nfreqre+nfreqim+1) then ! We are right on the imaginary axis
574       if (ploc(ip)==nfreqre+nfreqim) then
575         MSG_ERROR(' Peak in upper left corner. This should never happen')
576       end if
577       b2   = AIMAG(omega(ploc(ip)+1))
578       val2 = AIMAG(fval(ploc(ip)+1))
579         write(std_out,*) '3: ploc:',ploc(ip),' b2:',b2,' val2:',val2
580     else ! We are in the complex plane
581       idx  = ploc(ip)+nfreqre-1
582       b2   = AIMAG(omega(idx))
583       val2 = AIMAG(fval(idx))
584       write(std_out,*) '4: ploc:',ploc(ip),' idx:',idx,' b2:',b2,' val2:',val2
585     end if
586     df = ABS(val2 - maxv)
587     dk = df/b2
588     gam = -ABS(val2/dk)
589     !temp1 = SQRT(-b2*b2*val2*(val2+maxv)+(maxv*maxv)**2)
590     !temp2 = b2*(two*pol*pol+b2*b2)*val2-b2*maxv*pol*pol
591     !temp3 = (pol*pol+b2*b2)*val2+maxv*pol*pol
592     !gam = -((temp1-temp2)/temp3)
593     if (gam>zero) gam = ((temp1+temp2)/temp3)
594     osc = maxv*gam*pol
595     idx = 3*(ip-1)+1
596     coeff(idx  ) = osc ! Oscillator strength
597     coeff(idx+1) = pol ! Position of maximum
598     coeff(idx+2) = gam ! Spread of function
599     if (prtvol>9) then
600       write(std_out,'(a,a,i0)') ch10,' Pole no,: ',ip
601       write(std_out,'(a,ES16.8)')    '  Osc. strength:',osc
602       write(std_out,'(a,ES16.8,a)')  '  Peak location:',pol*Ha_eV,' eV'
603       write(std_out,'(a,ES16.8,a)')  '     Peak width:',gam*Ha_eV,' eV'
604       val2 = gam*half
605       val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
606       write(std_out,'(a,ES16.8,a)')  ' Re[z] for pole:',val1*Ha_eV,' eV'
607       write(std_out,'(a,ES16.8,a)')  ' Im[z] for pole:',val2*Ha_eV,' eV'
608       write(std_out,'(a,ES16.8)')    '      Amplitude:',osc*half/ABS(val1)
609     end if
610   end do
611 
612 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

PARENTS

      m_model_screening

CHILDREN

SOURCE

642 subroutine init_single_peak(omega,refval,imfval,nomega,nfreqre,coeff,prtvol)
643 
644 
645 !This section has been created automatically by the script Abilint (TD).
646 !Do not modify the following lines by hand.
647 #undef ABI_FUNC
648 #define ABI_FUNC 'init_single_peak'
649 !End of the abilint section
650 
651   implicit none
652 
653 !Arguments ------------------------------------
654 !scalars
655   integer,intent(in)   :: nomega,nfreqre,prtvol
656 !arrays
657   complex(dpc),intent(in)  :: omega(nomega)
658   real(gwp)  ,intent(out) :: coeff(3)
659   real(gwp)  ,intent(in) :: refval(nomega),imfval(nomega)
660 
661 !Local variables-------------------------------
662 !scalars
663   integer :: maxpos,idx
664   real(gwp) :: pol,osc,gam,val1,val2,pol_sq
665 
666 ! *********************************************************************
667 
668   maxpos = MAXLOC(ABS(imfval(1:nfreqre)),1)
669   if (maxpos==1) maxpos = MAXLOC(ABS(imfval(2:nfreqre)),1)
670   pol = REAL(omega(maxpos))
671   pol_sq = pol*pol
672   osc = -refval(1)*pol_sq
673   idx = nfreqre+1
674   val2 = refval(idx)
675   val1 = osc+val2*pol_sq+val2*AIMAG(omega(idx))*AIMAG(omega(idx))
676   gam  = -ABS(val1/(AIMAG(omega(idx))*val2))
677 
678   coeff(1) = osc
679   coeff(2) = pol
680   coeff(3) = gam
681 
682   if (prtvol>9) then
683     write(std_out,'(a,ES16.8)')    '  Osc. strength:',osc
684     write(std_out,'(a,ES16.8,a)')  '  Peak location:',pol*Ha_eV,' eV'
685     write(std_out,'(a,ES16.8,a)')  '     Peak width:',gam*Ha_eV,' eV'
686     val2 = gam*half
687     val1 = SIGN(1.0_gwp,pol*pol-val2*val2)*SQRT(ABS(pol*pol-val2*val2))
688     write(std_out,'(a,ES16.8,a)')  ' Re[z] for pole:',val1*Ha_eV,' eV'
689     write(std_out,'(a,ES16.8,a)')  ' Im[z] for pole:',val2*Ha_eV,' eV'
690     write(std_out,'(a,ES16.8)')    '      Amplitude:',osc*half/ABS(val1)
691   end if
692 
693 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

PARENTS

      m_model_screening

CHILDREN

SOURCE

815 subroutine print_peaks(omega,fval,nomega,nfreqre,nfreqim)
816 
817 
818 !This section has been created automatically by the script Abilint (TD).
819 !Do not modify the following lines by hand.
820 #undef ABI_FUNC
821 #define ABI_FUNC 'print_peaks'
822 !End of the abilint section
823 
824  implicit none
825 
826 !Arguments ------------------------------------
827 !scalars
828   integer,intent(in)   :: nomega,nfreqre,nfreqim
829 !arrays
830   complex(dpc) ,intent(in)  :: omega(nomega)
831   complex(gwpc),intent(in)  :: fval(nomega)
832 
833 !Local variables-------------------------------
834 !scalars
835   integer :: ire,iim,unt_tmp
836   integer :: idx1,idx2,idx3
837   real(gwp) :: rez,imz,val1,val2,val3
838   character(len=500) :: msg
839 
840 ! *********************************************************************
841 
842   if (open_file("grid_peak_tree.dat", msg, newunit=unt_tmp) /= 0) then
843     MSG_ERROR(msg)
844   end if
845 
846   do iim=nfreqim,1,-1
847 !    write(std_out,*) ' iim:',iim
848 !   Check first points
849     idx1 = nfreqre+iim
850     idx2 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)
851 !    write(std_out,*) ' idx1:',idx1,' idx2:',idx2
852     val1 = AIMAG(fval(idx1))
853     val2 = AIMAG(fval(idx2))
854     if (ABS(val1)>ABS(val2)) then
855       rez = REAL(omega(idx1))
856       imz = AIMAG(omega(idx1))
857       write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1
858     end if
859 !   Do all but the last
860     do ire=1,nfreqre-4
861        idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+ire
862        idx2 = idx1+1
863        idx3 = idx1+2
864 !       write(std_out,*) ' idx1:',idx1,' idx2:',idx2,' idx3:',idx3
865        rez = REAL(omega(idx2))
866        imz = AIMAG(omega(idx2))
867        val1 = AIMAG(fval(idx1))
868        val2 = AIMAG(fval(idx2))
869        val3 = AIMAG(fval(idx3))
870        if (((val1<val2).AND.(val2>val3))) then
871          if (sign(1.0_gwp,val2)<zero) CYCLE
872          write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
873        else if (((val1>val2).AND.(val2<val3))) then
874          if (sign(1.0_gwp,val2)>zero) CYCLE
875          write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
876        end if
877     end do
878 !   Check last point
879     idx1 = nfreqre+nfreqim+1+(nfreqre-1)*(iim-1)+nfreqre-3
880     idx2 = idx1 + 1
881 !   write(std_out,*) ' idx1:',idx1,' idx2:',idx2
882     rez = REAL(omega(idx2))
883     imz = AIMAG(omega(idx2))
884     val1 = AIMAG(fval(idx1))
885     val2 = AIMAG(fval(idx2))
886     if (ABS(val1)<ABS(val2)) then
887       write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
888     end if
889   end do
890 ! finally, do the purely real axis
891 ! Check first points
892   idx1 = 1; idx2 = 2
893   val1 = AIMAG(fval(idx1)); val2 = AIMAG(fval(idx2))
894   if (ABS(val1)>ABS(val2)) then
895     rez = REAL(omega(idx1))
896     imz = AIMAG(omega(idx1))
897     write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val1
898   end if
899   do ire=2,nfreqre-3
900     idx1 = ire; idx2 = idx1+1; idx3 = idx1+2
901     rez = REAL(omega(idx2))
902     imz = AIMAG(omega(idx2))
903     val1 = AIMAG(fval(idx1))
904     val2 = AIMAG(fval(idx2))
905     val3 = AIMAG(fval(idx3))
906     if (((val1<val2).AND.(val2>val3))) then
907       if (sign(1.0_gwp,val2)<zero) CYCLE
908       write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
909     else if (((val1>val2).AND.(val2<val3))) then
910       if (sign(1.0_gwp,val2)>zero) CYCLE
911       write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
912     end if
913   end do
914 ! Check last point
915   idx1 = nfreqre-2
916   idx2 = idx1 + 1
917 !  write(std_out,*) ' idx1:',idx1,' idx2:',idx2
918   rez = REAL(omega(idx2))
919   imz = AIMAG(omega(idx2))
920   val1 = AIMAG(fval(idx1))
921   val2 = AIMAG(fval(idx2))
922   if (ABS(val1)<ABS(val2)) then
923     write(unt_tmp,'(2f8.2,4x,ES16.8)') rez*Ha_eV,imz*Ha_eV,val2
924   end if
925 
926   close(unt_tmp)
927 
928 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

PARENTS

      m_model_screening

CHILDREN

SOURCE

248 subroutine re_and_im_screening(omega,fval,nomega,coeff,ncoeff)
249 
250 
251 !This section has been created automatically by the script Abilint (TD).
252 !Do not modify the following lines by hand.
253 #undef ABI_FUNC
254 #define ABI_FUNC 're_and_im_screening'
255 !End of the abilint section
256 
257   implicit none
258 
259 !Arguments ------------------------------------
260 !scalars
261   integer,intent(in)   :: nomega,ncoeff
262 !arrays
263   complex(dpc) ,intent(in)  :: omega(nomega)
264   real(gwp)   ,intent(in)  :: coeff(ncoeff)
265   complex(gwp),intent(out) :: fval(nomega)
266 
267 !Local variables-------------------------------
268 !scalars
269   integer :: io,ip
270   real(gwp) :: rez,imz,realp,imagp,refval,imfval
271   real(gwp) :: fn,wn,gamman
272 ! *********************************************************************
273 
274 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
275 !    /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
276 
277   do io=1,nomega
278     fval(io) = 0.0
279     rez =  REAL(omega(io))
280     imz = AIMAG(omega(io))
281     do ip=1,ncoeff,3
282       fn     = coeff(ip)
283       wn     = coeff(ip+1)
284       gamman = coeff(ip+2)
285       realp  = -imz*gamman+wn*wn+imz*imz-rez*rez
286       imagp  = rez*(two*imz-gamman)
287 
288         refval   = fn*realp/((realp*realp)+(imagp*imagp))
289         imfval   = fn*imagp/((realp*realp)+(imagp*imagp))
290 
291         fval(io) = fval(io)-CMPLX(refval,imfval)
292 
293     end do
294   end do
295 
296 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

PARENTS

      calc_sigc_pole_cd

CHILDREN

SOURCE

331 subroutine re_and_im_screening_with_phase(omega,fval,nomega,coeff,ncoeff)
332 
333 
334 !This section has been created automatically by the script Abilint (TD).
335 !Do not modify the following lines by hand.
336 #undef ABI_FUNC
337 #define ABI_FUNC 're_and_im_screening_with_phase'
338 !End of the abilint section
339 
340   implicit none
341 
342 !Arguments ------------------------------------
343 !scalars
344   integer,intent(in)   :: nomega,ncoeff
345 !arrays
346   complex(dpc) ,intent(in)  :: omega(nomega)
347   real(gwp)   ,intent(in)  :: coeff(ncoeff)
348   complex(gwpc),intent(out) :: fval(nomega)
349 
350 !Local variables-------------------------------
351 !scalars
352   integer :: io,ip,npoles
353   real(gwp) :: rez,imz,realp,imagp,refval,imfval,retemp,imtemp
354   real(gwp) :: fn,wn,gamman,imrot,rerot
355 ! *********************************************************************
356 
357 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
358 !    /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
359   npoles = (ncoeff-1)/3
360 
361   do io=1,nomega
362     fval(io) = 0.0
363     rez =  REAL(omega(io))
364     imz = AIMAG(omega(io))
365     do ip=1,(ncoeff-1),3
366       fn     = coeff(ip)
367       wn     = coeff(ip+1)
368       gamman = coeff(ip+2)
369       realp  = -imz*gamman+wn*wn+imz*imz-rez*rez
370       imagp  = rez*(two*imz-gamman)
371 
372         refval   = fn*realp/((realp*realp)+(imagp*imagp))
373         imfval   = fn*imagp/((realp*realp)+(imagp*imagp))
374 
375         fval(io) = fval(io)-CMPLX(refval,imfval)
376 
377     end do
378     ! Restore phase
379     rerot = COS(coeff(npoles*3+1))
380     imrot = SIN(coeff(npoles*3+1))
381     retemp = REAL(fval(io))
382     imtemp = AIMAG(fval(io))
383     fval(io) = CMPLX(rerot*retemp-imrot*imtemp,rerot*imtemp + imrot*retemp)
384   end do
385 
386 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

PARENTS

CHILDREN

SOURCE

168 subroutine re_screening(omega,fval,nomega,coeff,ncoeff)
169 
170 
171 !This section has been created automatically by the script Abilint (TD).
172 !Do not modify the following lines by hand.
173 #undef ABI_FUNC
174 #define ABI_FUNC 're_screening'
175 !End of the abilint section
176 
177   implicit none
178 
179 !Arguments ------------------------------------
180 !scalars
181   integer,intent(in)   :: nomega,ncoeff
182 !arrays
183   complex(dpc),intent(in)  :: omega(nomega)
184   real(gwp)  ,intent(in)  :: coeff(ncoeff)
185   real(gwp)  ,intent(out) :: fval(nomega)
186 
187 !Local variables-------------------------------
188 !scalars
189   integer :: io,ip
190   real(gwp) :: rez,imz,realp,imagp
191   real(gwp) :: fn,wn,gamman
192 ! *********************************************************************
193 
194 ! The expression is: fn*(-imz*gamma+w_n^2+imz^2-rez^2)
195 !    /( (-imz*gamma+w_n^2+imz^2-rez^2)^2 + (2*rez*imz-rez*gamma)^2 )
196 
197   do io=1,nomega
198     fval(io) = 0.0
199     rez =  REAL(omega(io))
200     imz = AIMAG(omega(io))
201     do ip=1,ncoeff,3
202       fn     = coeff(ip)
203       wn     = coeff(ip+1)
204       gamman = coeff(ip+2)
205       realp  = -imz*gamman+wn*wn+imz*imz-rez*rez
206       imagp  = rez*(two*imz-gamman)
207 
208         fval(io) = fval(io)-fn*realp/((realp*realp)+(imagp*imagp))
209 
210     end do
211   end do
212 
213 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

PARENTS

      mrgscr

CHILDREN

SOURCE

1145 subroutine remove_phase(fval,nomega,phase)
1146 
1147 
1148 !This section has been created automatically by the script Abilint (TD).
1149 !Do not modify the following lines by hand.
1150 #undef ABI_FUNC
1151 #define ABI_FUNC 'remove_phase'
1152 !End of the abilint section
1153 
1154   implicit none
1155 
1156 !Arguments ------------------------------------
1157 !scalars
1158   integer,    intent(in)  :: nomega
1159   real(gwp), intent(out) :: phase
1160 !arrays
1161   complex(gwpc), intent(inout) :: fval(nomega)
1162 
1163 !Local variables-------------------------------
1164 !scalars
1165   integer       :: io
1166   real(gwp)    :: a,b,retemp,imtemp
1167 
1168 ! *********************************************************************
1169 
1170 ! The phase can be found by checking when the function is
1171 !  identically zero along the imaginary axis
1172   if (ABS(AIMAG(fval(1)))<tol14) then ! Phase is zero
1173     phase = zero
1174     RETURN
1175   else if (ABS(REAL(fval(1)))<tol14) then ! Phase is exactly pi/2
1176     phase = pi*half
1177     a = zero
1178     b = -1.0_gwp
1179   else
1180     phase = ATAN(AIMAG(fval(1))/REAL(fval(1)))
1181     a = COS(phase)
1182     b = -SIN(phase)
1183   end if
1184 ! Rotate values
1185   do io=1,nomega
1186     retemp = REAL(fval(io))
1187     imtemp = AIMAG(fval(io))
1188     fval(io) = CMPLX(a*retemp-b*imtemp,a*imtemp+b*retemp)
1189   end do
1190 
1191 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

PARENTS

CHILDREN

SOURCE

419 subroutine sequential_fitting(omega,refval,imfval,nomega,nfreqre,coeff,&
420 & ncoeff,prtvol,startcoeff)
421 
422 
423 !This section has been created automatically by the script Abilint (TD).
424 !Do not modify the following lines by hand.
425 #undef ABI_FUNC
426 #define ABI_FUNC 'sequential_fitting'
427 !End of the abilint section
428 
429   implicit none
430 
431 !Arguments ------------------------------------
432 !scalars
433   integer,intent(in)   :: nomega,nfreqre,ncoeff,prtvol
434 !arrays
435   complex(dpc),intent(in)     :: omega(nomega)
436   real(gwp)  ,intent(out)    :: coeff(ncoeff)
437   real(gwp)  ,intent(inout)  :: refval(nomega),imfval(nomega)
438   real(gwp),optional,intent(out) :: startcoeff(ncoeff)
439 
440 !Local variables-------------------------------
441 !scalars
442   integer :: ip,npoles,idx
443   real(gwp) :: thiscoeff(3),norm,invnorm
444   real(dpc)  :: re_zvals(nomega),im_zvals(nomega)
445 !  real(dpc)  :: orig_refval(nomega),orig_imfval(nomega)
446   complex(gwp) :: pole_func(nomega)
447 ! *********************************************************************
448 
449   npoles = ncoeff/3
450   re_zvals(:) = REAL(omega(:))
451   im_zvals(:) = AIMAG(omega(:))
452 
453   ! Normalise
454   norm    = MAXVAL(ABS(imfval))
455   invnorm = 1.0_gwp/norm
456   refval  = invnorm*refval
457   imfval  = invnorm*imfval
458 
459   ! Loop over poles to fit
460   do ip=1,npoles
461     idx = 3*(ip-1)+1
462     ! Initialise pole
463     call init_single_peak(omega,refval,imfval,nomega,nfreqre,thiscoeff,prtvol)
464     if (present(startcoeff)) then
465       startcoeff(idx:idx+2) = thiscoeff(1:3)
466     end if
467     ! Make fit
468 #ifdef HAVE_LEVMAR
469     call dfit_re_and_im_screening(re_zvals,im_zvals,imfval,refval,&
470 &    nomega,3,thiscoeff,prtvol)
471 #else
472     MSG_ERROR(' ABINIT was not compiled with the levmar library!')
473 #endif
474     ! Remove current fit
475     call re_and_im_screening(omega,pole_func,nomega,thiscoeff,3)
476     refval(:) = refval(:) -  REAL(pole_func(:))
477     imfval(:) = imfval(:) - AIMAG(pole_func(:))
478     coeff(idx:idx+2) = thiscoeff(1:3)
479     coeff(idx) = norm*coeff(idx)
480   end do
481 
482 end subroutine sequential_fitting