Skip to content

Introduce Histogram observables #529

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 1 commit into
base: master
Choose a base branch
from
Open

Introduce Histogram observables #529

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
12 changes: 12 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 @@ -189,6 +189,7 @@ Module Hamiltonian_main
Type (Obser_Vec ), dimension(:), allocatable :: Obs_scal
Type (Obser_Latt), dimension(:), allocatable :: Obs_eq
Type (Obser_Latt), dimension(:), allocatable :: Obs_tau
Type (Obser_hist), dimension(:), allocatable :: Obs_hist


#include "Hamiltonians_interface.h"
Expand Down Expand Up @@ -481,6 +482,11 @@ Subroutine Pr_obs_base(LTAU)
Call Print_bin_Latt(Obs_tau(I), Group_Comm)
enddo
endif
if ( allocated(Obs_hist) ) then
Do I = 1,Size(Obs_hist,1)
Call Obs_hist(I)%print_bin(Group_Comm)
enddo
endif

end Subroutine Pr_obs_base

Expand Down Expand Up @@ -518,6 +524,12 @@ Subroutine Init_obs_base(Ltau)
Call Obser_Latt_Init(Obs_tau(I))
Enddo
Endif

if ( allocated(Obs_hist) ) then
Do I = 1,Size(Obs_hist,1)
Call Obs_hist(I)%init()
Enddo
endif

end Subroutine Init_obs_base

Expand Down
239 changes: 239 additions & 0 deletions Prog/observables_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,29 @@ Module Observables
procedure :: print_bin => print_bin_latt
end type Obser_Latt


type :: Obser_hist
!> Defines histogram observable
Integer :: N ! Number of measurements
real (Kind=Kind(0.d0)) :: Ave_Sign ! Averarge sign
Character (len=64) :: name ! Name of file in which the bins will be written out

!real (Kind=Kind(0.d0)) :: val ! temporary storage for value

integer :: N_classes ! Number of classes, in which the observable gets counted
real (Kind=Kind(0.d0)) :: upper, lower ! Range, in which the observable should lie.
real (Kind=Kind(0.d0)) :: d_class ! width of one class

integer, allocatable :: counts(:) !
integer :: N_below, N_above ! Counts occurences outside of range

contains
procedure :: make => Obser_hist_make
procedure :: init => Obser_hist_init
procedure :: print_bin => print_bin_hist
procedure :: measure => Obser_hist_measure
end type Obser_hist

Contains

Subroutine Obser_Latt_make(Obs, Nt, Filename, Latt, Latt_unit, Channel, dtau)
Expand Down Expand Up @@ -622,4 +645,220 @@ Subroutine Print_bin_Vec(Obs,Group_Comm)

End Subroutine Print_bin_Vec

!--------------------------------------------------------------------

Subroutine Obser_hist_make(obs, N_classes, upper, lower ,name)
Implicit none
class (Obser_hist) , intent(INOUT) :: obs
Integer , Intent(IN) :: N_classes
real (Kind=Kind(0.d0)), Intent(IN) :: upper, lower
Character (len=64) , Intent(IN) :: name

obs%N_classes = N_classes
obs%upper = upper
obs%lower = lower

obs%d_class = (upper - lower) / real(N_classes, Kind=Kind(0.d0))

Allocate ( obs%counts(N_classes) )
Obs%name = name
end subroutine Obser_hist_make
!--------------------------------------------------------------------

Subroutine Obser_hist_Init(obs)
Implicit none
class (Obser_hist), intent(INOUT) :: obs
obs%N = 0
obs%Ave_Sign= 0.d0

obs%counts(:) = 0
obs%N_below = 0
obs%N_above = 0
end subroutine Obser_hist_Init

!--------------------------------------------------------------------

Subroutine Obser_hist_measure(obs, value, sign_in)
Implicit none

class (Obser_hist), Intent(Inout) :: Obs
Real (Kind=Kind(0.d0)), Intent(In) :: value, sign_in

Integer :: n_x

obs%N = obs%N + 1
obs%Ave_sign = obs%Ave_sign + sign_in

n_x = ceiling( (value - obs%lower) / obs%d_class )
if ( n_x < 1 ) then
obs%N_below = obs%N_below + 1
elseif ( n_x > obs%N_classes ) then
obs%N_above = obs%N_above + 1
else
obs%counts(n_x) = obs%counts(n_x) + 1
endif

end Subroutine Obser_hist_measure

!--------------------------------------------------------------------

Subroutine Print_bin_hist(obs,Group_Comm)
#if defined MPI
Use mpi
#endif
#if defined HDF5
Use hdf5
Use alf_hdf5
#endif
Implicit none

class (Obser_hist), Intent(Inout) :: obs
Integer , Intent(In) :: Group_Comm

! Local
Integer :: I
Character (len=64) :: File_pr
real(Kind=Kind(0.d0)), pointer :: counts_out(:)
real(Kind=Kind(0.d0)), target :: N_above_out, N_below_out

