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