diff --git a/Prog/Hamiltonian_main_mod.F90 b/Prog/Hamiltonian_main_mod.F90 index 2d4844cd0..997945345 100644 --- a/Prog/Hamiltonian_main_mod.F90 +++ b/Prog/Hamiltonian_main_mod.F90 @@ -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 @@ -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 diff --git a/Prog/Langevin_HMC_mod.F90 b/Prog/Langevin_HMC_mod.F90 index 47f461ab1..2e0bb23e0 100644 --- a/Prog/Langevin_HMC_mod.F90 +++ b/Prog/Langevin_HMC_mod.F90 @@ -36,6 +36,7 @@ Module Langevin_HMC_mod Use runtime_error_mod + use observables Use Hamiltonian_main Use UDV_State_mod Use Control @@ -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 @@ -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(:,:) @@ -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.) @@ -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) @@ -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) @@ -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 @@ -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 @@ -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 diff --git a/Prog/Makefile b/Prog/Makefile index 20e227ca2..e0925f576 100644 --- a/Prog/Makefile +++ b/Prog/Makefile @@ -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 \ diff --git a/Prog/control_mod.F90 b/Prog/control_mod.F90 index d99e2b535..3f2a6b788 100644 --- a/Prog/control_mod.F90 +++ b/Prog/control_mod.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/Prog/main.F90 b/Prog/main.F90 index 686fa3214..50c67d8a3 100644 --- a/Prog/main.F90 +++ b/Prog/main.F90 @@ -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 @@ -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