fig

figalt

figalt

figalt

figalt

figalt

figalt

code

[[file:./2019-12-19-physics-dwaveRF/p-WaveSinglePole.f90][fortran code download]]

#+BEGIN_SRC fortran ! ========================== ! program : pWaveSinglePole ! version : 02-Dec-2019 ! author : Zeqing Wang ! purpose : calculate the RF spectrum of p-wave sue single pole approximation ! reference: ! ==========================

! ========================== ! program : param ! type : module ! version : 05-Dec-2019 ! purpose : set the precision ! comment : ! ========================== module prec implicit none integer, parameter :: rkind=16 ! the precision end module prec ! ==========================

! ========================== ! program : param ! type : module ! version : 05-Dec-2019 ! purpose : set the precision ! comment : ! ========================== module param use prec implicit none real(kindrkind), parameter :: piacos(-1.0_rkind) real(kindrkind), parameter :: beta1_rkind, mu-1_rkind, rkv& &-1.0_rkind, rp=1.0/30.0 real(kindrkind), parameter :: phi_k0.0 ! phi_k 是无关的量, 不影响结果 end module param ! ==========================

! ========================== ! program : pWavesinglePole ! type : program ! version : 19-Dec-2019 ! comment : the main program ! ========================== program pWaveSinglePole use prec use param implicit none

integer, parameter :: dim=50 ! 画图取点的个数 integer :: i ! 循环变量 real(kind=rkind) :: x(dim) ,y(dim) ! 纵坐标 real(kind=rkind) :: sigma, spectrum_k_thetak ! 画图用到的函数 real(kindrkind) :: k2.0, theta_k = pi/2 ! 给定 k 和 theta_k 画图 real(kind=rkind) :: start, end ! 用于计时的变量

interface function linspace(a, b, dim) result(list) use prec integer :: dim real(kind=rkind) :: a, b, list(dim) end function linspace end interface

! 计时开始 call cpu_time(start) x = linspace(-20.0_rkind, 10.0_rkind, dim) ! x = linspace(0.0_rkind, pi, dim) do i = 1, dim y(i) = sigma(x(i), k, theta_k) ! y(i) = spectrum_k_thetak(x(i), k, theta_k) print *, 'calculating...', i, 'of', dim end do

! 计时结束 call cpu_time(end) print *, 'calculate finish! time is ', end - start, 'seconds'

! 将数据保存到文件, 用于画图 open(1, file='./data/datax.csv') write(1, *) x close(1) open(2, file='./data/datay.csv') write(2, *) y close(2) print *, 'save data finish!' end program pWaveSinglePole

! ==========================

! ========================== ! program : spectrum ! type : function ! version : 19-Dec-2019 ! purpose : calculate spectrum function ! comment : 还是只有虚部 ! ========================== function spectrum_k_thetak(omega, k, theta_k) use prec use param implicit none

real(kind=rkind) :: omega, k, theta_k real(kind=rkind) :: spectrum_k_thetak real(kind=rkind) :: sigma, sigma_real, sigma_imag sigma_imag = sigma(omega, k, theta_k) ! print *, sigma_imag sigma_real = 0.0_rkind

spectrum_k_thetak = -2*sigma_imag / ((omega - k**2 - sigma_real)**2 +& & sigma_imag**2) end function spectrum_k_thetak ! ==========================

! ========================== ! program : sigma ! type : function ! version : 18-Dec-2019 ! purpose : calculate self energy ! comment : 积掉 phi_q . 还是只有虚部 ! ========================== function sigma(omega, k, theta_k) use prec use param implicit none

real(kind=rkind) :: sigma

real(kind=rkind) :: omega, k, theta_k, theta_q, phi_q

real(kind=rkind) :: gaussQuad, r1, r2, x, cos_theta_kq, deltapart integer, parameter :: n=50

! 将没有积掉 phi_q 时的自能接口进来 interface function sigma_phiq(omega, k, theta_k, phi_q) use prec implicit none real(kind=rkind) :: omega, k, theta_k, phi_q real(kind=rkind) :: sigma_phiq end function sigma_phiq end interface sigma = gaussQuad(fun, 0.0_rkind, 50.0_rkind, n)

contains function fun(phi_q) use prec implicit none real(kind=rkind) :: phi_q, fun fun = sigma_phiq(omega, k, theta_k, phi_q) end function fun end function sigma ! ==========================

! ========================== ! program : sigma_phiq ! type : function ! version : 18-Dec-2019 ! purpose : calculate self energy ! comment : 这里把 q 积掉了, 还乘下 theta_q, phi_q 没有积 ! 再积掉 theta_q, 还乘下 phi_q 没有积 ! ========================== function sigma_phiq(omega, k, theta_k, phi_q) use prec use param implicit none

