Skip to content

Instantly share code, notes, and snippets.

@DSCF-1224

DSCF-1224/README.md

Last active January 6, 2024 05:02
Show Gist options
  • Save DSCF-1224/3fadaac4e22efadb5c3b48e9c55e6e4d to your computer and use it in GitHub Desktop.
Save DSCF-1224/3fadaac4e22efadb5c3b48e9c55e6e4d to your computer and use it in GitHub Desktop.
Download ZIP
vmc の fortran 2008 実装
Raw
README.md

参考文献

実行例

Rev1

$ time build/test.exe
  0.99472140592013003       0.80000001192092896     
   1.0000596745980763       0.99472140592013003     
   1.0000573806829209        1.0000596745980763     
   1.0000570092412857        1.0000573806829209     
   1.0000570092412857      1.0 になってたらOK
 -0.50001904165407052      -0.5 になってたらOK

real    0m10.736s
user    0m10.736s
sys     0m0.000s

Rev2

  • 各部の SUBROUTINE 化
  • n_samples に関するループを n_samples_burn_in 前後で2分割し，if 分岐を削除
  • 正規分布に従う乱数の生成手法の変更
$ time build/test.exe
  0.98263620904471616       0.80000001192092896     
   1.0004793997956174       0.98263620904471616     
   1.0000004783215848        1.0004793997956174     
  0.99999999912233239        1.0000004783215848     
  0.99999999912233239      1.0 になってたらOK
 -0.50000000000001021      -0.5 になってたらOK

real    0m12.340s
user    0m12.340s
sys     0m0.001s
  • Rev1 より遅くなった （収束が遅い？）．

Rev3

  • 正規分布に従う乱数の生成手法の変更
    • Polar 法 → Box-Muller 法
    • random_number を必要毎に call → random_number 1回で予め必要な分だけ生成
$ time build/test.exe
  0.97654018253803376       0.80000000000000004     
  0.99963143237072738       0.97654018253803376     
   1.0000001321911660       0.99963143237072738     
   1.0000000002128571        1.0000001321911660     
   1.0000000002128571      1.0 になってたらOK
 -0.50000000000024758      -0.5 になってたらOK

