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