real(kind=rkind) :: omega, k, theta_k, theta_q, phi_q real(kind=rkind) :: sigma_phiq, sp real(kind=rkind) :: gaussQuad, r1, r2, x, cos_theta_kq, deltapart integer, parameter :: n=50

! 将没有积掉 theta_q 时的自能接口进来 interface function sigma_thetaq_phiq(omega, k, theta_k, theta_q, phi_q) use prec implicit none real(kind=rkind) :: omega, k, theta_k, theta_q, phi_q real(kind=rkind) :: sigma_thetaq_phiq end function sigma_thetaq_phiq end interface sigma_phiq = gaussQuad(fun, 0.0_rkind, 50.0_rkind, n)

! 将没有积掉 theta_q 的自能定义成一个 theta_q 的函数 contains function fun(theta_q) use prec implicit none real(kind=rkind) :: theta_q, fun fun = sigma_thetaq_phiq(omega, k, theta_k, theta_q, phi_q) end function fun end function sigma_phiq ! ==========================

! ========================== ! program : sigma_thetaq_phiq ! type : function ! version : 18-Dec-2019 ! purpose : calculate self energy ! comment : 这里把 q 积掉了, 还乘下 theta_q, phi_q 没有积 ! ========================== function sigma_thetaq_phiq(omega, k, theta_k, theta_q, phi_q) use prec use param implicit none

real(kind=rkind) :: omega, k, theta_k, theta_q, phi_q real(kind=rkind) :: sigma_thetaq_phiq, sp real(kind=rkind) :: cauthyGaussQuad, r1, r2, x, cos_theta_kq, deltapart real(kind=rkind) :: delta ! b^2 - 4*a*c

! 将没有积掉 q 时的自能接口进来 interface function numerator_of_sigma(omega, k, q, theta_k, theta_q, phi_q) use prec implicit none real(kind=rkind) :: omega, k, q, theta_k, theta_q, phi_q real(kind=rkind) :: numerator_of_sigma end function numerator_of_sigma end interface x = cos_theta_kq(theta_k, theta_q, phi_k, phi_q) ! x=0 delta = 4*k**2*x**2 - 2*(omega + k**2 +mu + rkv) if (delta .gt. 0.0_rkind) then r1 = 2*k*x - sqrt(delta) r2 = 2*k*x + sqrt(delta) sigma_thetaq_phiq = deltapart(fun, r1, r2) ! 计算出虚部 else sigma_thetaq_phiq = 0 ! 计算出虚部 end if

! 将没有积掉 q 的自能定义成一个 q 的函数(便于计算虚部, 去掉了分母) contains function fun(q) use prec implicit none real(kind=rkind) :: q, fun fun = numerator_of_sigma(omega, k, q, theta_k, theta_q, phi_q) end function fun end function sigma_thetaq_phiq ! ==========================

! ========================== ! program : sigma_q_thetaq_phiq ! type : function ! version : 05-Dec-2019 ! purpose : calculate self energy ! comment : 这里有三个待积变量, q, theta_q, phi_q ! ========================== function sigma_q_thetaq_phiq(omega, k, q, theta_k, theta_q, phi_q) use prec use param implicit none

real(kind=rkind) :: omega, k, theta_k, q, theta_q, phi_q ! real(kindrkind) :: phi_k0_rkind real(kind=rkind) :: real, imag real(kind=rkind) :: sigma_q_thetaq_phiq real(kind=rkind) :: kp2 !square of k' real(kind=rkind) :: cos_theta_kq, bose, cos_theta_kprime, sylm10

kp2 = 4*k**2 + q**2 - 4*k*q*cos_theta_kq(theta_k, theta_q, phi_k, phi_q) !k'^2

sigma_q_thetaq_phiq = bose(k**2 + q**2 - 2*cos_theta_kq(theta_k, theta_q, phi_k& &, phi_q) - mu, beta)& & - bose(q**2/2 - 2*mu - rkv, beta) sigma_q_thetaq_phiq = sigma_q_thetaq_phiq / (omega + k**2 + q**2/2 -2*k*q& &*cos_theta_kq(theta_k, theta_q, phi_k, phi_q) + mu +& & rkv)

q**2 * kp2 * sin(theta_q) *&