#if defined HDF5
Character (len=7), parameter :: File_h5 = "data.h5"
Character (len=64) :: filename, groupname, obs_dsetname, sgn_dsetname, above_dsetname, below_dsetname
INTEGER(HID_T) :: file_id, group_id
logical :: link_exists
INTEGER :: hdferr
INTEGER(HSIZE_T), allocatable :: dims(:)
TYPE(C_PTR) :: dat_ptr
real(Kind=Kind(0.d0)), target :: sgn
#endif
#if defined MPI
Integer :: Ierr, Isize, Irank, No
INTEGER :: irank_g, isize_g, igroup
real (Kind=Kind(0.d0)), allocatable :: Tmp(:)
Real (Kind=Kind(0.d0)) :: X

CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR)
call MPI_Comm_rank(Group_Comm, irank_g, ierr)
call MPI_Comm_size(Group_Comm, isize_g, ierr)
igroup = irank/isize_g
#endif
allocate( counts_out(obs%N_classes) )
counts_out = dble(Obs%counts ) / dble(Obs%N)
N_above_out = dble(Obs%N_above) / dble(Obs%N)
N_below_out = dble(Obs%N_below) / dble(Obs%N)
Obs%Ave_sign = Obs%Ave_sign/dble(Obs%N)
write(File_pr, '(A,A)') trim(Obs%name), "_hist"
#if defined HDF5
groupname = File_pr
filename = File_h5
#endif

#if defined MPI
No = size(counts_out, 1)
Allocate (Tmp(No) )
Tmp = cmplx(0.d0,0.d0,kind(0.d0))
CALL MPI_REDUCE(counts_out,Tmp,No,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR)
counts_out = Tmp/DBLE(ISIZE_g)
deallocate (Tmp )

I = 1
X = 0.d0
CALL MPI_REDUCE(Obs%Ave_sign,X,I,MPI_REAL8,MPI_SUM, 0,Group_comm,IERR)
Obs%Ave_sign = X/DBLE(ISIZE_g)

X = 0.d0
CALL MPI_REDUCE(N_above_out,X,I,MPI_REAL8,MPI_SUM, 0,Group_comm,IERR)
N_above_out = X/DBLE(ISIZE_g)

X = 0.d0
CALL MPI_REDUCE(N_below_out,X,I,MPI_REAL8,MPI_SUM, 0,Group_comm,IERR)
N_below_out = X/DBLE(ISIZE_g)

if (Irank_g == 0 ) then
#endif

#if defined TEMPERING
write(File_pr,'(A,I0,A,A,A)') "Temp_",igroup,"/",trim(Obs%name), "_hist"
#if defined HDF5
write(filename ,'(A,I0,A,A)') "Temp_",igroup,"/",trim(File_h5)
#endif
#endif

#if defined HDF5
write(obs_dsetname,'(A,A,A)') trim(groupname), "/obser"
write(sgn_dsetname,'(A,A,A)') trim(groupname), "/sign"
write(above_dsetname,'(A,A,A)') trim(groupname), "/above"
write(below_dsetname,'(A,A,A)') trim(groupname), "/below"

CALL h5fopen_f (filename, H5F_ACC_RDWR_F, file_id, hdferr)

!Check if observable already exists in hdf5 file
CALL h5lexists_f(file_id, groupname, link_exists, hdferr)

if ( .not. link_exists ) then
!Create Group for observable
CALL h5gcreate_f (file_id, groupname, group_id, hdferr)
call write_attribute(group_id, '.', "N_classes", obs%N_classes, hdferr)
call write_attribute(group_id, '.', "upper" , obs%upper , hdferr)
call write_attribute(group_id, '.', "lower" , obs%lower , hdferr)
CALL h5gclose_f (group_id, hdferr)

!Create Dataset for data
allocate( dims(2) )
dims = (/size(counts_out,1), 0/)
CALL init_dset(file_id, obs_dsetname, dims, .false.)
deallocate( dims )

!Create Dataset for sign
allocate( dims(1) )
dims = (/0/)
CALL init_dset(file_id, sgn_dsetname, dims, .false.)
deallocate( dims )

!Create Dataset
allocate( dims(1) )
dims = (/0/)
CALL init_dset(file_id, above_dsetname, dims, .false.)
deallocate( dims )

!Create Dataset
allocate( dims(1) )
dims = (/0/)
CALL init_dset(file_id, below_dsetname, dims, .false.)
deallocate( dims )
endif

!Write data
dat_ptr = C_LOC(counts_out(1))
CALL append_dat(file_id, obs_dsetname, dat_ptr)

!Write sign
sgn = Obs%Ave_sign
dat_ptr = C_LOC(sgn)
CALL append_dat(file_id, sgn_dsetname, dat_ptr)

!Write
dat_ptr = C_LOC(N_above_out)
CALL append_dat(file_id, above_dsetname, dat_ptr)

!Write
dat_ptr = C_LOC(N_below_out)
CALL append_dat(file_id, below_dsetname, dat_ptr)

CALL h5fclose_f(file_id, hdferr)
#else
Open (Unit=10,File=File_pr, status="unknown", position="append")
WRITE(10,*) obs%N_classes, obs%upper, obs%lower, N_above_out, N_below_out, &
& (counts_out(I), I=1,size(counts_out,1)), Obs%Ave_sign
close(10)
#endif
#if defined MPI
endif
#endif
deallocate( counts_out )


end Subroutine Print_bin_hist

end Module Observables