TABLE OF CONTENTS


ABINIT/nfourier [ Functions ]

[ Top ] [ Functions ]

NAME

  nfourier

FUNCTION

COPYRIGHT

  Copyright (C) 2014-2018 ABINIT group (XG)
  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

 20 #if defined HAVE_CONFIG_H
 21 #include "config.h"
 22 #endif
 23 
 24 #include "abi_common.h"
 25 
 26 MODULE m_lib_four
 27 
 28  use defs_basis
 29  use m_errors
 30  use m_abicore
 31 
 32 contains
 33 
 34 ! This routine contains direct and inverse fourier transformation
 35 ! It is a modification of a routine of the GNU GPL
 36 ! code available on http://dmft.rutgers.edu/ and
 37 ! described in the RMP 2006 paper written by
 38 ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
 39 !=======+=========+=========+=========+=========+=========+=========+=$
 40 !       TYPE   : SUBROUTINE
 41 !       PROGRAM: nfourier
 42 !       PURPOSE: fourier-transform the natural-spline interpolation
 43 !                of function Green(tau)
 44 !                calculate function Green(omega)
 45 !       I/O    :
 46 !       VERSION: 2-16-92
 47 !                29-Nov-95 removal of minimal bug concerning
 48 !                          DIMENSION of rindata
 49 !       COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
 50 !                analysis (Springer, New York, 1980)
 51 !=======+=========+=========+=========+=========+=========+=========+=$
 52 !
 53       SUBROUTINE nfourier(rindata,coutdata,iflag,Iwmax,L,Beta)
 54 
 55 !include 'param.dat'
 56 
 57 !This section has been created automatically by the script Abilint (TD).
 58 !Do not modify the following lines by hand.
 59 #undef ABI_FUNC
 60 #define ABI_FUNC 'nfourier'
 61 !End of the abilint section
 62 
 63        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 64        IMPLICIT INTEGER(I-N)
 65        integer  Iwmax,L,iflag
 66        DIMENSION rindata(L)
 67        DIMENSION rincopy(L+1),a(L),b(L),c(L),d(L),u(L+1), q(L+1),XM(L+1)
 68        complex*16 :: coutdata(Iwmax+1)
 69        complex*16 cdummy,explus,ex
 70        xpi = ACOS(-One)
 71        delta = Beta/L
 72        DO i = 1,L
 73           rincopy(i) = rindata(i)
 74        ENDDO
 75        if(iflag==1) then
 76          rincopy(L+1) = -1-rindata(1)
 77        else
 78          rincopy(L+1) = -rindata(1)
 79        endif
 80 !       Three = Two+One
 81 !       six = Two*Three
 82 
 83 !c
 84 !c     spline interpolation:  the spline is given by
 85 !c     G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
 86 !c     The following formulas are taken directly from  Stoer and
 87 !c     Bulirsch p. 102
 88 !c
 89        q(1) = Zero
 90        u(1) = Zero
 91        DO k = 2,L
 92           p = q(k-1)/Two+Two
 93           q(k)=-One/Two/p
 94           u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
 95           u(k)=(u(k)-u(k-1)/Two)/p
 96        ENDDO
 97        XM(L+1) = 0
 98        DO k = L,1,-1
 99           XM(k) = q(k)*XM(k+1)+u(k)
100        ENDDO
101 !c
102 !c     The following formulas are taken directly from  Stoer and
103 !c     Bulirsch p. 98
104 !c
105        DO j = 1, L
106           a(j) = rincopy(j)
107           c(j) = XM(j)/Two
108           b(j) = (rincopy(j+1)-rincopy(j))/delta - &
109      &       (Two*XM(j)+XM(j+1))*delta/6.
110           d(j) = (XM(j+1)-XM(j))/(6.*delta)
111        ENDDO
112 
113 !c
114 !c     The Spline multiplied by the exponential can now be exlicitely
115 !c     integrated. The following formulas were obtained using
116 !c     MATHEMATICA
117 !c
118         DO i = 0,Iwmax
119            om = (Two*(i)+One)*xpi/Beta
120            coutdata(i+1) = czero
121            DO j = 1,L
122               cdummy = j_dpc*om*delta*j
123               explus = exp(cdummy)
124               cdummy = j_dpc*om*delta*(j-1)
125               ex = exp(cdummy)
126               coutdata(i+1) = coutdata(i+1) + explus*(&
127      &         ( -six* d(j) )/om**4 + &
128      &         ( Two*j_dpc*c(j) + six*delta*j_dpc*d(j)  )/om**3 +&
129      &         ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +&
130      &         (- j_dpc*a(j) - delta*j_dpc*b(j) - delta**2*j_dpc*c(j) -&
131      &         delta**3*j_dpc*d(j))/om)
132 
133               coutdata(i+1) = coutdata(i+1) + ex*(&
134      &        six*d(j)/om**4 - Two*j_dpc*c(j)/om**3 &
135      &        -b(j)/om**2 + j_dpc*a(j)/om)
136            ENDDO
137         ENDDO
138         end subroutine nfourier