Skip to content

Instantly share code, notes, and snippets.

@terasakisatoshi

terasakisatoshi/main_vmc.f90

Last active January 5, 2024 05:08
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save terasakisatoshi/6baa636f926e4afccbfff9b2dffbcd03 to your computer and use it in GitHub Desktop.
Save terasakisatoshi/6baa636f926e4afccbfff9b2dffbcd03 to your computer and use it in GitHub Desktop.
Download ZIP
main_vmc.f90
Raw
main_vmc.f90
! gfortran -O3 -Wall -std=f2008 main_vmc.f90 && time ./a.out
module vmc
use, intrinsic :: iso_fortran_env, only: dp => real64
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 :: 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
real(dp), parameter, private :: sigma = 0.4_dp
!! A `PARAMETER` for this `MODULE`
!! variance of the normal distribution
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`
real(dp), dimension(5) :: u !! A support variable for this `SUBROUTINE`
real(dp), dimension(3) :: v !! A support variable for this `SUBROUTINE`
!allocate( sum_d_ln_psi (n_walkers), mold=0.0_dp )
!allocate( sum_E_loc (n_walkers), mold=0.0_dp )
!allocate( sum_E_loc_d_ln_psi (n_walkers), mold=0.0_dp )
allocate( sum_d_ln_psi (n_walkers))
allocate( sum_E_loc (n_walkers))
allocate( sum_E_loc_d_ln_psi (n_walkers))
do i_w = 1, n_walkers
sum_E_loc_d_ln_psi(i_w) = 0.0_dp
sum_E_loc(i_w) = 0.0_dp
sum_d_ln_psi(i_w) = 0.0_dp
end do
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_sampler: &!
do i_s = 1, n_samples
loop_walker: &!
do i_w = 1, n_walkers
call random_number( u(1:5) )
u(1) = sqrt( -2.0_dp * log( u(1) ) )
u(2) = pi2 * u(2)
v(1) = u(1) * cos( u(2) )
v(2) = u(1) * sin( u(2) )
v(3) = sqrt( -2.0_dp * log( u(3) ) ) * cos( pi2 * u(4) )
vec_next(:,i_w) = vec(:,i_w) + sigma * v(:)
rr_next ( i_w) = sqrt( dot_product( vec_next(:,i_w), vec_next(:,i_w) ) )
if ( u(5) < exp( 2.0_dp * alpha * ( rr(i_w) - rr_next(i_w) ) ) ) then
vec (:,i_w) = vec_next(:,i_w)
rr ( i_w) = rr_next ( i_w)
end if
end do &!
loop_walker
if (i_s >= n_samples_burn_in) then
do i_w = 1, n_walkers
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 if
end do &!
loop_sampler
associate( inv_denominator => 1.0_dp / real( n_samples - n_samples_burn_in + 1, dp ) )
do concurrent (i_w = 1 : n_walkers)
sum_E_loc (i_w) = sum_E_loc (i_w) * inv_denominator
sum_d_ln_psi (i_w) = sum_d_ln_psi (i_w) * inv_denominator
sum_E_loc_d_ln_psi (i_w) = sum_E_loc_d_ln_psi (i_w) * inv_denominator
end do
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
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
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
@terasakisatoshi
Copy link
Author

参考にしたコード

https://gist.github.com/DSCF-1224/3fadaac4e22efadb5c3b48e9c55e6e4d

@terasakisatoshi
Copy link
Author

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment