Skip to content

Draft: Resolve "measure acceptance of HMC moves as an observable" #540

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ALF-Import-Bot wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Draft: Resolve "measure acceptance of HMC moves as an observable" #540

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
89675eb
first shot at HMC acc as obser
johanneshofmann87 Sep 2, 2024
a447b46
fix call arguments
johanneshofmann87 Sep 3, 2024
6b3eaa9
implementing EJS
johanneshofmann87 Sep 25, 2024
3b38350
adding randomness to number of leapfrog steps, i.e., specitied number…
johanneshofmann87 Oct 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 19 additions & 0 deletions Prog/Hamiltonian_main_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,7 @@ Module Hamiltonian_main
procedure, nopass :: GR_reconstruction => GR_reconstruction_base
procedure, nopass :: GRT_reconstruction => GRT_reconstruction_base
procedure, nopass :: Apply_B_HMC => Apply_B_HMC_base
procedure, nopass :: Jump_dist_HMC => Jump_dist_HMC_base
#ifdef HDF5
procedure, nopass :: write_parameters_hdf5 => write_parameters_hdf5_base
#endif
Expand Down Expand Up @@ -694,6 +695,24 @@ Subroutine Apply_B_HMC_base(Forces_HMC, ltrans )
Logical , Intent(in) :: ltrans

end Subroutine Apply_B_HMC_base

!-------------------------------------------------------------------
real(kind=kind(0.d0)) function Jump_dist_HMC_base(fields_new, fields_old )

Implicit none

complex (Kind=Kind(0.d0)), Intent(in), allocatable :: fields_new(:,:), fields_old(:,:)

Integer :: I, J

Jump_dist_HMC_base=0.0d0
do I=1, size(fields_new,2)
do J=1, size(fields_new,1)
Jump_dist_HMC_base = Jump_dist_HMC_base + dble(fields_new(I,J) - fields_old(I,J))**2
enddo
enddo

end function Jump_dist_HMC_base

!--------------------------------------------------------------------
!> @brief
Expand Down
32 changes: 24 additions & 8 deletions Prog/Langevin_HMC_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
Module Langevin_HMC_mod

Use runtime_error_mod
use observables
Use Hamiltonian_main
Use UDV_State_mod
Use Control
Expand Down Expand Up @@ -70,6 +71,7 @@ Module Langevin_HMC_mod
Real (Kind=Kind(0.d0)), allocatable :: Det_vec_old(:,:)
Complex (Kind=Kind(0.d0)), allocatable :: Phase_Det_old(:)
Complex (Kind=Kind(0.d0)), allocatable :: Forces (:,:)
Type (Obser_Vec ), allocatable, public :: Obs_ESJD_HMC(:)

Real (Kind=Kind(0.d0)), allocatable :: Forces_0(:,:)
CONTAINS
Expand Down Expand Up @@ -352,7 +354,7 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab
Integer, intent(in) :: LOBS_ST, LOBS_EN, LTAU

!Local
Integer :: N_op, n, nt, n1, n2, i, j, t_leap, nf
Integer :: N_op, n, nt, n1, n2, i, j, t_leap, nf, LF_steps
Real (Kind=Kind(0.d0)) :: X, Xmax,E_kin_old, E_kin_new,T0_Proposal_ratio, weight, cluster_size
Logical :: Calc_Obser_eq, toggle
Real (Kind=Kind(0.d0)), allocatable :: Det_vec_old(:,:), Det_vec_new(:,:)
Expand Down Expand Up @@ -451,8 +453,9 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab
#endif

leap_frog_bulk=.true.
LF_steps=nranf(this%Leapfrog_Steps)
!Start Leapfrog loop (Leapfrog_steps)
Do t_leap=1,this%Leapfrog_Steps
Do t_leap=1,LF_steps
! update phi by delta t
this%Forces_0 = p_tilde
call ham%Apply_B_HMC(this%Forces_0 ,.True.)
Expand All @@ -462,11 +465,11 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab
#endif

! reset storage
if (t_leap == this%Leapfrog_Steps) leap_frog_bulk=.false.
if (t_leap == LF_steps) leap_frog_bulk=.false.
Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst)
! if last step
X=1.0d0
if (t_leap == this%Leapfrog_Steps) then
if (t_leap == LF_steps) then
! calc det[phi']
Call Compute_Fermion_Det(Phase_det_new,Det_Vec_new, udvl, udvst, Stab_nt, storage)
! LATER (optimization idea): save Phase, GR, udvr, udvl (move calc det out of the loop)
Expand All @@ -489,21 +492,21 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab
p_tilde=-p_tilde
!Do half step update of p
p_tilde=p_tilde - 0.5*this%Delta_t_Langevin_HMC*this%Forces_0
write(2000+this%Leapfrog_Steps,*) nsigma%f
write(2000+LF_steps,*) nsigma%f

!Start Leapfrog loop (Leapfrog_steps)
Do t_leap=1,this%Leapfrog_Steps
Do t_leap=1,LF_steps
! update phi by delta t
this%Forces_0 = p_tilde
call ham%Apply_B_HMC(this%Forces_0 ,.True.)
nsigma%f=nsigma%f + this%Delta_t_Langevin_HMC*this%Forces_0
write(2000+this%Leapfrog_Steps-t_leap,*) nsigma%f
write(2000+LF_steps-t_leap,*) nsigma%f

! reset storage
Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst)
! if last step
X=1.0d0
if (t_leap == this%Leapfrog_Steps) then
if (t_leap == LF_steps) then
! calc det[phi']
Call Compute_Fermion_Det(Phase_det_new,Det_Vec_new, udvl, udvst, Stab_nt, storage)
! LATER (optimization idea): save Phase, GR, udvr, udvl (move calc det out of the loop)
Expand Down Expand Up @@ -540,6 +543,11 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab
Phase_new = Phase_new*Phase_det_new(nf)
Enddo
Call Op_phase(Phase_new,OP_V,Nsigma,N_SUN)
X=ham%Jump_dist_HMC(nsigma%f, nsigma_old%f)
! write(*,*) LF_steps, x, x/size(nsigma%f), weight
this%Obs_ESJD_HMC(LF_steps)%N = this%Obs_ESJD_HMC(LF_steps)%N + 1
this%Obs_ESJD_HMC(LF_steps)%Ave_sign = this%Obs_ESJD_HMC(LF_steps)%Ave_sign + 1.d0
this%Obs_ESJD_HMC(LF_steps)%Obs_vec(1) = this%Obs_ESJD_HMC(LF_steps)%Obs_vec(1) + x*weight

TOGGLE = .false.
if ( Weight > ranf_wrap() ) Then
Expand Down Expand Up @@ -605,6 +613,8 @@ SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, Delta_t_Langevin_HMC, Max_Forc
!Local
Integer :: IERR
Logical :: lexist
Integer :: N
Character (len=64) :: Filename
#ifdef MPI
Real (Kind=Kind(0.d0)) :: X
INTEGER :: STATUS(MPI_STATUS_SIZE), ISIZE, IRANK
Expand Down Expand Up @@ -673,6 +683,12 @@ SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, Delta_t_Langevin_HMC, Max_Forc
this%L_Forces = .False.
this%Leapfrog_Steps = Leapfrog_steps
this%Delta_t_running = 1.0d0
allocate(this%Obs_ESJD_HMC(0:Leapfrog_steps))
N=1
do i=0,Leapfrog_steps
write(Filename,'(A,I0)') "ESJD_HMC_",i
Call Obser_Vec_make(this%Obs_ESJD_HMC(i),N,Filename)
enddo
! WRITE(error_unit,*) 'HMC step is not yet implemented'
! CALL Terminate_on_error(ERROR_GENERIC,__FILE__,__LINE__)
else
Expand Down
2 changes: 1 addition & 1 deletion Prog/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

OHAM = `./parse_ham.py --print_obj_list`

OBJS = Hamiltonians/LRC_mod.o control_mod.o Fields_mod.o Operator_mod.o WaveFunction_mod.o observables_mod.o \
OBJS = Hamiltonians/LRC_mod.o observables_mod.o control_mod.o Fields_mod.o Operator_mod.o WaveFunction_mod.o \
ContainerElementBase_mod.o DynamicMatrixArray_mod.o OpTTypes_mod.o OpT_time_dependent.o\
Predefined_Int_mod.o Predefined_Obs_mod.o Predefined_Latt_mod.o Predefined_Hop_mod.o Predefined_Trial_mod.o \
Hamiltonian_main_mod.o QDRP_decompose_mod.o udv_state_mod.o Hop_mod.o UDV_WRAP_mod.o \
Expand Down
10 changes: 10 additions & 0 deletions Prog/control_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ module Control
use files_mod
Use MyMats
use iso_fortran_env, only: output_unit, error_unit
Use Observables
Implicit none

real (Kind=Kind(0.d0)), private, save :: XMEANG, XMAXG, XMAXP, Xmean_tau, Xmax_tau
Expand All @@ -63,6 +64,7 @@ module Control

real (Kind=Kind(0.d0)), private, save :: XMAXP_HMC, XMEANP_HMC
Integer (Kind=Kind(0.d0)), private, save :: NC_Phase_HMC
Type (Obser_Vec ), public :: Obs_acc_HMC


real (Kind=Kind(0.d0)), private, save :: size_clust_Glob_up, size_clust_Glob_ACC_up
Expand All @@ -82,6 +84,8 @@ subroutine control_init(Group_Comm)
#endif
Implicit none
Integer , INTENT(IN) :: Group_Comm
Integer :: N
Character (len=64) :: Filename

XMEANG = 0.d0
XMEAN_tau = 0.d0
Expand All @@ -107,6 +111,9 @@ subroutine control_init(Group_Comm)
NC_Phase_HMC = 0
NC_HMC_up = 0
ACC_HMC_up = 0
N=1
Filename="Acc_HMC"
Call Obser_Vec_make(Obs_acc_HMC,N,Filename)

NC_Temp_up = 0
ACC_Temp_up = 0
Expand Down Expand Up @@ -198,8 +205,11 @@ Subroutine Control_upgrade_HMC(toggle)
Implicit none
Logical :: toggle
NC_HMC_up = NC_HMC_up + 1
Obs_acc_HMC%N = Obs_acc_HMC%N + 1
Obs_acc_HMC%Ave_sign = Obs_acc_HMC%Ave_sign + 1.d0
if (toggle) then
ACC_HMC_up = ACC_HMC_up + 1
Obs_acc_HMC%Obs_vec(1) = Obs_acc_HMC%Obs_vec(1) + 1.d0
endif
end Subroutine Control_upgrade_HMC

Expand Down
10 changes: 10 additions & 0 deletions Prog/main.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -764,6 +764,10 @@ Program Main
call system_clock(count_bin_start)

Call ham%Init_obs(Ltau)
Call Obser_vec_Init(Obs_acc_HMC)
do n=0, size(Langevin_HMC%Obs_ESJD_HMC,1)-1
Call Obser_vec_Init(Langevin_HMC%Obs_ESJD_HMC(n))
enddo
#if defined(TEMPERING)
Call Global_Tempering_init_obs
#endif
Expand Down Expand Up @@ -996,6 +1000,12 @@ Program Main

ENDDO
Call ham%Pr_obs(Ltau)
Call Print_bin_Vec(Obs_acc_HMC, Group_Comm)
do n=0, size(Langevin_HMC%Obs_ESJD_HMC,1)-1
if (Langevin_HMC%Obs_ESJD_HMC(n)%N > 0) then
Call Print_bin_Vec(Langevin_HMC%Obs_ESJD_HMC(n), Group_Comm)
endif
enddo
#if defined(TEMPERING)
Call Global_Tempering_Pr
#endif
Expand Down