& sylm10(cos_theta_kprime(k, q, cos_theta_kq(theta_k, theta_q, phi_k,&

rp * 2 / (pi**2)

end function sigma_q_thetaq_phiq ! ==========================

! ========================== ! program : numerator_of_sigma ! type : function ! version : 24-Dec-2019 ! purpose : calculate self energy ! comment : 这里有三个待积变量, q, theta_q, phi_q ! function sigma_q_thetaq_phiq 去掉分母, 用于计算虚部 ! ========================== function numerator_of_sigma(omega, k, q, theta_k, theta_q, phi_q) use prec use param implicit none

real(kind=rkind) :: omega, k, theta_k, q, theta_q, phi_q ! real(kindrkind) :: phi_k0_rkind real(kind=rkind) :: real, imag real(kind=rkind) :: numerator_of_sigma real(kind=rkind) :: kp2 !square of k' real(kind=rkind) :: cos_theta_kq, bose, cos_theta_kprime, sylm10

kp2 = k**2 + q**2/4 - k*q*cos_theta_kq(theta_k, theta_q, phi_k, phi_q) !k'^2

numerator_of_sigma = bose(k**2 + q**2 - 2*cos_theta_kq(theta_k, theta_q, phi_k& &, phi_q) - mu, beta)& & - bose(q**2/2 - 2*mu - rkv, beta)

q**2 * kp2 * sin(theta_q) *&

& sylm10(cos_theta_kprime(k, q, cos_theta_kq(theta_k, theta_q, phi_k,&

rp * 2 / (pi**2)

end function numerator_of_sigma ! ==========================

! ========================== ! program : cos_theta_kprime ! type : function ! version : 05-Dec-2019 ! purpose : known k, q, and the angle between k and q, calculate cos(k'), ! where k' is the angle between k and -k+q ! comment : ! ========================== function cos_theta_kprime(k, q, cos_theta_kq) use prec implicit none

real(kind=rkind), intent(in) :: k, q, cos_theta_kq real(kind=rkind) :: cos_theta_kprime

cos_theta_kprime = (2*k*q*cos_theta_kq - 2*k**2) / (2*k*sqrt(k**2 + q**2 -& & 2*k*q*cos_theta_kq)) end function cos_theta_kprime ! ==========================

! ========================== ! program : cos_theta_kq ! type : function ! version : 05-Dec-2019 ! purpose : known theta_k, theta_q, phi_k, phi_q, calculate the cosine of ! angle between k and q ! comment : ! ========================== function cos_theta_kq(theta_k, theta_q, phi_k, phi_q) use prec implicit none

real(kind=rkind), intent(in) :: theta_k, theta_q, phi_k, phi_q real(kind=rkind) :: cos_theta_kq

cos_theta_kq = sin(theta_k)*sin(theta_q)*cos(phi_k - phi_q) + cos(theta_k)& &*cos(theta_q) end function cos_theta_kq ! ==========================

! ========================== ! program : bose ! type : function ! version : 05-Dec-2019 ! purpose : Bose-Einstein function ! comment : n(x, beta) = 1 / (e^(beta*x) - 1) ! beta: 1/Temperature x: energy ! ========================== function bose(energy, beta) use prec implicit none

real(kind=rkind), intent(in) :: energy, beta real(kind=rkind) :: bose if (energy .lt. 0.0_rkind) then print *, '==Energy is negtive, Unphysical!===' stop end if

bose = exp(-beta*energy) / (1 - exp(- beta*energy)) end function bose ! ==========================

! ========================== ! program : sylm10 ! type : function ! version : 05-Dec-2019 ! purpose : square of spherical harmonica function ! comment : |Y_l1 m0(x)|^2 ! ========================== function sylm10(x) use prec use param implicit none

real(kind=rkind), intent(in) :: x real(kind=rkind) :: sylm10

cos(x)**2

end function sylm10 ! ==========================

! ========================== ! program : cauthyGaussQuad ! type : function ! version : 05-Dec-2019 ! purpose : calculate the cauthy principal value integral of function "fun" ! from a to b with singlarity sp ! comment : ! ========================== function cauthyGaussQuad(fun, a, b, n, sp) use prec implicit none

integer, intent(in) :: n real(kind=rkind), intent(in) :: a, b, sp

! 被积函数接口 interface function fun(x) use prec implicit none real(kind=rkind) :: x, fun end function fun end interface integer :: j real(kind=rkind) :: cauthyGaussQuad, diffa, diffb, ff, sgq1, sgq2, gaussQuad

diffa = sp - a diffb = b - sp

! 计算积分 if (diffa .lt. diffb) then sgq1 = gaussQuad(ffun, 0.0_rkind, sp-a, n) sgq2 = gaussQuad(fun, 2*sp-a, b, n) else sgq1 = gaussQuad(ffun, 0.0_rkind, b-sp, n) sgq2 = gaussQuad(fun, a, 2*sp-b, n) end if

cauthyGaussQuad = sgq1 + sgq2 contains function ffun(t) use prec implicit none real(kind=rkind) :: t, ffun ffun = fun(t+sp) + fun(-t+sp) end function ffun end function cauthyGaussQuad ! ==========================

! ========================== ! program : gaussQuad ! type : function ! version : 03-Dec-2019 ! purpose : calculate the integral of function "fun" from a to b ! comment : ! ========================== function gaussQuad(fun, a, b, n) use prec implicit none

integer, intent(in) :: n real(kind=rkind), intent(in) :: a, b

! 被积函数接口 interface function fun(x) use prec implicit none real(kind=rkind) :: x, fun end function fun end interface integer :: j real(kind=rkind) :: gaussQuad, r(2, n), fxi(n)

! 根与权重接口 interface function node_weight(n) result(r) use prec implicit none integer :: n real(kind=rkind) :: r(2, n) end function node_weight end interface

! 计算积分 r = node_weight(n) do j = 1, n fxi(j) = fun((r(1, j)*(b - a) + a + b) / 2) end do gaussQuad = dot_product(r(2, :), fxi)

(b - a) / 2

end function gaussQuad ! ==========================

! ========================== ! program : node_weight ! type : function ! version : 02-Dec-2019 ! purpose : calculate the root of n-order Legendre polynomial, and weight ! comment : use the method in reference ! ========================== function node_weight(n) result(r) use prec implicit none integer :: n, i, iter, k real(kind=rkind) :: r(2, n), x, f, df, dx real(kindrkind), parameter :: pi acos(-1._rkind) real(kind=rkind), allocatable :: p0(:), p1(:) real(kind=rkind), allocatable :: tmp(:)

! 利用递推公式求 n 阶 Legendre 多项式的系数, 幂次从高到低排列, 结果就是数组 p1 p0 = [1.0_rkind] p1 = [1._rkind, 0._rkind] do i = 2, n tmp = ((2*i - 1)*[p1, 0._rkind] - (i - 1)*[0._rkind, 0._rkind, p0]) / i p0 p1; p1 tmp end do ! 这个函数将 n 阶的情况的 根存在 r(1,:) 中, 权重存在 r(2, :) 中 do i = 1, n x = cos(pi*(i - .25_rkind) / (n + .5_rkind)) do iter = 1, 10 f p1(1); df 0._rkind do k = 2, size(p1) df = f + x*df ! 得到的是 P_n'(x_0) 的值 f = p1(k) + x*f ! 得到的是 P_n(x_0) 的值 end do dx = f/df x = x - dx if (abs(dx)<10*epsilon(dx)) exit end do r(1, i) = x r(2, i) = 2/((1 - x**2)*df**2) end do end function node_weight ! ==========================

! ========================== ! program : linspace ! type : function ! version : 19-Dec-2019 ! purpose : similar to 'numpy.linspace' in python ! comment : ! ========================== function linspace(a, b, dim) result(list) use prec implicit none

real(kind=rkind) :: a, b integer :: dim real(kind=rkind) :: list(dim) integer :: i real(kind=rkind) :: diff

do i = 1, dim list(i) = i - 1 ! i-1 instead of i, in order to be the same as ! numpy.linspace end do

diff = b - a

diff/dim

end function linspace ! ==========================

! ========================== ! program : deltapart ! type : function ! version : 18-Dec-2019 ! purpose : get the delta part of \int_0 ^\inf f(x) / ((x-r1)(x-r2) + i0) dx

(f(r1) - f(r2)) / (r1 - r2)

! 注意, 这个程序没有得到验证! ! ========================== function deltapart(fun, r1, r2) use prec use param implicit none real(kind=rkind) :: r1, r2, deltapart ! 函数 f(x) 的接口 interface function fun(x) use prec implicit none real(kind=rkind) :: x, fun end function fun end interface

if ((r1 .lt. 0) .and. (r2 .lt. 0)) then ! 两个根都小于零, 没有虚部 deltapart = 0 ! print *, 'no root' else if (r1 .lt. 0) then ! r1<0, r1>0, 取 r2

fun(r2) / (r1-r2)

! print *, 'root2' ! print *, fun(r2) else if (r2 .lt. 0) then ! r1<0, r2>0, 取 r1

fun(r1) / (r1-r2)

! print *, 'root1' else ! 两个根都大于0

(fun(r1) - fun(r2)) / (r1-r2)

! print *, '2 roots' end if end function deltapart ! ========================== #+END_SRC