real    0m11.511s
user    0m11.511s
sys     0m0.000s
Raw
main.f90
module prng
use, intrinsic :: iso_fortran_env, only: dp => real64
implicit none
private
public :: box_muller_method
public :: polar_method
real(dp), parameter, private :: pi = acos(-1.0_dp)
!! A `PARAMETER` for this `MODULE`
!! mathematical constant: \pi
real(dp), parameter, private :: pi2 = 2.0_dp * pi
!! A `PARAMETER` for this `MODULE`
!! mathematical constant: 2\pi
contains
subroutine box_muller_method(harvest)
real(dp), dimension(3), intent(out) :: harvest
!! A dummy variable for this `SUBROUTINE`
real(dp) :: u(4)
!! A variable for this `SUBROUTINE`
call random_number( u(:) )
u(1:2) = sqrt( -2.0_dp * log( u(1:2) ) )
u(3:4) = pi2 * u(3:4)
harvest(1) = u(1) * cos( u(3) )
harvest(2) = u(1) * sin( u(3) )
harvest(3) = u(2) * cos( u(4) )
end subroutine box_muller_method
subroutine polar_method(harvest)
real(dp), dimension(3), intent(out) :: harvest
!! A dummy variable for this `SUBROUTINE`
real(dp) :: v1 !! A variable for this `SUBROUTINE`
real(dp) :: v2 !! A variable for this `SUBROUTINE`
real(dp) :: w !! A variable for this `SUBROUTINE`
call polar_method_kernel(v1, v2, w)
harvest(1) = v1 * w
harvest(2) = v2 * w
call polar_method_kernel(v1, v2, w)
harvest(3) = v1 * w
end subroutine polar_method
subroutine polar_method_kernel(v1, v2, w)
real(dp), intent(out) :: v1 !! A dummy argument for this `SUBROUTINE`
real(dp), intent(out) :: v2 !! A dummy argument for this `SUBROUTINE`
real(dp), intent(out) :: w !! A dummy argument for this `SUBROUTINE`
real(dp) :: v !! A variable for this `SUBROUTINE`
do
call random_number(v1)
call random_number(v2)
v1 = 2.0_dp * v1 - 1.0_dp
v2 = 2.0_dp * v2 - 1.0_dp
v = v1 * v1 + v2 * v2
if ( (0.0_dp < v) .and. (v < 1.0_dp) ) then
w = sqrt( - 2.0_dp * log(v) / v )
return
end if
end do
end subroutine polar_method_kernel
end module prng
module vmc
use, intrinsic :: iso_fortran_env, only: dp => real64
use, non_intrinsic :: prng
implicit none
private
public :: metropolis_test
integer, parameter, private :: n_walkers = 400
!! A `PARAMETER` for this `MODULE`
integer, parameter, private :: n_samples_burn_in = 4000
!! A `PARAMETER` for this `MODULE`
integer, parameter, private :: n_samples = 30000
!! A `PARAMETER` for this `MODULE`
real(dp), parameter, private :: sigma = 0.4_dp
!! A `PARAMETER` for this `MODULE`
!! variance of the normal distribution
real(dp), parameter, private :: pi = acos(-1.0_dp)
!! A `PARAMETER` for this `MODULE`
!! mathematical constant: \pi
real(dp), parameter, private :: pi2 = 2.0_dp * pi
!! A `PARAMETER` for this `MODULE`
!! mathematical constant: 2\pi
contains
subroutine metropolis_test(alpha, dE_E)
real(dp), intent(in) :: alpha
!! A dummy argument for this `SUBROUTINE`
real(dp), dimension(2), intent(out) :: dE_E
!! A dummy argument for this `SUBROUTINE`
real(dp) :: d_ln_psi !! A variable for this `SUBROUTINE`
real(dp) :: E_loc !! A variable for this `SUBROUTINE`
real(dp) :: E_loc_d_ln_psi !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:) :: sum_d_ln_psi !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:) :: sum_E_loc !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:) :: sum_E_loc_d_ln_psi !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:) :: rr !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:,:) :: vec !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:) :: rr_next !! A variable for this `SUBROUTINE`
real(dp), allocatable, dimension(:,:) :: vec_next !! A variable for this `SUBROUTINE`
integer :: i_s !! A support variable for this `SUBROUTINE`
integer :: i_w !! A support variable for this `SUBROUTINE`
allocate( sum_d_ln_psi (n_walkers) , source=0.0_dp )
allocate( sum_E_loc (n_walkers) , source=0.0_dp )
allocate( sum_E_loc_d_ln_psi (n_walkers) , source=0.0_dp )
allocate( rr ( n_walkers) )
allocate( vec (3,n_walkers) )
allocate( rr_next ( n_walkers) )
allocate( vec_next (3,n_walkers) )
call random_number( vec(:,:) )
do concurrent (i_w = 1 : n_walkers)
vec (:,i_w) = 4.0_dp * ( vec(:,i_w) - 0.5_dp )
rr ( i_w) = sqrt( dot_product( vec(:,i_w), vec(:,i_w) ) )
end do
loop_samples_burn_in: &!
do i_s = 1, (n_samples_burn_in - 1)
do i_w = 1, n_walkers
call metropolis_test_kernel( alpha, vec(:,i_w), rr(i_w), vec_next(:,i_w), rr_next(i_w) )
end do
end do &!
loop_samples_burn_in
loop_samples_main: &!
do i_s = n_samples_burn_in, n_samples
do i_w = 1, n_walkers
call metropolis_test_kernel( alpha, vec(:,i_w), rr(i_w), vec_next(:,i_w), rr_next(i_w) )
E_loc = - 1.0_dp / rr(i_w) - 0.5_dp * alpha * ( alpha - 2.0_dp / rr(i_w) )
d_ln_psi = - rr(i_w)
E_loc_d_ln_psi = E_loc * d_ln_psi
sum_E_loc (i_w) = sum_E_loc (i_w) + E_loc
sum_d_ln_psi (i_w) = sum_d_ln_psi (i_w) + d_ln_psi
sum_E_loc_d_ln_psi (i_w) = sum_E_loc_d_ln_psi (i_w) + E_loc_d_ln_psi
end do
end do &!
loop_samples_main
associate( inv_denominator => 1.0_dp / real( n_samples - n_samples_burn_in + 1, dp ) )
sum_E_loc (:) = sum_E_loc (:) * inv_denominator
sum_d_ln_psi (:) = sum_d_ln_psi (:) * inv_denominator
sum_E_loc_d_ln_psi (:) = sum_E_loc_d_ln_psi (:) * inv_denominator
end associate
E_loc = sum( sum_E_loc (:) ) / n_walkers
d_ln_psi = sum( sum_d_ln_psi (:) ) / n_walkers
E_loc_d_ln_psi = sum( sum_E_loc_d_ln_psi (:) ) / n_walkers
dE_E(1) = 2.0 * (E_loc_d_ln_psi - E_loc * d_ln_psi)
dE_E(2) = E_loc
end subroutine metropolis_test
subroutine metropolis_test_kernel(alpha, vec, rr, vec_next, rr_next)
real(dp), intent(in) :: alpha
!! A dummy argument for this `SUBROUTINE`
real(dp), dimension(3), intent(inout) :: vec
!! A dummy argument for this `SUBROUTINE`
real(dp), intent(inout) :: rr
!! A dummy argument for this `SUBROUTINE`
real(dp), dimension(3), intent(out) :: vec_next
!! A dummy argument for this `SUBROUTINE`
real(dp), intent(out) :: rr_next
!! A dummy argument for this `SUBROUTINE`
real(dp) :: u
!! A support variable for this `SUBROUTINE`
real(dp), dimension(3) :: v
!! A support variable for this `SUBROUTINE`
call box_muller_method( v(:) )
! call polar_method( v(:) )
vec_next (:) = vec(:) + sigma * v(:)
rr_next = sqrt( dot_product( vec_next(:), vec_next(:) ) )
call random_number( u )
if ( u < exp( 2.0_dp * alpha * ( rr - rr_next ) ) ) then
vec (:) = vec_next (:)
rr = rr_next
end if
end subroutine metropolis_test_kernel
end module vmc
program main
use, intrinsic :: iso_fortran_env, only: dp=>real64
use, non_intrinsic :: vmc
implicit none
real(dp) :: alpha
real(dp) :: alpha_ini
real(dp) :: alpha_next
real(dp) :: alpha_plus
real(dp) :: alpha_minus
real(dp) :: E, dE, d_dE, dE_plus, dE_minus
real(dp) :: dE_E(2), dE_E_plus(2), dE_E_minus(2)
real(dp), parameter :: h = 1e-3
integer :: i
alpha_ini = 0.8_dp
alpha = alpha_ini
do i = 1, 10
call metropolis_test(alpha, dE_E)
dE = dE_E(1)
alpha_plus = alpha + h
call metropolis_test(alpha_plus, dE_E_plus)
dE_plus = dE_E_plus(1)
alpha_minus = alpha - h
call metropolis_test(alpha_minus, dE_E_minus)
dE_minus = dE_E_minus(1)
d_dE = (dE_plus - dE_minus) / (2.0_dp * h)
alpha_next = alpha - dE / d_dE
print *, alpha_next, alpha
if (abs(alpha_next - alpha) < 1e-6) then
alpha = alpha_next
exit
end if
alpha = alpha_next
end do
print *, alpha, "1.0 になってたらOK"
call metropolis_test(alpha, dE_E)
E = dE_E(2)
print *, E, "-0.5 になってたらOK"
end program main
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment