diff --git a/Prog/Hamiltonian_main_mod.F90 b/Prog/Hamiltonian_main_mod.F90 index 2d4844cd0..b496f68ce 100644 --- a/Prog/Hamiltonian_main_mod.F90 +++ b/Prog/Hamiltonian_main_mod.F90 @@ -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" @@ -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 @@ -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 diff --git a/Prog/observables_mod.F90 b/Prog/observables_mod.F90 index 8a2631488..9eed98f4b 100644 --- a/Prog/observables_mod.F90 +++ b/Prog/observables_mod.F90 @@ -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) @@ -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