どこだ。。。どこにバグが。。。
!========== V_G potential ==========
subroutine VGpotential(dVG,il,dr)
include 'head.inc'
real(dp), dimension(mrmax,mrmax), intent(out) :: dVG
integer, intent(in) :: il
real(dp), dimension(mrmax), intent(in) :: dr
real(dp), parameter :: dsigma = 0.1_dp
integer, parameter :: ithe_max = 3000
real(dp), parameter :: dthestep = dpi/dble(ithe_max)
integer :: ir1, ir2
integer :: ithe
real(dp) :: dthe
real(dp) :: dsum
write(55,*) il
write(55,*) dr
dVG = 0._dp
do ir1 = 1, mrmax
do ir2 = 1, mrmax
dsum = 0._dp
do ithe = 1, ithe_max
dthe = dble(ithe) * dthestep
dsum = dsum&
+ cos(dble(il)*dthe) * exp(2d0*dr(ir1)*dr(ir2)*cos(dthe)/dsigma**2)
end do
dVG(ir1,ir2) = 2._dp * dg_G * exp(- (dr(ir1)**2+dr(ir2)**2)/dsigma**2)&
* dsum * dthestep
write(55,*) ir1, ir2, dVG(ir1,ir2)
end do
end do
end subroutine VGpotential
!---
mrmax = 140
il = 0, 1, 2, ...
dr(i) = 0, ..., 16
発見に一時間以上かかった。。。
答え
手計算で整理したときに、exp を二つに分けたんだけど、
実際の計算ではその指数が大きくなりすぎて、オーバーフロー
起こしてたみたい。NaN って形で計算結果がでたけど、
0/0 (ゼロ割るゼロ) 以外にオーバーフローでもその様に
計算結果を返すらしい。
exp の中身が大きくなり過ぎないように、do ループの外の
exp を do ループ内に入れて、足し引きを exp の肩でやれば
解決。