diff --git a/Analysis/Max_SAC.F90 b/Analysis/Max_SAC.F90 index acec43762..0baf09b73 100644 --- a/Analysis/Max_SAC.F90 +++ b/Analysis/Max_SAC.F90 @@ -42,6 +42,7 @@ Program MaxEnt_Wrapper ! !-------------------------------------------------------------------- Use runtime_error_mod + Use Natural_Constants, only: pi, Eps_small Use MaxEnt_stoch_mod Use MaxEnt_mod use MaxEnt_Wrapper_mod @@ -56,10 +57,10 @@ Program MaxEnt_Wrapper Real (Kind=Kind(0.d0)), Dimension(:) , allocatable :: XQMC, XQMC_st, XTAU, Xtau_st, & & Alpha_tot, om_bf, alp_bf, xom, A Real (Kind=Kind(0.d0)), Dimension(:,:), allocatable :: XCOV, XCOV_st - Real (Kind=Kind(0.d0)) :: X_moments(2), Xerr_moments(2), ChiSq + Real (Kind=Kind(0.d0)) :: ChiSq Character (Len=64) :: command, File1, File2 Complex (Kind=Kind(0.d0)) :: Z - Logical :: Test =.false. + Logical, parameter :: Test =.false. Integer :: Ngamma, Ndis, NBins, NSweeps, Nwarm, N_alpha, N_cov @@ -73,9 +74,10 @@ Program MaxEnt_Wrapper Real (Kind=Kind(0.d0)), allocatable :: Xker_classic(:,:), A_classic(:), Default(:) Integer :: nt, nt1, io_error, n,nw, nwp, ntau, N_alpha_1, i, nbin_qmc - Integer :: ntau_st, ntau_en, ntau_new, Ntau_old - Real (Kind=Kind(0.d0)) :: dtau, pi, xmom1, x,x1,x2, tau, omp, om, Beta,err, delta, Dom - Real (Kind=Kind(0.d0)) :: Zero, Alpha_classic_st=100000.d0 + Integer :: ntau_st, ntau_en, Ntau_old + Real (Kind=Kind(0.d0)) :: dtau, xmom1, x,x1,x2, tau, omp, om, Beta,err, delta, Dom + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small + Real (Kind=Kind(0.d0)) :: Alpha_classic_st=100000.d0 Integer :: N_BZ_Zones = 1 Logical :: Extended_Zone = .false. @@ -159,8 +161,6 @@ Program MaxEnt_Wrapper Write(50,13) "Alpha_st",alpha_st Write(50,13) "R", R endif - Zero= 1.D-10 - pi = acos(-1.d0) Ntau_st = 1 Ntau_en = Ntau Select Case (str_to_upper(Channel)) @@ -420,7 +420,6 @@ Program MaxEnt_Wrapper ! Compute the real frequency Green function. delta = Dom - pi = acos(-1.d0) x = 0.d0 x1 = 0.d0 x2 = 0.d0 diff --git a/Analysis/ana_hdf5.F90 b/Analysis/ana_hdf5.F90 index 8191bd002..3babb08c0 100644 --- a/Analysis/ana_hdf5.F90 +++ b/Analysis/ana_hdf5.F90 @@ -48,7 +48,7 @@ Program ana_hdf5 #endif implicit none Integer :: i, hdferr, nargs, storage_type, nlinks, max_corder - Character (len=64) :: File_out, File_in, name + Character (len=64) :: File_in, name Character (len=64), allocatable :: names(:) INTEGER(HID_T) :: file_id, group_id INTEGER(HSIZE_T) :: n diff --git a/Analysis/ana_mod.F90 b/Analysis/ana_mod.F90 index 5aa98f1d7..38b7973f2 100644 --- a/Analysis/ana_mod.F90 +++ b/Analysis/ana_mod.F90 @@ -38,6 +38,7 @@ module ana_mod !> Collection of routines for postprocessing the Monte Carlo bins ! !-------------------------------------------------------------------- + Use Natural_Constants, only: Eps_small, Eps_convergence use iso_fortran_env, only: output_unit, error_unit use Files_mod Use Errors @@ -171,9 +172,8 @@ Subroutine read_vec_hdf5(filename, groupname, sgn, bins, analysis_mode) Character (len=64) :: file, obs_dsetname, sgn_dsetname INTEGER :: rank, hdferr - INTEGER(HSIZE_T) :: mem_dims(1) INTEGER(HSIZE_T), allocatable :: dims(:), maxdims(:) - INTEGER(HID_T) :: file_id, dset_id, dataspace, memspace + INTEGER(HID_T) :: file_id, dset_id, dataspace TYPE(C_PTR) :: dat_ptr file = 'data.h5' @@ -267,7 +267,7 @@ Subroutine read_local(file, sgn, bins, Latt, Latt_unit, dtau, Channel) Character (len=64) :: file_aux, str_temp1, str_temp2 Integer, allocatable :: List(:,:), Invlist(:,:) ! For orbital structure of Unit cell - Integer :: no, no1, n, nt, nb, Ntau, Ndim, Nbins, stat, Ndim_unit + Integer :: no, n, nt, nb, Ntau, Ndim, Nbins, stat, Ndim_unit Real(Kind=Kind(0.d0)) :: X Real(Kind=Kind(0.d0)), allocatable :: Xr_p(:,:), Orb_pos_temp(:) Real(Kind=Kind(0.d0)) :: x_p(2), a1_p(2), a2_p(2), L1_p(2), L2_p(2) @@ -1051,11 +1051,11 @@ Subroutine ana_tau(name_obs, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Ch Logical :: PartHole, L_Back, Extended_Zone Character (len=64) :: File_out, command, xk1_str, xk2_str - Real (Kind=Kind(0.d0)), parameter :: Zero=1.D-8 + Real (Kind=Kind(0.d0)), parameter :: Zero=Eps_small Integer :: N_skip, N_rebin, N_Cov, N_Back, N_auto, N_BZ_Zones Integer :: Nbins, LT, Lt_eff, n_mk Integer :: nb, no, no1, no2, n,i, nt, nt1, ierr, Norb, NBZ_1, NBZ_2 - Complex (Kind=Kind(0.d0)) :: Z, Zmean, Zerr + Complex (Kind=Kind(0.d0)) :: Z Real (Kind=Kind(0.d0)), allocatable :: Phase(:) Complex (Kind=Kind(0.d0)), allocatable :: V_help_loc(:,:,:,:), Bins_help(:,:,:,:) Real (Kind=Kind(0.d0)), allocatable :: Xk_p(:,:), Xk_p1(:), Xk_extended_p(:), X @@ -1156,8 +1156,7 @@ Subroutine ana_tau(name_obs, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Ch Xk_Extended_p(:) = Xk_p(:,n) + NBZ_1 * Latt%BZ1_p + NBZ_2 * Latt%BZ2_p if ( Xk_Extended_p(1) >= -zero .and. XK_Extended_p(2) >= -zero ) then L_back = .false. - if ( sqrt(Xk_extended_p(1)**2 + Xk_extended_p(2)**2) < 1.D-6 .and. N_Back == 1 ) L_back = .true. - ! Set weights + if ( sqrt(Xk_extended_p(1)**2 + Xk_extended_p(2)**2) < Eps_convergence .and. N_Back == 1 ) L_back = .true. If (Extended_Zone) then do no = 1,Norb X = 0.d0 @@ -1263,7 +1262,7 @@ Subroutine ana_tau(name_obs, sgn, bins_raw, bins0_raw, Latt, Latt_unit, dtau, Ch enddo if (PartHole) V_help_suscep = V_help_suscep * cmplx(2.d0,0.d0,Kind(0.d0)) L_back = .false. - if ( sqrt(Xk_p(1,n)**2 + Xk_p(2,n)**2) < 1.D-6 .and. N_Back == 1 ) L_back = .true. + if ( sqrt(Xk_p(1,n)**2 + Xk_p(2,n)**2) < Eps_convergence .and. N_Back == 1 ) L_back = .true. call COV(V_help_suscep, phase, Xcov, Xmean, background, L_back, N_rebin, Weights ) Xmean_st(n) = Xmean(1)*dtau Xerr_st(n) = Sqrt(Xcov(1,1))*dtau @@ -1333,12 +1332,11 @@ Subroutine ana_local(name, sgn, bins_raw, Latt, Latt_unit) Integer :: N_skip, N_rebin, N_Cov, N_Back, N_auto Integer :: Nbins, N_BZ_Zones Logical :: Extended_Zone - Integer :: i, n, nb, no, no1, ierr + Integer :: n, nb, no, ierr Complex (Kind=Kind(0.d0)), allocatable :: Phase(:) Complex (Kind=Kind(0.d0)), allocatable :: V_help(:) Real (Kind=Kind(0.d0)) :: Xr_p(2) - Complex (Kind=Kind(0.d0)) :: Z, Xmean, Xerr, Xmean_r, Xerr_r - Real (Kind=Kind(0.d0)) :: Xm,Xe + Complex (Kind=Kind(0.d0)) :: Xmean, Xerr NAMELIST /VAR_errors/ N_skip, N_rebin, N_Cov, N_Back, N_auto, N_BZ_Zones, Extended_Zone @@ -1542,7 +1540,7 @@ Subroutine ana_eq(name, sgn, bins_raw, bins0_raw, Latt, Latt_unit) do n = 1, Latt%N Xk_Extended_p(:) = Xk_p_s(:,n) + NBZ_1 * Latt%BZ1_p + NBZ_2 * Latt%BZ2_p L_back = .false. - if ( sqrt(Xk_Extended_p(1)**2 + Xk_Extended_p(2)**2) < 1.D-6 .and. N_Back == 1 ) L_back = .true. + if ( sqrt(Xk_Extended_p(1)**2 + Xk_Extended_p(2)**2) < Eps_convergence .and. N_Back == 1 ) L_back = .true. do no = 1,Latt_unit%Norb X = 0.d0 do i = 1,size(Latt%BZ1_p,1) @@ -1575,7 +1573,7 @@ Subroutine ana_eq(name, sgn, bins_raw, bins0_raw, Latt, Latt_unit) do nb = 1,Nbins V_help(1,nb) = bins (n,nb)%el(no,no1) enddo - if ( sqrt(Xk_p(1)**2 + Xk_p(2)**2) < 1.D-6 ) then + if ( sqrt(Xk_p(1)**2 + Xk_p(2)**2) < Eps_convergence ) then do nb = 1,Nbins V_help(2,nb) = Bins0(nb,no)*Latt%N V_help(3,nb) = Bins0(nb,no1) @@ -1865,7 +1863,7 @@ subroutine write_obs_vec_hdf5(filename, groupname, sgn, bins, analysis_mode) Integer :: Nobs, Nbins - INTEGER :: i, ierr + INTEGER :: ierr Character (len=64) :: obs_dsetname, sgn_dsetname INTEGER(HID_T) :: file_id, group_id logical :: file_exists, link_exists @@ -1981,7 +1979,7 @@ subroutine write_obs_latt_hdf5(filename, groupname, sgn, bins, bins0, Latt, Latt Integer :: Norb, Nunit, Ntau, Nbins - INTEGER :: i, ierr + INTEGER :: ierr Character (len=64) :: obs_dsetname, bak_dsetname, sgn_dsetname INTEGER(HID_T) :: file_id, group_id logical :: file_exists, link_exists @@ -2107,8 +2105,8 @@ subroutine write_obs_local_hdf5(filename, groupname, sgn, bins, Latt, Latt_unit, Integer :: Norb, Nunit, Ntau, Nbins - INTEGER :: i, ierr - Character (len=64) :: obs_dsetname, bak_dsetname, sgn_dsetname + INTEGER :: ierr + Character (len=64) :: obs_dsetname, sgn_dsetname INTEGER(HID_T) :: file_id, group_id logical :: file_exists, link_exists INTEGER(HSIZE_T), allocatable :: dims(:) diff --git a/Analysis/cov_eq.F90 b/Analysis/cov_eq.F90 index 6f139f370..b7eddd3a8 100644 --- a/Analysis/cov_eq.F90 +++ b/Analysis/cov_eq.F90 @@ -54,16 +54,15 @@ Program Cov_eq Integer :: Nunit, Norb, ierr - Integer :: no, no1, n, n1,m, nbins, n_skip, nb, N_rebin, N_cov, N_Back + Integer :: no, no1, n, m, nbins, n_skip, nb, N_rebin, N_cov, N_Back real (Kind=Kind(0.d0)):: X, Y Complex (Kind=Kind(0.d0)), allocatable :: Phase(:) Type (Mat_C), allocatable :: Bins (:,:), Bins_R(:,:) Complex (Kind=Kind(0.d0)), allocatable :: Bins0(:,:) Complex (Kind=Kind(0.d0)) :: Z, Xmean,Xerr, Xmean_r, Xerr_r Real (Kind=Kind(0.d0)) :: Xm,Xe - Real (Kind=Kind(0.d0)) :: Xk_p(2), XR_p(2) , XR1_p(2) - Complex (Kind=Kind(0.d0)), allocatable :: V_help(:), V_help_R(:) - Real (Kind=Kind(0.d0)) :: Pi, a1_p(2), a2_p(2), L1_p(2), L2_p(2), del_p(2) + Real (Kind=Kind(0.d0)) :: Xk_p(2), XR_p(2) + Complex (Kind=Kind(0.d0)), allocatable :: V_help(:) Real (Kind=Kind(0.d0)), allocatable :: AutoCorr(:),En(:) Integer :: L1, L2, I, N_auto, Ndim @@ -73,7 +72,6 @@ Program Cov_eq Type(Unit_cell) :: Latt_Unit Character (len=64) :: File_out - NAMELIST /VAR_Lattice/ L1, L2, Lattice_type, Model NAMELIST /VAR_errors/ n_skip, N_rebin, N_Cov, N_Back, N_auto @@ -101,7 +99,6 @@ Program Cov_eq ! Determine the number of bins. - Pi = acos(-1.d0) Open ( Unit=10, File="ineq", status="unknown" ) nbins = 0 do @@ -131,7 +128,7 @@ Program Cov_eq endif ! Allocate space - Allocate ( bins(Nunit,Nbins), bins_r(Nunit,Nbins), Phase(Nbins), V_help(Nbins), V_help_R(Nbins), Bins0(Nbins,Norb)) + Allocate ( bins(Nunit,Nbins), bins_r(Nunit,Nbins), Phase(Nbins), V_help(Nbins), Bins0(Nbins,Norb)) Do n = 1,Nunit do nb = 1,nbins Call Make_Mat(bins (n,nb),Norb) @@ -294,7 +291,6 @@ end Program Cov_eq !!$ Integer :: m !!$ !!$ Zero = 1.D-4 -!!$ pi = acos(-1.d0) !!$ X1_p(1) = Xk_p(2,n) !!$ X1_p(2) = -Xk_p(1,n) !!$ if (X1_p(1) < -pi + Zero ) X1_p(1) = X1_p(1) + 2.0*pi diff --git a/Analysis/cov_tau.F90 b/Analysis/cov_tau.F90 index a6a26ba5e..4cf6cba6d 100644 --- a/Analysis/cov_tau.F90 +++ b/Analysis/cov_tau.F90 @@ -42,6 +42,7 @@ Program Cov_tau Use Errors Use MyMats Use Matrix + Use Natural_Constants, only: Eps_small use iso_fortran_env, only: output_unit, error_unit #ifdef _OPENMP use check_omp_num_threads_mod @@ -53,11 +54,11 @@ Program Cov_tau Integer :: Nunit, Norb, N_auto Integer :: no, no1, n, nbins, n_skip, nb, NT, NT1, Lt, N_rebin, N_cov, ierr, N_Back Integer :: Lt_eff - real (Kind=Kind(0.d0)):: X, Y, dtau, X_diag + real (Kind=Kind(0.d0)):: X, Y, dtau Complex (Kind=Kind(0.d0)), allocatable :: Xmean(:), Xcov(:,:) Complex (Kind=Kind(0.d0)) :: Zmean, Zerr Complex (Kind=Kind(0.d0)) :: Z, Z_diag - Real (Kind=Kind(0.d0)) :: Zero=1.D-8 + Real (Kind=Kind(0.d0)), parameter :: Zero=Eps_small Real (Kind=Kind(0.d0)), allocatable :: Phase(:) Complex (Kind=Kind(0.d0)), allocatable :: PhaseI(:) Complex (Kind=Kind(0.d0)), allocatable :: Bins(:,:,:), Bins_chi(:,:), OneBin(:,:,:) @@ -309,11 +310,10 @@ end Program Cov_tau !!$ Real (Kind=Kind(0.d0)), INTENT(IN) :: Xk_p(2,Nunit) !!$ !!$ !Local -!!$ real (Kind=Kind(0.d0)) :: X1_p(2), Zero, pi, X +!!$ real (Kind=Kind(0.d0)) :: X1_p(2), X +!!$ real (Kind=Kind(0.d0)), parameter :: Zero = 1.D-4, pi = acos(-1.d0) !!$ Integer :: m !!$ -!!$ Zero = 1.D-4 -!!$ pi = acos(-1.d0) !!$ X1_p(1) = Xk_p(2,n) !!$ X1_p(2) = -Xk_p(1,n) !!$ if (X1_p(1) < -pi + Zero ) X1_p(1) = X1_p(1) + 2.0*pi diff --git a/Analysis/maxent_wrapper_mod.F90 b/Analysis/maxent_wrapper_mod.F90 index 178452c94..fa009d75f 100644 --- a/Analysis/maxent_wrapper_mod.F90 +++ b/Analysis/maxent_wrapper_mod.F90 @@ -32,6 +32,7 @@ module MaxEnt_Wrapper_mod Use MyMats + Use Natural_Constants, only: pi, Eps_small implicit none Real (Kind=Kind(0.d0)), allocatable, private :: Ra(:), ba(:) @@ -39,9 +40,7 @@ module MaxEnt_Wrapper_mod Real (Kind=Kind(0.d0)) function XKER_ph(tau,om, beta) Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta - - pi = 3.1415927 + real (Kind=Kind(0.d0)) :: tau, om, beta if ( beta*om > 600 ) then XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi @@ -56,9 +55,7 @@ end function XKER_ph Real (Kind=Kind(0.d0)) function XKER_ph_c(tau,om, beta) ! Kernal for A_c(om), same as XKER_ph Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta - - pi = 3.1415927 + real (Kind=Kind(0.d0)) :: tau, om, beta if ( beta*om > 600 ) then XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi @@ -73,9 +70,7 @@ end function XKER_ph_c Real (Kind=Kind(0.d0)) function XKER_pp(tau,om, beta) Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta - - pi = 3.1415927 + real (Kind=Kind(0.d0)) :: tau, om, beta if ( beta*om > 600 ) then XKER_pp = exp(-tau*om) / pi @@ -90,9 +85,7 @@ end function XKER_pp Real (Kind=Kind(0.d0)) function XKER_p(tau,om, beta) Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta - - pi = 3.1415927 + real (Kind=Kind(0.d0)) :: tau, om, beta if ( beta*om > 600 ) then XKER_p = exp(-tau*om) / pi @@ -107,9 +100,7 @@ end function XKER_p Real (Kind=Kind(0.d0)) function XKER_T0(tau,om, beta) Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta - - pi = 3.1415927 + real (Kind=Kind(0.d0)) :: tau, om, beta XKER_T0 = exp(-tau*om) / pi @@ -124,8 +115,6 @@ end function F Real (Kind=Kind(0.d0)) function F_QFI_ph(om, beta) Implicit None real (Kind=Kind(0.d0)) :: om, beta - real (Kind=Kind(0.d0)) :: pi - pi = 3.1415927 if ( abs(beta*om) > 600 ) then F_QFI_ph = 4.d0/pi else @@ -138,8 +127,6 @@ Real (Kind=Kind(0.d0)) function F_QFI_ph_c(om, beta) ! will improve Implicit None real (Kind=Kind(0.d0)) :: om, beta - real (Kind=Kind(0.d0)) :: pi - pi = 3.1415927 if ( abs(beta*om) > 600 ) then F_QFI_ph_c = 4.d0/pi else @@ -151,8 +138,6 @@ end function F_QFI_ph_c Real (Kind=Kind(0.d0)) function F_QFI_pp(om, beta) Implicit None real (Kind=Kind(0.d0)) :: om, beta - real (Kind=Kind(0.d0)) :: pi - pi = 3.1415927 if ( abs(beta*om) > 600 ) then F_QFI_pp = 4.d0/pi else @@ -169,14 +154,14 @@ Subroutine Set_Ra_ba(N) Integer, Intent(In) :: N Real (Kind=Kind(0.d0)), allocatable :: Mat(:,:), U(:,:), W(:) - Real (Kind=Kind(0.d0)) :: X, Y + Real (Kind=Kind(0.d0)) :: X, y Integer :: I, J , m , nc - Logical :: Test=.false. + Logical, parameter :: Test=.false. allocate(Mat(N,N), U(N,N), W(N)) allocate(Ra(N/2),ba(N/2)) - If (Test) Write(6,*) "Setting Ra and ba using the method of Karrasch of Phys. Rev. B 82 (2010), 125114" + if (Test) Write(6,*) "Setting Ra and ba using Karrasch, Phys. Rev. B 82 (2010), 125114" Mat = 0.d0 do i = 1,N-1 @@ -200,7 +185,7 @@ Subroutine Set_Ra_ba(N) X = X + Mat(m,j)*U(j,i) enddo X = X -W(i)*U(m,i) - if (Abs(X) >= 1.d-10) then + if (Abs(X) >= Eps_small) then Write(6,*) ABS(X) write(error_unit,*) "Issue with eigenvalue in subroutine Set_Ra_ba of mod maxent_wrapper_mod" CALL Terminate_on_error(ERROR_MAXENT,__FILE__,__LINE__) @@ -208,7 +193,7 @@ Subroutine Set_Ra_ba(N) enddo enddo - If (Test) then + if (Test) then Open(Unit=10,File="Ra_ba.dat", status="Unknown") Do i = 1,size(ba,1) write(10,*) Ra(i),ba(i) @@ -315,7 +300,7 @@ Real (Kind=Kind(0.d0)) function Back_trans_ph_c(Aom, om, beta) real (Kind=Kind(0.d0)) :: Zero ! same as Back_trans_pp, since Back_trans_pp gives = chi(q,om)/omega ! = A(q,om)*tanh(beta om/2)/om - Zero = 1.D-8 + Zero = Eps_small if ( abs(om) < zero ) then Back_trans_ph_c = beta * Aom/2.d0 elseif ( beta*om > 600 ) then @@ -339,7 +324,7 @@ Real (Kind=Kind(0.d0)) function Back_trans_pp(Aom, om, beta) real (Kind=Kind(0.d0)), intent(in) :: Aom, beta, om real (Kind=Kind(0.d0)) :: Zero - Zero = 1.D-8 + Zero = Eps_small if ( abs(om) < Zero ) then Back_trans_pp = beta * Aom/2.d0 elseif ( beta*om > 600 ) then @@ -360,9 +345,8 @@ end function BACK_TRANS_PP Real (Kind=Kind(0.d0)) function XKER_p_ph(tau,om, beta) Implicit None - real (Kind=Kind(0.d0)) :: tau, om, pi, beta + real (Kind=Kind(0.d0)) :: tau, om, beta - pi = 3.1415927 if ( beta*om > 600 ) then XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi @@ -488,7 +472,8 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e Character (Len=*), INTENT(IN) :: Channel Logical, INTENT(IN) :: Default_model_exists, Stochastic Integer :: Ndis, Nw - Real (Kind = Kind(0.d0)) :: Dom, X, Om, Zero = 1.D-8 + Real (Kind = Kind(0.d0)) :: Dom, X, Om + Real (Kind = Kind(0.d0)), parameter :: Zero = Eps_small Ndis = size(Default,1) Dom = (OM_en - Om_st)/dble(Ndis) @@ -528,7 +513,7 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e ! See Back_trans_ph_c/Back_trans_pp ! Default(om) = (1 - exp(-beta*om))/(1 + exp(-beta*om))*A(om)/om Om = Om_st + dble(nw)*dom - if ( abs(om) < zero ) then + if ( abs(om) < Zero ) then Default(nw) = Default(nw)*2.d0/ beta elseif ( beta*om > 600 ) then Default(nw) = Default(nw) * om @@ -557,7 +542,7 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e X = 0.d0 Do nw = 1, Ndis Om = Om_st + dble(nw)*dom - if ( abs(om) < zero ) then + if ( abs(om) < Zero ) then Default(nw) = Default(nw)*2.d0/ beta elseif ( beta*om > 600 ) then Default(nw) = Default(nw) * om diff --git a/Libraries/Modules/Makefile b/Libraries/Modules/Makefile index 6f094b9cf..708acc4cf 100644 --- a/Libraries/Modules/Makefile +++ b/Libraries/Modules/Makefile @@ -67,6 +67,9 @@ libfakeintel.so: fakeintel.c echo "No C compiler found, libfakeintel.so not created." 1>&2; \ fi +entanglement_mod.o: entanglement_mod.F90 + $(ALF_FC) -c $(shell echo $(ALF_FLAGS_MODULES) | sed 's| -pedantic||') $< + $(DEPSFILE): $(SRCS) $(DEPSGEN) python3 $(DEPSGEN) $(SRCS) > $(DEPSFILE) diff --git a/Libraries/Modules/alf_hdf5_mod.F90 b/Libraries/Modules/alf_hdf5_mod.F90 index 7ff13eb88..0a00ee147 100644 --- a/Libraries/Modules/alf_hdf5_mod.F90 +++ b/Libraries/Modules/alf_hdf5_mod.F90 @@ -47,6 +47,7 @@ Module alf_hdf5 Use hdf5 use h5lt + Use Natural_Constants, only: Eps_small Use Lattices_v3 implicit none @@ -738,7 +739,7 @@ Subroutine test_attribute_double(loc_id, obj_name, attr_name, attr_value, ierr) INTEGER, INTENT(OUT) :: ierr LOGICAL :: attr_exists - real(Kind=Kind(0.d0)), parameter :: ZERO = 10D-8 + real(Kind=Kind(0.d0)), parameter :: ZERO = Eps_small real(Kind=Kind(0.d0)) :: test_double, diff call h5aexists_by_name_f(loc_id, obj_name, attr_name, attr_exists, ierr) @@ -1002,5 +1003,9 @@ Subroutine write_comment(loc_id, obj_name, attr_name, comment, ierr) end Subroutine write_comment - end Module alf_hdf5 + end Module alf_hdf5 +#else + module alf_hdf5_disabled_mod + implicit none + end module alf_hdf5_disabled_mod #endif diff --git a/Libraries/Modules/entanglement_mod.F90 b/Libraries/Modules/entanglement_mod.F90 index a44edad67..102167d60 100644 --- a/Libraries/Modules/entanglement_mod.F90 +++ b/Libraries/Modules/entanglement_mod.F90 @@ -64,9 +64,8 @@ Subroutine Init_Entanglement_replicas(Group_Comm) #endif Implicit none Integer, INTENT(IN) :: Group_Comm - - Integer :: ISIZE, IRANK, IERR #ifdef MPI + Integer :: ISIZE, IRANK, IERR ! Create subgroups of two replicas each CALL MPI_COMM_SIZE(Group_Comm,ISIZE,IERR) CALL MPI_COMM_RANK(Group_Comm,IRANK,IERR) @@ -96,7 +95,7 @@ Subroutine Calc_Mutual_Inf_indep(GRC,List_A,Nsites_A,List_B,Nsites_B,N_SUN,Renyi Complex (kind=kind(0.d0)), INTENT(OUT) :: Renyi_A, Renyi_B, Renyi_AB Integer, Dimension(:), Allocatable :: List_AB - Integer :: I, J, IERR, INFO, Nsites_AB + Integer :: I, Nsites_AB Nsites_AB=Nsites_A+Nsites_B @@ -132,7 +131,7 @@ Subroutine Calc_Mutual_Inf_gen_fl(GRC,List_A,Nsites_A,List_B,Nsites_B,N_SUN,Reny Complex (kind=kind(0.d0)), INTENT(OUT) :: Renyi_A, Renyi_B, Renyi_AB Integer, Allocatable :: List_AB(:,:), Nsites_AB(:) - Integer :: I, J, IERR, INFO, N_FL, Nsites_AB_max + Integer :: I, J, N_FL, Nsites_AB_max N_FL=size(GRC,3) Allocate(Nsites_AB(N_FL)) @@ -176,7 +175,7 @@ Subroutine Calc_Mutual_Inf_gen_all(GRC,List_A,Nsites_A,List_B,Nsites_B,Renyi_A,R Complex (kind=kind(0.d0)), INTENT(OUT) :: Renyi_A, Renyi_B, Renyi_AB Integer, Allocatable :: List_AB(:,:,:), Nsites_AB(:,:) - Integer :: I, J, IERR, INFO, N_FL, Nsites_AB_max, nc, num_nc + Integer :: I, J, N_FL, Nsites_AB_max, nc, num_nc N_FL=size(GRC,3) num_nc=size(List_B,3) @@ -248,17 +247,17 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_indep(GRC,List,Nsites,N_SUN) Integer, INTENT(IN) :: List(:) Integer, INTENT(IN) :: Nsites, N_SUN +#ifdef MPI Complex (kind=kind(0.d0)), Dimension(:,:), Allocatable :: GreenA, GreenA_tmp, IDA ! Integer, Dimension(:), Allocatable :: PIVOT - Complex (kind=kind(0.d0)) :: DET, PRODDET, alpha, beta - Integer :: J, IERR, INFO, N_FL, nf, N_FL_half + Complex (kind=kind(0.d0)) :: PRODDET, alpha, beta + Integer :: J, N_FL, nf, N_FL_half Integer , Dimension(:,:), Allocatable :: List_tmp Integer , Dimension(2) :: Nsites_tmp,nf_list,N_SUN_tmp - Logical, save :: First_call = .True. EXTERNAL ZGEMM EXTERNAL ZGETRF - + N_FL = size(GRC,3) N_FL_half = N_FL/2 @@ -271,7 +270,6 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_indep(GRC,List,Nsites,N_SUN) return endif -#ifdef MPI ! Check if entanglement replica group is of size 2 such that the second renyi entropy can be calculated if(ENT_SIZE==2) then @@ -312,11 +310,14 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_indep(GRC,List,Nsites,N_SUN) ! At this point, each task of the temepering group / world returns the same averaged value of the pairs, including the possible "free"/ unpaired one. ! This mechanisms leads to some syncronization, but I (Johannes) am lacking a better way to treat odd number of tasks. #else + Block + Logical, save :: First_call = .True. Calc_Renyi_Ent_indep=0.0d0 if (First_call) then write(error_unit,*) "Entanglement module compiled without MPI, no Renyi entropy results possible!" First_call = .false. endif + End Block #endif End function Calc_Renyi_Ent_indep @@ -361,15 +362,14 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_fl(GRC,List,Nsites,N_SUN) Integer, INTENT(IN) :: List(:,:) Integer, INTENT(IN) :: Nsites(:), N_SUN(:) ! new +#ifdef MPI Complex (kind=kind(0.d0)), Dimension(:,:), Allocatable :: GreenA, GreenA_tmp, IDA ! Integer, Dimension(:), Allocatable :: PIVOT - Complex (kind=kind(0.d0)) :: DET, PRODDET, alpha, beta - Integer :: I, J, IERR, INFO, N_FL, nf, N_FL_half, x, dim, dim_eff, nf_eff, start_flav + Complex (kind=kind(0.d0)) :: PRODDET, alpha, beta + Integer :: I, J, N_FL, nf, N_FL_half, x, dim, nf_eff, start_flav Integer , Dimension(:), Allocatable :: SortedFlavors ! new Integer , Dimension(:,:), Allocatable :: List_tmp Integer , Dimension(2) :: Nsites_tmp,nf_list,N_SUN_tmp - Logical, save :: First_call = .True. - EXTERNAL ZGEMM EXTERNAL ZGETRF @@ -401,8 +401,7 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_fl(GRC,List,Nsites,N_SUN) Calc_Renyi_Ent_gen_fl=CMPLX(1.d0,0.d0,kind(0.d0)) alpha=CMPLX(2.d0,0.d0,kind(0.d0)) beta =CMPLX(1.d0,0.d0,kind(0.d0)) - -#ifdef MPI + ! Check if entanglement replica group is of size 2 such that the second reny entropy can be calculated if(ENT_SIZE==2) then @@ -447,11 +446,14 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_fl(GRC,List,Nsites,N_SUN) ! At this point, each task of the temepering group / world returns the same averaged value of the pairs, including the possible "free"/ unpaired one. ! This mechanisms leads to some syncronization, but I (Johannes) am lacking a better way to treat odd number of tasks. #else + Block + Logical, save :: First_call = .True. Calc_Renyi_Ent_gen_fl=0.0d0 if (First_call) then write(error_unit,*) "Entanglement module compiled without MPI, no Renyi entropy results possible!" First_call = .false. endif + End Block #endif End function Calc_Renyi_Ent_gen_fl @@ -493,16 +495,15 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_all(GRC,List,Nsites) Integer, Dimension(:,:,:), INTENT(IN) :: List Integer, INTENT(IN) :: Nsites(:,:) +#ifdef MPI Complex (kind=kind(0.d0)), Dimension(:,:), Allocatable :: GreenA, GreenA_tmp, IDA ! Integer, Dimension(:), Allocatable :: PIVOT - Complex (kind=kind(0.d0)) :: DET, PRODDET, alpha, beta - Integer :: I, J, IERR, INFO, N_FL, nf, N_FL_half, x, dim, dim_eff, nf_eff, start_flav + Complex (kind=kind(0.d0)) :: PRODDET, alpha, beta + Integer :: I, J, N_FL, nf, N_FL_half, x, dim, start_flav Integer :: nc, num_nc - Integer , Dimension(:), Allocatable :: SortedFlavors,N_SUN_fl,df_list + Integer , Dimension(:), Allocatable :: SortedFlavors,N_SUN_fl Integer , Dimension(:,:), Allocatable :: List_tmp, eff_ind, eff_ind_inv Integer , Dimension(2) :: Nsites_tmp,nf_list,N_SUN_tmp - Logical, save :: First_call = .True. - EXTERNAL ZGEMM EXTERNAL ZGETRF @@ -548,8 +549,7 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_all(GRC,List,Nsites) Calc_Renyi_Ent_gen_all=CMPLX(1.d0,0.d0,kind(0.d0)) alpha=CMPLX(2.d0,0.d0,kind(0.d0)) beta =CMPLX(1.d0,0.d0,kind(0.d0)) - -#ifdef MPI + ! Check if entanglement replica group is of size 2 such that the second reny entropy can be calculated if(ENT_SIZE==2) then @@ -599,11 +599,14 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_gen_all(GRC,List,Nsites) ! At this point, each task of the temepering group / world returns the same averaged value of the pairs, including the possible "free"/ unpaired one. ! This mechanisms leads to some syncronization, but I (Johannes) am lacking a better way to treat odd number of tasks. #else + Block + Logical, save :: First_call = .True. Calc_Renyi_Ent_gen_all=0.0d0 if (First_call) then write(error_unit,*) "Entanglement module compiled without MPI, no Renyi entropy results possible!" First_call = .false. endif + End Block #endif @@ -673,8 +676,7 @@ Complex (kind=kind(0.d0)) function Calc_Renyi_Ent_pair(GRC,List,Nsites,nf_list,N Integer, Dimension(:), Allocatable :: PIVOT Complex (kind=kind(0.d0)) :: DET, PRODDET, alpha, beta - Integer :: I, J, IERR, INFO, N_FL, nf, N_FL_half, x, dim, dim_eff, nf_eff, start_flav, dim_sq - Integer , Dimension(:), Allocatable :: SortedFlavors ! new + Integer :: I, J, IERR, INFO, dim, dim_eff, nf_eff, dim_sq Calc_Renyi_Ent_pair=CMPLX(1.d0,0.d0,kind(0.d0)) @@ -796,8 +798,7 @@ Complex (Kind=8) function Calc_Renyi_Ent_single(GRC,List,Nsites,nf_eff,N_SUN,Gre Integer, Dimension(:), Allocatable :: PIVOT Complex (kind=kind(0.d0)) :: DET, PRODDET, alpha, beta - Integer :: I, J, IERR, INFO, N_FL, nf, N_FL_half, x, dim, dim_eff, start_flav, dim_sq - Integer , Dimension(:), Allocatable :: SortedFlavors ! new + Integer :: I, J, IERR, INFO, dim, dim_eff, dim_sq Calc_Renyi_Ent_single=CMPLX(1.d0,0.d0,kind(0.d0)) alpha=CMPLX(2.d0,0.d0,kind(0.d0)) diff --git a/Libraries/Modules/fourier_mod.F90 b/Libraries/Modules/fourier_mod.F90 index 922c05eea..086fcad0c 100644 --- a/Libraries/Modules/fourier_mod.F90 +++ b/Libraries/Modules/fourier_mod.F90 @@ -1,6 +1,7 @@ Module Fourier use runtime_error_mod + Use Natural_Constants, only: Eps_small Use MaxEnt_mod Use MaxEnt_stoch_mod Use Matrix @@ -650,7 +651,7 @@ subroutine Tau_Matz_T_Bose(griom, xiom, grtau, xtau, beta, A, xom, cov) Real (Kind=Kind(0.d0)) :: Alpha_st, Chisq, x, Zero Nom = Size(Xom ,1) - Zero = 1.D-10 + Zero = Eps_small Do Nw = 1,Nom if ( xom(Nw) .lt. -Zero ) then Write(error_unit,*) 'Tau_Matz_T_Bose: Frequencies should be larger than zero' diff --git a/Libraries/Modules/lattices_v3_mod.F90 b/Libraries/Modules/lattices_v3_mod.F90 index 083e81f9f..6e8c76a0d 100644 --- a/Libraries/Modules/lattices_v3_mod.F90 +++ b/Libraries/Modules/lattices_v3_mod.F90 @@ -42,6 +42,7 @@ Module Lattices_v3 !-------------------------------------------------------------------- Use Matrix Use runtime_error_mod + Use Natural_Constants, only: twopi, Eps_small use iso_fortran_env, only: output_unit, error_unit implicit none private @@ -100,15 +101,16 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) Real (Kind=Kind(0.d0)), dimension(:), allocatable :: x_p, x1_p, a_p,d_p Real (Kind=Kind(0.d0)), allocatable :: Mat(:,:), Mat_inv(:,:) - Integer :: ndim, L, L1, nc, i, i1,i2, L_f, LQ, n,m, nd1,nd2,nr, nnr1, nnr2, nnr, nr1, imj_1, imj_2 + Integer :: L, L1, nc, i, i1,i2, L_f, LQ, n,m, nd1,nd2,nr, nnr1, nnr2, nnr, nr1, imj_1, imj_2 Integer :: imj - Real (Kind=Kind(0.d0)) :: Zero,pi, X + Real (Kind=Kind(0.d0)) :: X + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small - ndim = size(L1_p) + Integer, parameter :: ndim = 2 + ! ndim = size(L1_p) TODO: Generalize to 3D allocate (Latt%L2_p(ndim), Latt%L1_p(ndim), Latt%a1_p(ndim) , Latt%a2_p(ndim), & & Latt%b1_p(ndim), Latt%b2_p(ndim), Latt%BZ1_p(ndim), Latt%BZ2_p(ndim) ) allocate (Latt%b1_perp_p(ndim), Latt%b2_perp_p(ndim) ) - Zero = 1.D-5 Latt%L1_p = L1_p Latt%L2_p = L2_p Latt%a1_p = a1_p @@ -121,8 +123,6 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) Allocate ( x_p(ndim), x1_p(ndim), d_p(ndim), a_p(ndim) ) - pi = acos(-1.d0) - ! Setup the 2X2 matrix to determine BZ1_p, BZ2_p Allocate ( Mat(2 , 2), Mat_inv( 2 , 2 ) ) Mat(1,1) = dble(a1_p(1)) @@ -134,10 +134,10 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) Mat_inv(2,2) = Mat(1,1)/X Mat_inv(1,2) = -Mat(1,2)/X Mat_inv(2,1) = -Mat(2,1)/X - BZ1_p(1) = 2.d0*pi*Mat_inv(1,1) - BZ1_p(2) = 2.d0*pi*Mat_inv(2,1) - BZ2_p(1) = 2.d0*pi*Mat_inv(1,2) - BZ2_p(2) = 2.d0*pi*Mat_inv(2,2) + BZ1_p(1) = twopi*Mat_inv(1,1) + BZ1_p(2) = twopi*Mat_inv(2,1) + BZ2_p(1) = twopi*Mat_inv(1,2) + BZ2_p(2) = twopi*Mat_inv(2,2) Latt%BZ1_p = BZ1_p Latt%BZ2_p = BZ2_p @@ -145,7 +145,7 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) ! K-space Quantization from periodicity in L1_p and L2_p - X = 2.d0*pi / ( Iscalar(BZ1_p,L1_p) * Iscalar(BZ2_p,L2_p) - & + X = twopi / ( Iscalar(BZ1_p,L1_p) * Iscalar(BZ2_p,L2_p) - & & Iscalar(BZ2_p,L1_p) * Iscalar(BZ1_p,L2_p) ) X = abs(X) b1_p = X*( Iscalar(BZ2_p,L2_p) * BZ1_p - Iscalar(BZ1_p,L2_p) * BZ2_p ) @@ -174,12 +174,12 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) ! Count the number of lattice points. - L = abs(nint ( Iscalar(Latt%BZ1_p,L1_p) / (2.d0*pi) )) - L1 = abs(nint ( Iscalar(Latt%BZ2_p,L1_p) / (2.d0*pi) )) + L = abs(nint ( Iscalar(Latt%BZ1_p,L1_p) / twopi )) + L1 = abs(nint ( Iscalar(Latt%BZ2_p,L1_p) / twopi )) if (L1 .gt. L) L = L1 - L1 = abs(nint ( Iscalar(Latt%BZ1_p,L2_p) / (2.d0*pi) )) + L1 = abs(nint ( Iscalar(Latt%BZ1_p,L2_p) / twopi )) if (L1 .gt. L) L = L1 - L1 = abs(nint ( Iscalar(Latt%BZ2_p,L2_p) / (2.d0*pi) )) + L1 = abs(nint ( Iscalar(Latt%BZ2_p,L2_p) / twopi )) if (L1 .gt. L) L = L1 nc = 0 do i1 = -L,L @@ -286,8 +286,8 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) call npbc(x_p , x1_p, Latt%L1_p, Latt%L2_p) call npbc(x1_p, x_p , Latt%L1_p, Latt%L2_p) call npbc(x_p , x1_p, Latt%L1_p, Latt%L2_p) - nnr1 = nint ( Iscalar(Latt%BZ1_p,x_p) / (2.d0*pi) ) - nnr2 = nint ( Iscalar(Latt%BZ2_p,x_p) / (2.d0*pi) ) + nnr1 = nint ( Iscalar(Latt%BZ1_p,x_p) / twopi ) + nnr2 = nint ( Iscalar(Latt%BZ2_p,x_p) / twopi ) nnr = Latt%invlist(nnr1,nnr2) Latt%nnlist(nr,nd1,nd2) = nnr if ( nnr < 1 .or. nnr > Latt%N ) then @@ -322,8 +322,8 @@ subroutine Make_lattice(L1_p, L2_p, a1_p, a2_p, Latt) d_p = x_p - x1_p call npbc(x1_p , d_p , Latt%L1_p, Latt%L2_p) call npbc(d_p , x1_p, Latt%L1_p, Latt%L2_p) - imj_1 = nint ( Iscalar(Latt%BZ1_p,d_p) / (2.d0*pi) ) - imj_2 = nint ( Iscalar(Latt%BZ2_p,d_p) / (2.d0*pi) ) + imj_1 = nint ( Iscalar(Latt%BZ1_p,d_p) / twopi ) + imj_2 = nint ( Iscalar(Latt%BZ2_p,d_p) / twopi ) imj = Latt%invlist(imj_1,imj_2) Latt%imj(nr,nr1) = imj enddo @@ -356,10 +356,10 @@ subroutine npbc_I(nr_p, n_p, L1_p, L2_p) integer, dimension(:), intent(out) :: nr_p integer, dimension(:), allocatable :: x_p - Real (Kind=Kind(0.d0)) :: Zero, X + Real (Kind=Kind(0.d0)) :: X + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small Integer :: Ndim, i - Zero = 1.D-8 nr_p = n_p ndim = size(n_p) @@ -389,10 +389,10 @@ subroutine npbc_R(nr_p, n_p, L1_p, L2_p) Real (Kind=Kind(0.d0)), dimension(:), allocatable :: x_p - Real (Kind=Kind(0.d0)) :: Zero, X + Real (Kind=Kind(0.d0)) :: X + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small Integer :: ndim, i - Zero = 1.D-8 nr_p = n_p ndim = size(n_p) allocate(x_p(ndim)) @@ -422,10 +422,10 @@ subroutine npbc_R_B(nr_p, n_p, L1_p, L2_p, N1, N2 ) Real (Kind=Kind(0.d0)), dimension(:), allocatable :: x_p - Real (Kind=Kind(0.d0)) :: Zero, X + Real (Kind=Kind(0.d0)) :: X + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small Integer :: ndim, i, Del_N1, Del_N2 - Zero = 1.D-8 nr_p = n_p ndim = size(n_p) allocate(x_p(ndim)) @@ -467,7 +467,8 @@ integer Function Inv_K(XK_P,Latt) Type (Lattice) :: Latt Integer :: nkx, nky, nk - Real (Kind=Kind(0.d0)) :: XK1_P(2), XK2_P(2), Zero + Real (Kind=Kind(0.d0)) :: XK1_P(2), XK2_P(2) + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small call npbc(xk1_p, xk_p , Latt%BZ1_p, Latt%BZ2_p) call npbc(xk2_p, xk1_p, Latt%BZ1_p, Latt%BZ2_p) @@ -477,7 +478,6 @@ integer Function Inv_K(XK_P,Latt) nk = Latt%Invlistk(nkx,nky) !Test - Zero = 1.D-10 XK1_P = Latt%listk(nk,1)*latt%b1_p + Latt%listk(nk,2)*latt%b2_p XK1_P = XK1_P - XK2_P if (Xnorm(XK1_P) < Zero ) then @@ -515,14 +515,11 @@ integer Function Inv_R(XR_P,Latt) Real (Kind=Kind(0.d0)) :: XR1_P(2), XR2_P(2) Integer :: n_1, n_2 - Real (Kind=Kind(0.d0)) :: pi - - pi = acos(-1.d0) call npbc(xr1_p, xr_p , Latt%L1_p, Latt%L2_p) call npbc(xr2_p, xr1_p, Latt%L1_p, Latt%L2_p) - n_1 = nint ( Iscalar(Latt%BZ1_p,XR2_p) / (2.d0*pi) ) - n_2 = nint ( Iscalar(Latt%BZ2_p,XR2_p) / (2.d0*pi) ) + n_1 = nint ( Iscalar(Latt%BZ1_p,XR2_p) / twopi ) + n_2 = nint ( Iscalar(Latt%BZ2_p,XR2_p) / twopi ) Inv_R = Latt%invlist(n_1,n_2) end Function Inv_R diff --git a/Libraries/Modules/maxent_mod.F90 b/Libraries/Modules/maxent_mod.F90 index 69d50b839..ebce0109f 100644 --- a/Libraries/Modules/maxent_mod.F90 +++ b/Libraries/Modules/maxent_mod.F90 @@ -47,6 +47,7 @@ Module MaxEnt_mod ! !-------------------------------------------------------------------- + Use Natural_Constants, only: Eps_small, Eps_convergence Use MyMats Use Errors use iso_fortran_env, only: output_unit, error_unit @@ -82,7 +83,7 @@ Subroutine MaxEnt_T( XQMC, COV, A, XKER, ALPHA_ST, CHISQ,DEFAULT) !WRITE(6,*) 'NTAU, Nom: ', NTAU,NOM Xmom1 = Xqmc(1) - ZERO = 1.0D-8 + ZERO = Eps_small ALLOCATE ( XLAM(NTAU), SIG1(NTAU), COVM1(NTAU,NTAU), UC(NTAU,NTAU), DEF(NOM) ) XLAM=0.D0; SIG1=0.D0; UC = 0.D0 @@ -268,7 +269,7 @@ Subroutine Maximize_Newton( XQMC, COV, A, XKER, XQ,XENT,CHISQ) ENDDO NITER = NITER + 1 !WRITE(6,*) 'Maximize: ', XNORM, NITER - IF (XNORM.LT.1.0D-6 .OR. NITER.GE.500) EXIT + IF (XNORM.LT.Eps_convergence .OR. NITER.GE.500) EXIT ENDDO CALL SETQ(A,XKER,XQMC, XQ,XENT,CHISQ) @@ -336,7 +337,7 @@ Subroutine Maximize_Self( XQMC, COV, A, XKER, XQ,XENT,CHISQ) ENDDO NITER = NITER + 1 !WRITE(6,*) 'Maximize_Self: ', XNORM, NITER - IF (XNORM.LT.1.0D-6 .OR. NITER.GE.1000) EXIT + IF (XNORM.LT.Eps_convergence .OR. NITER.GE.1000) EXIT ENDDO CALL SETQ(A,XKER,XQMC, XQ,XENT,CHISQ) @@ -685,7 +686,8 @@ Subroutine MaxEnt_T_Bryan( XQMC, COV, A, XKER, ALPHA_ST, ALPHA_EN, CHISQ ) Implicit None Real (Kind=Kind(0.d0)), Dimension(:) :: XQMC, A Real (Kind=Kind(0.d0)), Dimension(:,:) :: COV, XKER - Real (Kind=Kind(0.d0)) :: ALPHA_ST, ALPHA_N, ALPHA_EN ! , PI + Real (Kind=Kind(0.d0)) :: ALPHA_ST, ALPHA_N, ALPHA_EN + ! Real (Kind=Kind(0.d0)), parameter :: PI = acos(-1.d0) Real (Kind=Kind(0.D0)), Intent(Out) :: CHISQ Integer :: NT, NT1, NT2, NW, NCOUNT, NTH @@ -700,9 +702,8 @@ Subroutine MaxEnt_T_Bryan( XQMC, COV, A, XKER, ALPHA_ST, ALPHA_EN, CHISQ ) NOM = SIZE(A, 1) ALLOCATE(A_ME(NOM)) !WRITE(6,*) 'NTAU, Nom: ', NTAU,NOM -! PI = ACOS(-1.d0) XMOM1= 1.0D0 !PI - ZERO = 1.0D-8 + ZERO = Eps_small ALLOCATE ( XLAM(NTAU), SIG1(NTAU), COVM1(NTAU,NTAU), UC(NTAU,NTAU), DEF(NOM) ) XLAM=0.D0; SIG1=0.D0; UC = 0.D0 diff --git a/Libraries/Modules/maxent_stoch_mod.F90 b/Libraries/Modules/maxent_stoch_mod.F90 index 2fad95ad6..232d86e04 100644 --- a/Libraries/Modules/maxent_stoch_mod.F90 +++ b/Libraries/Modules/maxent_stoch_mod.F90 @@ -50,7 +50,7 @@ Module MaxEnt_stoch_mod Integer, private :: NTAU, nt, Ngamma, ng, Ndis, nd, L_seed Integer, private, allocatable:: Iseed_vec(:) - Real (Kind=Kind(0.d0)), private :: Delta, Delta2, OM_st_1, Om_en_1, DeltaXMAX, Beta, Pi, Dom_table, Dom_spectral, & + Real (Kind=Kind(0.d0)), private :: Delta, Delta2, OM_st_1, Om_en_1, DeltaXMAX, Beta, Dom_table, Dom_spectral, & & Dx_spectral, Dx_table Real (Kind=Kind(0.d0)), allocatable, private :: XQMC1(:) Integer, allocatable, private :: Phim1_func(:), Phi_func(:) @@ -88,10 +88,9 @@ Subroutine MaxEnt_stoch(XQMC, Xtau, COV,Xmom1, XKER, Back_Trans_Aom, Beta_1, Alp Real (Kind=Kind(0.d0)), allocatable :: F_A_m(:), F_A_e(:) - Pi = acos(-1.d0) NDis = Ndis_1 - DeltaXMAX = 0.01 - delta = 0.001 + DeltaXMAX = 1.d-2 + delta = 1.d-3 delta2 = delta*delta Ngamma = Ngamma_1 Beta = Beta_1 ! Physical temperature for calculation of the kernel. @@ -459,7 +458,7 @@ Subroutine Set_default_table(Default, Default_table, Xmom1) Integer :: nw, nw1, nw_d, nx Real (Kind= Kind(0.d0)) :: om, dom, a, b, x, x1, f1,f2 - Logical :: Test=.false. + Logical, parameter :: Test=.false. dom = (Om_en_1 - Om_st_1)/dble(Ndis) @@ -506,7 +505,7 @@ Subroutine Set_default_table(Default, Default_table, Xmom1) nx = Int(x1/Dx_table) + 1 Phi_func(nw) = nx enddo - If (Test) then + If (Test) then ! Check the Phim1 function do nx = 1,Size(Default_table,1) x = dble(nx)*dx_table @@ -561,7 +560,7 @@ Real (Kind=Kind(0.d0)) Function Phim1(x) Implicit None ! Flat Default with sum xmom1. ! D(om) = xmom1/(Om_en_1 - Om_st_1) - Real (Kind=Kind(0.d0)) :: x, test + Real (Kind=Kind(0.d0)) :: x Integer :: nx nx = int(x/Dx_table) + 1 @@ -604,7 +603,6 @@ Subroutine Sum_Xn_Boxes(Xn_m,Xn) Implicit none Real (Kind=Kind(0.d0)), Dimension(:,:) :: Xn Real (Kind=Kind(0.d0)), Dimension(:) :: Xn_m - Real (Kind=Kind(0.d0)) :: om Integer :: nd, ng do ng = 1,Ngamma diff --git a/Libraries/Modules/mpi_shared_mem_mod_v2.F90 b/Libraries/Modules/mpi_shared_mem_mod_v2.F90 index 7b8dfba11..a15fd861a 100644 --- a/Libraries/Modules/mpi_shared_mem_mod_v2.F90 +++ b/Libraries/Modules/mpi_shared_mem_mod_v2.F90 @@ -101,8 +101,7 @@ subroutine mpi_shared_memory_init(mpi_communicator, chunk_size) real(Kind=Kind(0.d0)), intent(in) :: chunk_size #ifdef MPI - integer :: ierr, tmp_int, status - real(Kind=Kind(0.d0)) :: dummy_real_dp + integer :: ierr chunk_size_gb=chunk_size if (chunk_size_gb > 0) then @@ -592,9 +591,9 @@ end subroutine allocate_shared_memory_4Dcmplx subroutine deallocate_all_shared_memory() Implicit none - integer :: ierr, i, mpi_win_loc #ifdef MPI + integer :: ierr, i, mpi_win_loc external :: MPI_Win_free ! This seems to be required by gfortran10 with OpenMPI on Fedora33 (should be part of MPI module) do i=1,num_chunks_real diff --git a/Libraries/Modules/natural_constants_mod.F90 b/Libraries/Modules/natural_constants_mod.F90 index 505620035..f5ebee2ad 100644 --- a/Libraries/Modules/natural_constants_mod.F90 +++ b/Libraries/Modules/natural_constants_mod.F90 @@ -1,16 +1,18 @@ Module Natural_Constants - Real (Kind=Kind(0.d0)) :: eV, amu, Ang, hbar, pi + ! Mathematical constants + Real (Kind=Kind(0.d0)), parameter :: pi = acos(-1.d0) + Real (Kind=Kind(0.d0)), parameter :: twopi = 2.d0 * pi - contains - - subroutine Set_NC - - pi = acos(-1.d0) - eV = (1.0/6.24150974) *( 10.0**(-18) ) - amu = 1.66053886 * (10.0**(-27)) - Ang = 10.0**(-10) - hbar = 6.6260755*(10.0**(-34))/(2.0*pi) + ! Tolerance constants + Real (Kind=Kind(0.d0)), parameter :: Eps_machine = 1.D-12 ! Near machine precision: eigenvalue zero-detection, floating-point identity + Real (Kind=Kind(0.d0)), parameter :: Eps_small = 1.D-10 ! Negligibility: parameter on/off gating, singularity protection + Real (Kind=Kind(0.d0)), parameter :: Eps_convergence = 1.D-6 ! Iterative convergence, coarse k-space detection + + ! Physical constants + Real (Kind=Kind(0.d0)), parameter :: eV = (1.0/6.24150974) *( 10.0**(-18) ) + Real (Kind=Kind(0.d0)), parameter :: amu = 1.66053886 * (10.0**(-27)) + Real (Kind=Kind(0.d0)), parameter :: Ang = 10.0**(-10) + Real (Kind=Kind(0.d0)), parameter :: hbar = 6.6260755*(10.0**(-34))/(2.0*pi) - end subroutine Set_NC end Module Natural_Constants diff --git a/Libraries/Modules/random_wrap_mod.F90 b/Libraries/Modules/random_wrap_mod.F90 index 6701998cd..ac1e81328 100644 --- a/Libraries/Modules/random_wrap_mod.F90 +++ b/Libraries/Modules/random_wrap_mod.F90 @@ -144,14 +144,14 @@ end function ranf_wrap real (kind=kind(0.D0)) function rang_wrap(iq) ! Random variable according to the distribution: exp(-x**2/2)/(sqrt(2*3.1415927)) + Use Natural_Constants, only: twopi Implicit none integer, optional :: iq - real (Kind=kind(0.D0)) :: pi, ranmod, theta + real (Kind=kind(0.D0)) :: ranmod, theta - PI = 3.1415926536D0 RANMOD = SQRT(-2.D0 * log(RANF_Wrap(iq))) - THETA = 2.D0 * PI * RANF_wrap(iq) + THETA = twopi * RANF_wrap(iq) rang_wrap = RANMOD * COS(THETA) end function rang_wrap diff --git a/Prog/DynamicMatrixArray_mod.F90 b/Prog/DynamicMatrixArray_mod.F90 index c0ac8d20f..9cac46842 100644 --- a/Prog/DynamicMatrixArray_mod.F90 +++ b/Prog/DynamicMatrixArray_mod.F90 @@ -100,12 +100,14 @@ subroutine DynamicMatrixArray_dealloc(this) !> If out of space the vector grows. !> !> @param[inout] this the vector -!> @param[in] itm the object that we like to store a pointer to. +!> @param[inout] itm the object that we like to store a pointer to. ! !-------------------------------------------------------------------- subroutine DynamicMatrixArray_pushback(this, itm) class(DynamicMatrixArray) :: this - class(ContainerElementBase), intent(in), target :: itm !Type(...) has to match exactly, class(...) allows for polymorphism + ! INTENT(INOUT) required: flang rejects INTENT(IN) when a pointer to the dummy + ! argument is stored (this%data(...)%dat => itm). itm is not modified in this routine. + class(ContainerElementBase), intent(inout), target :: itm !Type(...) has to match exactly, class(...) allows for polymorphism type(OpTbasePtrWrapper), allocatable, dimension(:) :: temp integer :: i diff --git a/Prog/Fields_mod.F90 b/Prog/Fields_mod.F90 index a63164560..80d7c210a 100644 --- a/Prog/Fields_mod.F90 +++ b/Prog/Fields_mod.F90 @@ -343,11 +343,13 @@ Subroutine Fields_in(this,Group_Comm,Initial_field) Complex (Kind=Kind(0.d0)), Dimension(:,:), Optional :: Initial_field ! LOCAL - Integer :: I, I1, IERR, SEED_IN, K, NT - Real (Kind=Kind(0.d0) ) :: X + Integer :: IERR, SEED_IN +#ifdef MPI + Integer :: I +#endif Integer, DIMENSION(:), ALLOCATABLE :: SEED_VEC Logical :: LCONF, LCONF_H5 - Character (LEN=64) :: FILE_SR, FILE_TG, FILE_seeds, FILE_info, File1, FILE_TG_H5, File1_h5 + Character (LEN=64) :: FILE_TG, FILE_seeds, FILE_info, File1, FILE_TG_H5, File1_h5 #ifdef MPI INTEGER :: STATUS(MPI_STATUS_SIZE), irank_g, isize_g, igroup, ISIZE, IRANK @@ -524,7 +526,7 @@ Subroutine Fields_test(this) Class (Fields), INTENT(INOUT) :: this - Integer :: nt, I, I1 + Integer :: nt, I !Write(6,*) "Fields_set", size(this%f,1), size(this%f,2) Do nt = 1,size(this%f,2) diff --git a/Prog/Global_mod.F90 b/Prog/Global_mod.F90 index 46f47e338..f8a3d8f75 100644 --- a/Prog/Global_mod.F90 +++ b/Prog/Global_mod.F90 @@ -123,10 +123,10 @@ Subroutine Exchange_Step(Phase,GR, udvr, udvl, Stab_nt, udvst, N_exchange_steps, !> Local variables. - Integer :: NST, NSTM, NF, nf_eff, NT, NT1, NVAR,N, N1,N2, I, NC, I_Partner, n_step, N_count, N_part + Integer :: NST, NSTM, NF, nf_eff, NT, NT1, NVAR,N, N1,N2, I, NC, n_step, N_count Type (Fields) :: nsigma_old - Real (Kind=Kind(0.d0)) :: log_T0_Proposal_ratio, Weight, Weight1, delta_S0_log, exp_delta_S0 - Complex (Kind=Kind(0.d0)) :: Z_ONE = cmplx(1.d0, 0.d0, kind(0.D0)), Z, Ratiotot, Ratiotot_p, Phase_old, Phase_new + Real (Kind=Kind(0.d0)) :: log_T0_Proposal_ratio, Weight, delta_S0_log + Complex (Kind=Kind(0.d0)) :: Z, Ratiotot, Phase_old, Phase_new Real (Kind=Kind(0.d0)), allocatable :: Det_vec_old(:,:), Det_vec_new(:,:) Complex (Kind=Kind(0.d0)), allocatable :: Phase_Det_new(:), Phase_Det_old(:) Complex (Kind=Kind(0.d0)) :: Ratio(2), Ratio_p(2),Phase_array(N_FL) @@ -464,9 +464,9 @@ Subroutine Global_Updates(Phase,GR, udvr, udvl, Stab_nt, udvst,N_Global) ! Local variables. - Integer :: NST, NSTM, NF, NT, NT1, NVAR,N, N1,N2, I, NC, N_part,j, nf_eff + Integer :: NST, NSTM, NF, NT, NT1, NVAR,N, N1,N2, I, NC, nf_eff Real (Kind=Kind(0.d0)) :: log_T0_Proposal_ratio, Weight - Complex (Kind=Kind(0.d0)) :: Z_ONE = cmplx(1.d0, 0.d0, kind(0.D0)), Z, Ratiotot, Phase_old, Phase_new + Complex (Kind=Kind(0.d0)) :: Z, Ratiotot, Phase_old, Phase_new Complex (Kind=Kind(0.d0)), allocatable :: Det_vec_test(:,:), Phase_Det_new(:), Phase_Det_old(:) Real (Kind=Kind(0.d0)), allocatable :: Det_vec_old(:,:), Det_vec_new(:,:) Type (Fields) :: nsigma_old @@ -684,8 +684,8 @@ Complex (Kind=Kind(0.d0)) Function Compute_Ratio_Global(Phase_Det_old, Phase_Det ! Local Integer :: Nf, i, nt, nf_eff - Complex (Kind=Kind(0.d0)) :: Z, Z1, Ratio_1_array(N_FL), Ratio_2_array(N_FL), g_loc - Real (Kind=Kind(0.d0)) :: X, Ratio_2, delta, log_delta + Complex (Kind=Kind(0.d0)) :: Z, Ratio_1_array(N_FL), Ratio_2_array(N_FL), g_loc + Real (Kind=Kind(0.d0)) :: X, log_delta Ratio = cmplx(0.d0,0.d0,kind(0.d0)) Ratio_2_array = 0.d0 diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 index 96fed7bc3..6c1068ac8 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 @@ -197,15 +197,14 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: nf, unit_info Character (len=64) :: file_info Logical :: Half_filling #ifdef MPI - Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) + Integer :: ierr, Isize, Irank, irank_g, isize_g, igroup 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) @@ -386,12 +385,13 @@ end Subroutine Ham_Hop Subroutine Ham_Trial() Use Predefined_Trial + Use Natural_Constants, only: pi Implicit none Integer :: nf, Ix, Iy, I, n Real (Kind=Kind(0.d0)), allocatable :: H0(:,:), U0(:,:), E0(:) - Real (Kind=Kind(0.d0)) :: Pi = acos(-1.d0), Delta = 0.01d0 + Real (Kind=Kind(0.d0)), parameter :: Delta = 0.01d0 Allocate(WF_L(N_FL),WF_R(N_FL)) do nf=1,N_FL @@ -597,10 +597,8 @@ subroutine Obser(GR,Phase,Ntau, Mc_step_weight) !Local Complex (Kind=Kind(0.d0)), allocatable :: GRC(:,:,:) - Complex (Kind=Kind(0.d0)) :: ZK - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY, ZDen + Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, ZP,ZS, ZZ, ZXY, ZDen Integer :: I,J, imj, nf, Ix, Iy - Real (Kind=Kind(0.d0)) :: X ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -723,9 +721,8 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE, Mc_step_weight ) Real (Kind=Kind(0.d0)), INTENT(IN) :: Mc_step_weight !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY, ZDEN - Real (Kind=Kind(0.d0)) :: X - Integer :: IMJ, I, J, I1, J1, no_I, no_J + Complex (Kind=Kind(0.d0)) :: ZP, ZS, ZZ, ZXY, ZDEN + Integer :: IMJ, I, J ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -844,7 +841,7 @@ Subroutine GRT_reconstruction(GT0, G0T) Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GT0(Ndim,Ndim,N_FL), G0T(Ndim,Ndim,N_FL) Integer :: I,J,imj - real (kind=kind(0.d0)) :: X, ZZ + real (kind=kind(0.d0)) :: X If (Ham_U >= 0.d0) then Do J = 1,Latt%N diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 index b94051eda..7f1ab23e5 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 @@ -129,6 +129,7 @@ Use Fields_mod Use Predefined_Hoppings Use LRC_Mod + Use Natural_Constants, only: Eps_small Implicit none @@ -210,7 +211,7 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: nf, unit_info Character (len=64) :: file_info Logical :: toggle @@ -221,8 +222,7 @@ Subroutine Ham_Set ! Simulation type --> Finite T or Projection Symmetrize Trotter. #ifdef MPI - Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) + Integer :: ierr, Isize, Irank, irank_g, isize_g, igroup 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) @@ -256,7 +256,7 @@ Subroutine Ham_Set Write(error_unit,*) 'Ham_Set: If N_FL = 2, N_SUN has to be even' CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) endif - If (abs(ham_h0) >= 1.D-10 .and. N_FL /= 2 ) then + If (abs(ham_h0) >= Eps_small .and. N_FL /= 2 ) then Write(error_unit,*) 'Ham_Set: Set N_fl=2 if you want to use the pinning field' CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) endif @@ -317,7 +317,7 @@ Subroutine Ham_Set Write(unit_info,*) 'N_FL : ', N_FL Write(unit_info,*) 't : ', Ham_T Write(unit_info,*) 'Ham_U : ', Ham_U - if (abs(ham_h0) >= 1.D-10 ) Write(unit_info,*) 'Pinning field : ', ham_h0 + if (abs(ham_h0) >= Eps_small ) Write(unit_info,*) 'Pinning field : ', ham_h0 if (Index(str_to_upper(Lattice_type),'BILAYER') > 0 ) then Write(unit_info,*) 't2 : ', Ham_T2 Write(unit_info,*) 'Ham_U2 : ', Ham_U2 @@ -370,7 +370,7 @@ Subroutine Ham_Hop Integer, allocatable :: N_Phi_vec(:) ! Use predefined stuctures or set your own hopping - Integer :: n,nth, i + Integer :: i ! Indices of pinned vertices. Shape [N_pinned_vertices, 2] Integer, allocatable :: pinned_vertices(:,:) ! Factor, by which the vertex matrix elements will get multiplied. Shape [N_pinned_vertices, N_FL] @@ -459,7 +459,7 @@ Subroutine Ham_Trial() Implicit none - Integer :: N_part, nf + Integer :: N_part ! Use predefined stuctures or set your own Trial wave function N_part = Ndim/2 Call Predefined_TrialWaveFunction(Lattice_type ,Ndim, List,Invlist,Latt, Latt_unit, & @@ -479,8 +479,8 @@ Subroutine Ham_V Use Predefined_Int Implicit none - Integer :: nf, I, I1, I2, nc, J, no, N_ops - Real (Kind=Kind(0.d0)) :: X, Zero = 1.D-10 + Integer :: nf, I, I1, nc, no, N_ops + Real (Kind=Kind(0.d0)) :: Zero = Eps_small Real (Kind=Kind(0.d0)), allocatable :: Ham_U_vec(:) @@ -578,7 +578,7 @@ Subroutine Alloc_obs(Ltau) Call Obser_Vec_make(Obs_scal(I),N,Filename) enddo ! Local quantities - If (abs(ham_h0) >= 1.D-10 ) then + If (abs(ham_h0) >= Eps_small ) then Allocate ( Obs_local(1) ) Filename = "SpinZ" Channel = "--" @@ -705,10 +705,8 @@ subroutine Obser(GR,Phase,Ntau, Mc_step_weight) !Local Complex (Kind=Kind(0.d0)), allocatable :: GRC(:,:,:) - Complex (Kind=Kind(0.d0)) :: ZK - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY - Integer :: I,J, imj, nf, dec, I1, J1, no_I, no_J,n - Real (Kind=Kind(0.d0)) :: X + Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, ZP,ZS + Integer :: I,J, nf, dec, I1, no_I ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -773,7 +771,7 @@ subroutine Obser(GR,Phase,Ntau, Mc_step_weight) Obs_scal(4)%Obs_vec(1) = Obs_scal(4)%Obs_vec(1) + (Zkin + Zpot)*ZP*ZS - If (abs(ham_h0) >= 1.D-10 ) then + If (abs(ham_h0) >= Eps_small ) then ! Compute local observables. Do I = 1,Size(Obs_local,1) Obs_local(I)%N = Obs_local(I)%N + 1 @@ -837,9 +835,7 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE, Mc_step_weight) Real (Kind=Kind(0.d0)), INTENT(IN) :: Mc_step_weight !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY - Real (Kind=Kind(0.d0)) :: X - Integer :: IMJ, I, J, I1, J1, no_I, no_J + Complex (Kind=Kind(0.d0)) :: ZP, ZS ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -880,7 +876,6 @@ Real (Kind=Kind(0.d0)) function S0(n,nt,Hs_new) !> New local field on time slice nt and operator index n Complex (Kind=Kind(0.d0)), Intent(In) :: Hs_new - Integer :: nt1,I if (Continuous) then S0 = exp( (-real(Hs_new,kind(0.d0))**2 + real(nsigma%f(n,nt),Kind(0.d0))**2 ) /2.d0 ) diff --git a/Prog/Hamiltonians/Hamiltonian_Kondo_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Kondo_smod.F90 index 2d66e53c3..984658484 100644 --- a/Prog/Hamiltonians/Hamiltonian_Kondo_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Kondo_smod.F90 @@ -127,6 +127,7 @@ Use Fields_mod Use Predefined_Hoppings Use LRC_Mod + Use Natural_Constants, only: Eps_small Implicit none @@ -207,7 +208,7 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: nf, unit_info Character (len=64) :: file_info @@ -218,8 +219,7 @@ Subroutine Ham_Set #ifdef MPI - Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) + Integer :: ierr, Isize, Irank, irank_g, isize_g, igroup 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) @@ -401,7 +401,6 @@ Subroutine Ham_Hop Integer, allocatable :: N_Phi_vec(:) ! Use predefined stuctures or set your own hopping - Integer :: n,nth Allocate (Ham_T_vec(N_FL), Ham_T2_vec(N_FL), Ham_Tperp_vec(N_FL), Ham_Chem_vec(N_FL), Phi_X_vec(N_FL), Phi_Y_vec(N_FL),& & N_Phi_vec(N_FL), Ham_Lambda_vec(N_FL) ) @@ -447,7 +446,7 @@ Subroutine Ham_Trial() Implicit none - Integer :: N_part, nf + Integer :: N_part ! Use predefined stuctures or set your own Trial wave function N_part = Ndim/2 Call Predefined_TrialWaveFunction(Lattice_type ,Ndim, List,Invlist,Latt, Latt_unit, & @@ -467,9 +466,8 @@ Subroutine Ham_V Use Predefined_Int Implicit none - Integer :: nf, I, I1, I2, nc, no, N_ops - Real (Kind=Kind(0.d0)) :: X, Zero=1.D-10 - Real (Kind=Kind(0.d0)), allocatable :: Ham_U_vec(:) + Integer :: I, I1, I2, nc, no, N_ops + Real (Kind=Kind(0.d0)), parameter :: Zero=Eps_small N_ops = 0 @@ -661,10 +659,8 @@ subroutine Obser(GR,Phase,Ntau, Mc_step_weight) !Local Complex (Kind=Kind(0.d0)), allocatable :: GRC(:,:,:) - Complex (Kind=Kind(0.d0)) :: ZK - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, Zhubc, ZCon, ZJ, Z, ZP,ZS, ZZ, ZXY - Integer :: I,J, no, n, I_c,I_f, nf, J_c, J_f, no_I, no_J, imj - Real (Kind=Kind(0.d0)) :: X + Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, Zhubc, ZCon, ZJ, Z, ZP,ZS + Integer :: I,J, no, I_c,I_f, nf, J_c, J_f, no_I, no_J, imj ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -807,8 +803,7 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE,Mc_step_weight) !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY - Real (Kind=Kind(0.d0)) :: X + Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS Integer :: IMJ, I_c, I_f, J_c, J_f, I,J, no_I, no_J ZP = PHASE/Real(Phase, kind(0.D0)) diff --git a/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 b/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 index c60b8cfc5..5a93a62e3 100644 --- a/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 @@ -224,7 +224,6 @@ Subroutine Ham_Set #ifdef MPI Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) 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) @@ -364,7 +363,6 @@ Subroutine Ham_Hop Integer, allocatable :: N_Phi_vec(:) ! Use predefined stuctures or set your own hopping - Integer :: n,nth Allocate (Ham_T_vec(N_FL), Ham_T2_vec(N_FL), Ham_Tperp_vec(N_FL), Ham_Chem_vec(N_FL), Phi_X_vec(N_FL), Phi_Y_vec(N_FL),& & N_Phi_vec(N_FL), Ham_Lambda_vec(N_FL) ) @@ -492,7 +490,7 @@ Real (Kind=Kind(0.d0)) function S0(n,nt,Hs_new) !> New local field on time slice nt and operator index n Complex (Kind=Kind(0.d0)), Intent(In) :: Hs_new - Integer :: nt1, I, ns + Integer :: I Real (Kind=Kind(0.d0)) :: Y Real (Kind=Kind(0.d0)), Allocatable :: V(:) @@ -617,10 +615,8 @@ subroutine Obser(GR,Phase,Ntau,Mc_step_weight) !Local Complex (Kind=Kind(0.d0)), allocatable :: GRC(:,:,:) - Complex (Kind=Kind(0.d0)) :: ZK - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY - Integer :: I,J, imj, nf, dec, I1, J1, no_I, no_J,n - Real (Kind=Kind(0.d0)) :: X + Complex (Kind=Kind(0.d0)) :: Z, Zrho, Zkin, ZPot, ZP,ZS + Integer :: I,J, nf ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -718,9 +714,7 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE, Mc_step_weight) !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY - Real (Kind=Kind(0.d0)) :: X - Integer :: IMJ, I, J, I1, J1, no_I, no_J + Complex (Kind=Kind(0.d0)) :: ZP, ZS ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -784,8 +778,7 @@ Subroutine Global_move_tau(T0_Proposal_ratio, S0_ratio, & ! Local - Integer :: n_op, n, ns - Real (Kind=Kind(0.d0)) :: T0_proposal + Integer :: n, ns Real (Kind=Kind(0.d0)), allocatable :: V(:), V1(:) n = Size(nsigma%f,1) diff --git a/Prog/Hamiltonians/Hamiltonian_Spin_Peierls_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Spin_Peierls_smod.F90 index cea12f4f7..466cc51e2 100644 --- a/Prog/Hamiltonians/Hamiltonian_Spin_Peierls_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Spin_Peierls_smod.F90 @@ -130,6 +130,7 @@ Use Predefined_Hoppings Use LRC_Mod Use Predefined_Lattices + Use Natural_Constants, only: Eps_small Implicit none @@ -218,13 +219,12 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: unit_info Character (len=64) :: file_info #ifdef MPI - Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) + Integer :: ierr, Isize, Irank, irank_g, isize_g, igroup 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) @@ -295,7 +295,7 @@ Subroutine Ham_Set CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) endif - If ( Ham_h <= 1.D-8 ) SU2_Symm = .true. + If ( Ham_h <= Eps_small ) SU2_Symm = .true. If (SU2_Symm .and. N_FL .ne. 1 .and. N_SUN .ne. 2 ) then write(error_unit,*) " SU(2) symmetry is present " write(error_unit,*) " N_FL has to be equal to 1 " @@ -335,8 +335,6 @@ end Subroutine Ham_Set Subroutine Ham_Latt Implicit none - Real (Kind=Kind(0.d0)) :: a1_p(2), a2_p(2), L1_p(2), L2_p(2) - Integer :: nc, no, I Call Predefined_Latt(Lattice_type, L1, L2, Ndim, List, Invlist, Latt, Latt_Unit ) @@ -399,9 +397,9 @@ Subroutine Ham_V Use Predefined_Int Implicit none - Integer :: nf, I, I1, I2, nc, nc1, J, N_op, Ix, Iy - Integer :: n, nst,nen, n_fam, no - Real (Kind=Kind(0.d0)) :: X, J_Heis, X_p(2) + Integer :: nf, I, I1, I2, nc, N_op, Ix, Iy + Integer :: n, nst, n_fam, no + Real (Kind=Kind(0.d0)) :: J_Heis, X_p(2) @@ -888,7 +886,7 @@ subroutine GR_reconstruction(GR) Implicit none Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) - Integer :: I,J,imj + Integer :: I,J real (kind=kind(0.d0)) :: XI,XJ, ZZ Do J = 1,Ndim @@ -924,8 +922,8 @@ Subroutine GRT_reconstruction(GT0, G0T) Implicit none Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GT0(Ndim,Ndim,N_FL), G0T(Ndim,Ndim,N_FL) - Integer :: I,J,imj - real (kind=kind(0.d0)) :: XI,XJ, ZZ + Integer :: I,J + real (kind=kind(0.d0)) :: XI,XJ Do J = 1,NDIM XJ = 1.d0 diff --git a/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 index 05a099c26..8f89e32e2 100644 --- a/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 @@ -54,6 +54,7 @@ Use Fields_mod Use Predefined_Hoppings Use Predefined_Obs + Use Natural_Constants, only: Eps_small Implicit none @@ -102,7 +103,7 @@ Type (Unit_cell), target :: Latt_unit Logical :: One_dimensional Integer, allocatable :: List(:,:), Invlist(:,:) ! For orbital structure of Unit cell - real (Kind=Kind(0.d0)) :: Zero = 1.D-10 + real (Kind=Kind(0.d0)) :: Zero = Eps_small !> Storage for the Ising action Real (Kind=Kind(0.d0)) :: DW_Ising_tau(-1:1), DW_Ising_Space(-1:1), DW_Ising_Flux(-1:1,-1:1) @@ -133,12 +134,12 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: nf, unit_info Character (len=64) :: file_info #ifdef MPI + Integer :: ierr Integer :: Isize, Irank, igroup, irank_g, isize_g - Integer :: STATUS(MPI_STATUS_SIZE) 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) @@ -314,8 +315,7 @@ Subroutine Ham_V Use Predefined_Int Implicit none - Integer :: nf, I, I1, I2, nc, nc1, J, N_Field_type, N_ops - Real (Kind=Kind(0.d0)) :: X + Integer :: nf, I, I1, nc, N_Field_type, N_ops N_ops = size(Field_list_inv,1) Allocate(Op_V(N_ops, N_FL)) @@ -380,12 +380,13 @@ end Subroutine Ham_V Subroutine Ham_Trial() Use Predefined_Trial + Use Natural_Constants, only: pi Implicit none Integer :: nf, Ix, Iy, I, n Real (Kind=Kind(0.d0)), allocatable :: H0(:,:), U0(:,:), E0(:) - Real (Kind=Kind(0.d0)) :: Pi = acos(-1.d0), Delta = 0.01d0 + Real (Kind=Kind(0.d0)), parameter :: Delta = 0.01d0 Allocate(WF_L(N_FL),WF_R(N_FL)) do nf=1,N_FL @@ -544,7 +545,7 @@ Subroutine Global_move_tau(T0_Proposal_ratio, S0_ratio, & !Local - Integer :: ns , nc, n_op, n_op1, ntau_p1, ntau_m1, I, n + Integer :: n_op, n_op1, ntau_p1, ntau_m1, I, n Integer, allocatable :: Isigma1(:),Isigma2(:),Isigma3(:) Real (Kind = Kind(0.d0)) :: S0_Matter, T0_Proposal @@ -663,7 +664,7 @@ Real (Kind=kind(0.d0)) Function Get_Delta_S0_global(Nsigma_old) !> Arguments type (Fields), Intent(IN) :: nsigma_old !> Local - Integer :: I,n,n1,n2,n3,n4,nt,nt1, nc_F, nc_J, nc_h_p, nc_h_m, n1_m, n4_m + Integer :: I,n1,n2,n3,n4,nt,nt1, nc_F, nc_J, nc_h_p, nc_h_m, n1_m, n4_m Real (Kind=kind(0.d0)) :: exp_delta_S0 @@ -739,9 +740,8 @@ Subroutine Setup_Ising_action_and_field_list ! This subroutine sets up lists and arrays so as to enable an ! an efficient calculation of S0(n,nt) - Integer :: nc, nth, n, n1, n2, n3, n4, I, I1, n_orientation, Ix, Iy, N_Field_type, N_Pos + Integer :: nc, n, n1, I, I1, n_orientation, Ix, Iy, N_Field_type, N_Pos Integer :: N_ops - Real (Kind=Kind(0.d0)) :: X_p(2) N_ops = 0 @@ -969,16 +969,12 @@ Subroutine Obser(GR,Phase,Ntau, Mc_step_weight) !Local - Complex (Kind=Kind(0.d0)) :: GRC(Ndim,Ndim,N_FL), ZK - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin_mat, ZPot_mat, Z, ZP,ZS, Z1, Z2, ZN - Complex (Kind=Kind(0.d0)) :: ZQ, ZSTAR, ZQT, ZQTT - Integer :: I,J, imj, nf, dec, I1, I2,I3,I4, J1,J2,J3,J4, no_I, no_J, iFlux_tot, & - & no, no1, ntau1, ntau2, L_Vison, L_Wilson, n, nx,ny - Real (Kind=Kind(0.d0)) :: X_ave, X, XI1,XI2,XI3,XI4, X_p(2) + Complex (Kind=Kind(0.d0)) :: GRC(Ndim,Ndim,N_FL) + Complex (Kind=Kind(0.d0)) :: Zrho, Z, ZP,ZS, Z1, ZN + Integer :: I,J, imj, nf, I1, J1, iFlux_tot, & + & no, ntau1 + Real (Kind=Kind(0.d0)) :: X_ave Integer, allocatable :: Isigma(:), Isigmap1(:) - Integer :: IB_x, IB_y, Ix, Iy - - Real (Kind=Kind(0.d0)) :: X_star_i, X_star_j, X_star_ij ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -1140,8 +1136,8 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE,Mc_step_weight) Real (Kind=Kind(0.d0)), INTENT(IN) :: Mc_step_weight !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS - Integer :: IMJ, I, J, I1, J1, no_I, no_J, NT1 + Complex (Kind=Kind(0.d0)) :: ZP, ZS + Integer :: NT1 NT1 = NT If (NT == 0 ) NT1 = LTROT @@ -1366,7 +1362,7 @@ Real (Kind=Kind(0.d0)) function tau_x(i,nt, Isigma, Isigmap1) Integer, Intent(IN) :: Isigma(:), Isigmap1(:) Real ( Kind =Kind(0.d0) ) :: X - Integer :: I3, I4, nt1 + Integer :: I3, I4 tau_x = 1.d0 If (Abs(Ham_T) > Zero ) then diff --git a/Prog/Hamiltonians/Hamiltonian_tV_smod.F90 b/Prog/Hamiltonians/Hamiltonian_tV_smod.F90 index 017f13474..3a6275ab7 100644 --- a/Prog/Hamiltonians/Hamiltonian_tV_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_tV_smod.F90 @@ -204,7 +204,7 @@ Subroutine Ham_Set #endif Implicit none - integer :: ierr, nf, unit_info + integer :: nf, unit_info Character (len=64) :: file_info ! L1, L2, Lattice_type, List(:,:), Invlist(:,:) --> Lattice information @@ -213,8 +213,7 @@ Subroutine Ham_Set ! Simulation type --> Finite T or Projection Symmetrize Trotter. #ifdef MPI - Integer :: Isize, Irank, irank_g, isize_g, igroup - Integer :: STATUS(MPI_STATUS_SIZE) + Integer :: ierr, Isize, Irank, irank_g, isize_g, igroup 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) @@ -324,7 +323,6 @@ Subroutine Ham_Hop Integer, allocatable :: N_Phi_vec(:) ! Use predefined stuctures or set your own hopping - Integer :: n,nth Allocate (Ham_T_vec(N_FL), Ham_T2_vec(N_FL), Ham_Tperp_vec(N_FL), Ham_Chem_vec(N_FL), Phi_X_vec(N_FL), Phi_Y_vec(N_FL),& & N_Phi_vec(N_FL), Ham_Lambda_vec(N_FL) ) @@ -389,7 +387,7 @@ Subroutine Ham_Trial() Implicit none - Integer :: N_part, nf + Integer :: N_part ! Use predefined stuctures or set your own Trial wave function N_part = Ndim/2 Call Predefined_TrialWaveFunction(Lattice_type ,Ndim, List,Invlist,Latt, Latt_unit, & @@ -415,9 +413,8 @@ Subroutine Ham_Vint Type (Hopping_Matrix_type), Allocatable :: Bond_Matrix(:) Integer :: I, J, I1, J1, no_I, no_J, nf - Integer :: n_1, n_2, Nb, n_f,l_f, n_l, N, nc + Integer :: n_1, n_2, Nb, n_f,l_f, N, nc Complex(Kind=Kind(0.d0)) :: Z - real(Kind=Kind(0.d0)) :: Zero = 1.0E-6 Allocate (Ham_V_vec(N_FL), Ham_V2_vec(N_FL), Ham_Vperp_vec(N_FL), Ham_Chem_vec(N_FL), Phi_X_vec(N_FL), Phi_Y_vec(N_FL),& & N_Phi_vec(N_FL), Ham_Lambda_vec(N_FL) ) @@ -635,10 +632,9 @@ subroutine Obser(GR,Phase,Ntau, Mc_step_weight) Real (Kind=Kind(0.d0)), INTENT(IN) :: Mc_step_weight !Local - Complex (Kind=Kind(0.d0)) :: GRC(Ndim,Ndim,N_FL), ZK, Zn, weight, delta - Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY, tmp - Integer :: I,J, imj, nf, dec, I1, J1, no_I, no_J,n, nf2, k, k1, l, l1 - Real (Kind=Kind(0.d0)) :: X + Complex (Kind=Kind(0.d0)) :: GRC(Ndim,Ndim,N_FL), Zn, weight, delta + Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, ZP,ZS, tmp + Integer :: I,J, nf, I1, J1,n, nf2, k, k1, l, l1 ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) @@ -760,9 +756,7 @@ Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE,Mc_step_weight) Real (Kind=Kind(0.d0)), INTENT(IN) :: Mc_step_weight !Locals - Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY - Real (Kind=Kind(0.d0)) :: X - Integer :: IMJ, I, J, I1, J1, no_I, no_J + Complex (Kind=Kind(0.d0)) :: ZP, ZS ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) diff --git a/Prog/Hamiltonians/LRC_mod.F90 b/Prog/Hamiltonians/LRC_mod.F90 index 836d7c588..abcabcfec 100644 --- a/Prog/Hamiltonians/LRC_mod.F90 +++ b/Prog/Hamiltonians/LRC_mod.F90 @@ -49,6 +49,7 @@ Module LRC_mod Use Lattices_v3 Use MyMats Use Random_wrap + Use Natural_Constants, only: Eps_small use iso_fortran_env, only: output_unit, error_unit Implicit none @@ -86,7 +87,7 @@ Real ( Kind=Kind(0.d0) ) function LRC_V_func(X_p, Uhub, alpha, d1) LRC_V_func = 0.d0 X = Xnorm(X_p) - if ( abs(X) < 1.D-10 ) then + if ( abs(X) < Eps_small ) then LRC_V_func = Uhub else LRC_V_func = Uhub*alpha*d1/X @@ -195,7 +196,7 @@ Subroutine LRC_Print(Latt, Latt_unit, list, invlist) Integer, intent(in), Dimension(:,:) :: List, Invlist ! Local - Integer :: I,J, no_J, Ju, no_I, Iu, I0, imj, Latt_dim + Integer :: I,J, no_J, Ju, no_I, Iu, Latt_dim Real (Kind=Kind(0.d0)), allocatable :: X_p(:), X0_p(:) Real (Kind=Kind(0.d0)), allocatable :: A1_p(:), A2_p(:), L1_p(:), L2_p(:) @@ -222,7 +223,6 @@ Subroutine LRC_Print(Latt, Latt_unit, list, invlist) Do Ju = 1, Latt%N do no_J = 1,Latt_unit%Norb J = invlist(Ju,no_J) - ImJ = Latt%imj(Iu,Ju) X_p(:) = dble(Latt%list(Iu,1))*A1_p(:) + dble(Latt%list(Iu,2))*A2_p(:) + & & Latt_unit%Orb_pos_p(no_i,:) - & & dble(Latt%list(Ju,1))*A1_p(:) - dble(Latt%list(Ju,2))*A2_p(:) - & @@ -294,11 +294,11 @@ Subroutine LRC_Set_VIJ(Latt, Latt_unit, Uhub, alpha, list, invlist) Integer, intent(in), Dimension(:,:) :: List, Invlist !Local - Integer :: I,J,no_i,no_j, n, m, no, imj, Latt_dim + Integer :: I,J,no_i,no_j, n, m, no, Latt_dim Real (Kind=Kind(0.d0)) ::d1, X, X_min, Xmean,Xmax, Xmax1 Real (Kind=Kind(0.d0)), allocatable :: M_Tmp(:,:), M_Tmp1(:,:), X_p(:), X0_p(:), X1_p(:) Real (Kind=Kind(0.d0)), allocatable :: A1_p(:), A2_p(:), L1_p(:), L2_p(:) - Logical :: L_test=.true. + Logical, parameter :: L_test=.true. Latt_dim = Size(Latt_unit%Orb_pos_p,2) Allocate ( X_p(Latt_dim), X0_p(Latt_dim), X1_p(Latt_dim), & @@ -358,7 +358,7 @@ Subroutine LRC_Set_VIJ(Latt, Latt_unit, Uhub, alpha, list, invlist) Do I = 1,size(E_int,1) !Write(25,*) E_int(I) - if ( E_int(i) < 1.D-10 ) then + if ( E_int(i) < Eps_small ) then Write(error_unit,*) 'LRC_Set_VIJ: V_int(i,j) is not positive definite ' CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) endif @@ -506,8 +506,8 @@ Subroutine LRC_draw_field(Percent_change, Dtau, A_old, A_new,N_SUN) Real (Kind=Kind(0.d0)), Intent(IN) , dimension(:) :: A_old Real (Kind=Kind(0.d0)), Intent(OUT), dimension(:) :: A_new - Integer :: n, n1, i, m - Real (Kind=Kind(0.d0)) :: X, Alpha, Beta + Integer :: n, n1, m + Real (Kind=Kind(0.d0)) :: Alpha, Beta M = Size(E_int,1) do n = 1,M diff --git a/Prog/Hamiltonians/Not_supported/Hamiltonian_Hub_Canonical_mod.F90 b/Prog/Hamiltonians/Not_supported/Hamiltonian_Hub_Canonical_mod.F90 index 951a536b0..2e14d34b6 100644 --- a/Prog/Hamiltonians/Not_supported/Hamiltonian_Hub_Canonical_mod.F90 +++ b/Prog/Hamiltonians/Not_supported/Hamiltonian_Hub_Canonical_mod.F90 @@ -168,6 +168,7 @@ end Subroutine Ham_Latt !=================================================================================== Subroutine Ham_hop + Use Natural_Constants, only: twopi Implicit none !Setup the hopping @@ -175,12 +176,10 @@ Subroutine Ham_hop ! e^{-dtau H_t} = Prod_{n=1}^{Ncheck} e^{-dtau_n H_{n,t}} Integer :: I, I1, I2, n, Ncheck,nc - Complex (Kind=8) :: Z - Real (Kind=8) :: Pi + Complex (Kind=kind(0.d0)) :: Z Ncheck = 1 - Pi = acos(-1.d0) - Z = exp( cmplx(0.d0,Phi_x*2.d0*pi/real(L1,kind(0.d0)),kind(0.d0))) + Z = exp( cmplx(0.d0,Phi_x*twopi/real(L1,kind(0.d0)),kind(0.d0))) allocate(Op_T(Ncheck,N_FL)) do n = 1,N_FL Do nc = 1,Ncheck diff --git a/Prog/Hamiltonians/Not_supported/Hamiltonian_KN_Kondo_chain.F90 b/Prog/Hamiltonians/Not_supported/Hamiltonian_KN_Kondo_chain.F90 index c63fe2c01..b756acac3 100644 --- a/Prog/Hamiltonians/Not_supported/Hamiltonian_KN_Kondo_chain.F90 +++ b/Prog/Hamiltonians/Not_supported/Hamiltonian_KN_Kondo_chain.F90 @@ -480,6 +480,7 @@ End Subroutine Ham_Latt !-------------------------------------------------------------------- Subroutine Ham_Hop + Use Natural_Constants, only: twopi Implicit none Real (Kind=Kind(0.d0)) :: X_p(2), Delta_K @@ -507,7 +508,7 @@ Subroutine Ham_Hop Enddo Enddo elseif (K_space) then - Delta_K = 2.d0*acos(-1.d0)/real(L,Kind=Kind(0.d0)) + Delta_K = twopi/real(L,Kind=Kind(0.d0)) allocate(Op_T(K*L*Latt%N,1)) nc = 0 Do i_0 = 1,Latt%N diff --git a/Prog/Hop_mod.F90 b/Prog/Hop_mod.F90 index 9785c2209..504726305 100644 --- a/Prog/Hop_mod.F90 +++ b/Prog/Hop_mod.F90 @@ -59,7 +59,6 @@ Module Hop_mod ! Private variables Type(DynamicMatrixArray), private, allocatable :: ExpOpT_vec(:) ! for now we have for simplicity for each flavour a vector Integer, private, save :: Ncheck - Real (Kind=Kind(0.d0)), private, save :: Zero Contains @@ -137,7 +136,6 @@ subroutine Hop_mod_init enddo enddo - Zero = 1.E-12 end subroutine Hop_mod_init !-------------------------------------------------------------------- diff --git a/Prog/Langevin_HMC_mod.F90 b/Prog/Langevin_HMC_mod.F90 index 83e8da2cb..276fa86e0 100644 --- a/Prog/Langevin_HMC_mod.F90 +++ b/Prog/Langevin_HMC_mod.F90 @@ -119,9 +119,8 @@ SUBROUTINE Langevin_HMC_Forces(this, Phase, GR, GR_Tilde, Test, udvr, udvl, Sta !Local - Integer :: NSTM, n, nf, nf_eff, NST, NTAU, nt, nt1, Ntau1, NVAR, N_Type, I, J + Integer :: NSTM, nf, nf_eff, NST, NTAU, nt1, Ntau1, NVAR Complex (Kind=Kind(0.d0)) :: Z, Z1, Phase_array(N_FL) - Real (Kind=Kind(0.d0)) :: spin NSTM = Size(Stab_nt,1) - 1 !Do n = 0,NSTM diff --git a/Prog/OpTTypes_mod.F90 b/Prog/OpTTypes_mod.F90 index d0b06a5ed..7bfc2d100 100644 --- a/Prog/OpTTypes_mod.F90 +++ b/Prog/OpTTypes_mod.F90 @@ -34,6 +34,7 @@ module OpTTypes_mod use ContainerElementBase_mod use Operator_mod use mat_subroutines + use Natural_Constants, only: Eps_machine implicit none private @@ -95,7 +96,7 @@ subroutine RealExpOpT_init(this, Op_T) Complex(kind=kind(0.d0)) :: cg Integer :: i, j - this%Zero = 1.E-12 + this%Zero = Eps_machine this%Ndim_hop = Op_T%N allocate(this%mat(this%Ndim_hop, this%Ndim_hop), this%invmat(this%Ndim_hop, this%Ndim_hop)) allocate(this%mat_1D2(this%Ndim_hop, this%Ndim_hop), this%invmat_1D2(this%Ndim_hop, this%Ndim_hop)) @@ -204,7 +205,7 @@ subroutine CmplxExpOpT_init(this, Op_T) Type(Operator), intent(in) :: Op_T Integer :: i, j - this%Zero = 1.E-12 + this%Zero = Eps_machine this%Ndim_hop = Op_T%N allocate(this%mat(this%Ndim_hop, this%Ndim_hop), this%invmat(this%Ndim_hop, this%Ndim_hop)) diff --git a/Prog/OpT_time_dependent.F90 b/Prog/OpT_time_dependent.F90 index 8fc4fa53f..cd2a38467 100644 --- a/Prog/OpT_time_dependent.F90 +++ b/Prog/OpT_time_dependent.F90 @@ -34,6 +34,7 @@ module OpT_time_dependent_mod use ContainerElementBase_mod use Operator_mod use mat_subroutines + use Natural_Constants, only: Eps_machine implicit none private @@ -74,7 +75,7 @@ subroutine OpT_time_dependent_init(this, Op_T, symm) Type(Operator), intent(in) :: Op_T logical, intent(in) :: symm - this%Zero = 1.E-12 + this%Zero = Eps_machine this%Ndim_hop = Op_T%N this%P = Op_T%P ! copy all data locally to be consistent and less error prone this%U = Op_T%U diff --git a/Prog/Operator_mod.F90 b/Prog/Operator_mod.F90 index 9ef790180..81388d0bc 100644 --- a/Prog/Operator_mod.F90 +++ b/Prog/Operator_mod.F90 @@ -47,6 +47,7 @@ Module Operator_mod Use Fields_mod Use runtime_error_mod Use iso_fortran_env, only: error_unit + Use Natural_Constants, only: Eps_machine, Eps_small Implicit none @@ -261,10 +262,13 @@ subroutine Op_set(Op) Complex (Kind=Kind(0.d0)), allocatable :: U(:,:), TMP(:, :) Real (Kind=Kind(0.d0)), allocatable :: E(:) - Real (Kind=Kind(0.d0)) :: Zero = 1.D-9 !, Phi(-2:2) - Real (Kind=Kind(0.d0)) :: herm_tol = 1.D-12 + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small + Real (Kind=Kind(0.d0)) :: herm_tol = Eps_machine Real (Kind=Kind(0.d0)) :: herm_dev - Integer :: N, I, J, np,nz, noderank, arrayshape2d(2), arrayshape(3), ierr + Integer :: N, I, J, np,nz, noderank, arrayshape2d(2), arrayshape(3) +#ifdef MPI + Integer :: ierr +#endif Complex (Kind=Kind(0.d0)) :: Z Type (Fields) :: nsigma_single @@ -571,7 +575,7 @@ subroutine Op_mmultL(Mat,Op,HS_Field,cop, nt, sign) g_loc = OP%g if (op%g_t_alloc) g_loc = Op%g_t(nt) - if ( abs(g_loc) < 1.D-12 ) return + if ( abs(g_loc) < Eps_machine ) return if ( op%type < 3 ) then sp = sign*nint(Real(HS_Field)) @@ -665,7 +669,7 @@ subroutine Op_mmultR(Mat,Op,HS_Field,cop,nt) ! quick return if possible g_loc = Op%g if (op%g_t_alloc) g_loc = Op%g_t(nt) - if ( abs(g_loc) < 1.D-12 ) return + if ( abs(g_loc) < Eps_machine ) return if ( op%type < 3 ) then sp = nint(Real(Hs_Field)) diff --git a/Prog/Predefined_Hop_mod.F90 b/Prog/Predefined_Hop_mod.F90 index e89f8e5b1..742cc863a 100644 --- a/Prog/Predefined_Hop_mod.F90 +++ b/Prog/Predefined_Hop_mod.F90 @@ -48,6 +48,7 @@ Module Predefined_Hoppings Use Operator_mod Use WaveFunction_mod Use MyMats + Use Natural_Constants, only: twopi, Eps_small use iso_fortran_env, only: output_unit, error_unit use Hamiltonian_main Implicit none @@ -128,7 +129,7 @@ Integer function inquire_hop(this) Real (Kind=Kind(0.d0)) :: Xmax_loc, Xmax_hop, Zero, X Integer :: nc, nf - Zero = 1.D-10 + Zero = Eps_small Xmax_loc = 0.d0 Xmax_hop = 0.d0 @@ -269,8 +270,9 @@ Subroutine Set_Default_hopping_parameters_square(this, Ham_T_vec, Ham_Chem_vec, ! Local - Integer :: nf,N_Bonds, nc, I, I1 - Real (Kind = Kind(0.d0) ) :: Zero = 1.0E-8, Ham_T_max + Integer :: nf, nc, I + Real (Kind = Kind(0.d0) ), parameter :: Zero = Eps_small + Real (Kind = Kind(0.d0) ) :: Ham_T_max Real (Kind = Kind(0.d0) ), allocatable :: Ham_T_perp_vec(:) @@ -385,13 +387,14 @@ Subroutine Set_Default_hopping_parameters_triangular(this, Ham_T_vec, Ham_Chem_v ! Local - Integer :: nf,N_Bonds, nc, I, I1 - Real (Kind = Kind(0.d0) ) :: Zero = 1.0E-8, Ham_T_max, x_p(2) + Integer :: nf, nc, I + Real (Kind = Kind(0.d0) ), parameter :: Zero = Eps_small + Real (Kind = Kind(0.d0) ) :: Ham_T_max !Write(6,*) Iscalar(Latt%L1_p,Latt%BZ1_p)/(2.d0*acos(-1.d0)),Iscalar(Latt%L2_p,Latt%BZ2_p)/(2.d0*acos(-1.d0)) - If ( mod(nint(Iscalar(Latt%L1_p,Latt%BZ1_p)/(2.d0*acos(-1.d0))),2) /= 0 .or. & - & mod(nint(Iscalar(Latt%L2_p,Latt%BZ2_p)/(2.d0*acos(-1.d0))),2) /= 0 ) then + If ( mod(nint(Iscalar(Latt%L1_p,Latt%BZ1_p)/twopi),2) /= 0 .or. & + & mod(nint(Iscalar(Latt%L2_p,Latt%BZ2_p)/twopi),2) /= 0 ) then Write(error_unit,*) '*** For the triangular lattice, our implementation of the checkerborad ' Write(error_unit,*) 'decomposition requires even values of L_1 and L_2 ***' CALL Terminate_on_error(ERROR_GENERIC,__FILE__,__LINE__) @@ -520,8 +523,9 @@ Subroutine Set_Default_hopping_parameters_kagome(this, Ham_T_vec, Ham_Chem_vec, ! Local - Integer :: nf,N_Bonds, nc, I, I1 - Real (Kind = Kind(0.d0) ) :: Zero = 1.0E-8, Ham_T_max, x_p(2) + Integer :: nf, nc, I + Real (Kind = Kind(0.d0) ), parameter :: Zero = Eps_small + Real (Kind = Kind(0.d0) ) :: Ham_T_max !Write(6,*) Iscalar(Latt%L1_p,Latt%BZ1_p)/(2.d0*acos(-1.d0)),Iscalar(Latt%L2_p,Latt%BZ2_p)/(2.d0*acos(-1.d0)) @@ -642,8 +646,8 @@ Subroutine Set_Default_hopping_parameters_N_Leg_Ladder & ! Local - Integer :: nf,N_Bonds, nc, I, I1, n, no - Real (Kind=Kind(0.d0)) :: Zero = 1.0E-8 + Integer :: nf, nc, I, n, no + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small select case (Latt%N) @@ -836,8 +840,9 @@ Subroutine Set_Default_hopping_parameters_honeycomb(this,Ham_T_vec, Ham_Lambda_v Type(Unit_cell),Intent(in) :: Latt_unit ! Local - Integer :: nf,N_Bonds, nc, I, I1, n, no - Real (Kind=Kind(0.d0)) :: Zero = 1.0E-8, Ham_Lambda_Max + Integer :: nf, nc, I + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small + Real (Kind=Kind(0.d0)) :: Ham_Lambda_Max !Write(6,*) Ham_T_vec, Ham_Chem_vec Ham_Lambda_Max = 0.d0 @@ -928,9 +933,9 @@ Subroutine Set_Default_hopping_parameters_Bilayer_square(this,Ham_T1_vec,Ham_T2_ ! Local - Integer :: nf,N_Bonds, nc, I, I1, No_Shift, n, nb - Real (Kind=Kind(0.d0)) :: Zero = 1.0E-8 - Logical :: Test=.false. + Integer :: nf,N_Bonds, nc, I, No_Shift, n, nb + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small + Logical, parameter :: Test=.false. Real (Kind=Kind(0.d0)) :: Ham_T1_max, Ham_T2_max, Ham_Tperp_max Ham_T1_max = 0.d0 @@ -1238,9 +1243,9 @@ Subroutine Set_Default_hopping_parameters_Bilayer_honeycomb(this,Ham_T1_vec,Ham_ Real (Kind=Kind(0.d0)) :: Ham_T1_max, Ham_T2_max, Ham_Tperp_max ! Local - Integer :: nf,N_Bonds, nc, I, I1, No_Shift, n, nb, no - Real (Kind=Kind(0.d0)) :: Zero = 1.0E-8 - Logical :: Test=.false. + Integer :: nf,N_Bonds, nc, I, No_Shift, n, nb, no + Real (Kind=Kind(0.d0)), parameter :: Zero = Eps_small + Logical, parameter :: Test=.false. Ham_T1_max = 0.d0 Ham_T2_max = 0.d0 @@ -1533,7 +1538,7 @@ complex(Kind=Kind(0.d0)) function get_pinning_factor(I, J, N_pinned_vertices, pi else get_pinning_factor = cmplx(1.d0,0.d0, kind(0.d0)) endif - if( abs(get_pinning_factor - cmplx(1.d0,0.d0, kind(0.d0))) >= 1.d-10 ) pinning_notice_issued = .true. + if( abs(get_pinning_factor - cmplx(1.d0,0.d0, kind(0.d0))) >= Eps_small ) pinning_notice_issued = .true. end function get_pinning_factor !-------------------------------------------------------------------- @@ -1556,7 +1561,7 @@ complex(Kind=Kind(0.d0)) function get_pinning_offset(I, J, N_pinned_vertices, pi else get_pinning_offset = cmplx(0.d0,0.d0, kind(0.d0)) endif - if( abs(get_pinning_offset) >= 1.d-10 ) pinning_notice_issued = .true. + if( abs(get_pinning_offset) >= Eps_small ) pinning_notice_issued = .true. end function get_pinning_offset !-------------------------------------------------------------------- @@ -1640,12 +1645,12 @@ Subroutine Predefined_Hoppings_set_OPT(this,List,Invlist,Latt, Latt_unit, Dtau ! Local Integer :: Ndim, N_FL, N_Phi, I, J, I1, J1, no_I, no_J, nf - Integer :: n_1, n_2, Nb, n_f,l_f, n_l, N, nc, orb - Real (Kind=Kind(0.d0)) :: Ham_T, Ham_Chem, Phi_X, Phi_Y + Integer :: n_1, n_2, Nb, n_f, l_f, N, nc, orb + Real (Kind=Kind(0.d0)) :: Phi_X, Phi_Y Logical :: Bulk Complex(Kind=Kind(0.d0)) :: Z, Z1, Z2 - Integer :: N_pinned_vertices, i_pinned_vertex + Integer :: N_pinned_vertices N_FL = size(this,1) @@ -1862,12 +1867,12 @@ Subroutine Predefined_Hoppings_Compute_Kin(this,List,Invlist, Latt, Latt_unit, !Local - Integer :: Ndim, N_FL, N_Phi, I, J, I1, J1, no_I, no_J, nf - Integer :: n_1, n_2, Nb, n_f,l_f, n_l, N, nc - Real (Kind=Kind(0.d0)) :: Ham_T, Ham_Chem, Phi_X, Phi_Y + Integer :: N_FL, N_Phi, I, J, I1, J1, no_I, no_J, nf + Integer :: n_1, n_2, Nb + Real (Kind=Kind(0.d0)) :: Phi_X, Phi_Y Logical :: Bulk Complex(Kind=Kind(0.d0)) :: Z - Integer :: N_pinned_vertices, i_pinned_vertex + Integer :: N_pinned_vertices if( (present(pinned_vertices) .and. .not.(present(pinning_factor ).and. present(pinning_offset))) .or. & & (present(pinning_factor) .and. .not.(present(pinned_vertices).and. present(pinning_offset))) .or. & @@ -2067,8 +2072,8 @@ complex (Kind=kind(0.d0)) function Generic_hopping(i,no_i, del_1, del_2, no_j, !Local - Integer :: j, N1, N2,n - real (Kind=Kind(0.d0)) :: xj_p(2), xi_p(2), xjp_p(2), del_p(2), A_p(2), pi, XB_p(2), V, B, Zero, x_p(2), x1_p(2) + Integer :: N1, N2, n + real (Kind=Kind(0.d0)) :: xj_p(2), xi_p(2), xjp_p(2), del_p(2), A_p(2), XB_p(2), V, B, Zero, x_p(2), x1_p(2) Complex (Kind=Kind(0.d0)) :: Z_hop @@ -2098,7 +2103,6 @@ complex (Kind=kind(0.d0)) function Generic_hopping(i,no_i, del_1, del_2, no_j, del_p = xj_p - xi_p !Twist - pi = acos(-1.d0) A_p(:) = Flux_1 * Xnorm(Latt%a1_p) * latt%bZ1_p(:) / Xnorm(Latt%L1_p) + & & Flux_2 * Xnorm(Latt%a2_p) * latt%bZ2_p(:) / Xnorm(Latt%L2_p) @@ -2111,18 +2115,18 @@ complex (Kind=kind(0.d0)) function Generic_hopping(i,no_i, del_1, del_2, no_j, endif !Orbital magnetic field (Landau gauge) - Zero = 1.0E-8 + Zero = Eps_small V = abs(Latt%L1_p(1) * Latt%L2_p(2) - Latt%L1_p(2) * Latt%L2_p(1) ) If ( V > Zero ) then B = real(N_Phi,kind(0.d0))/V - Z_hop = Z_hop*exp(cmplx(0.d0, -2.d0*pi* B * del_p(1) * ( xj_p(2) + xi_p(2) )/2.d0,kind(0.d0) ) ) + Z_hop = Z_hop*exp(cmplx(0.d0, -twopi* B * del_p(1) * ( xj_p(2) + xi_p(2) )/2.d0,kind(0.d0) ) ) ! Boundary x_p = Real(N2,Kind(0.d0))*Latt%L2_p x1_p = Xjp_p + Real(N1,Kind(0.d0))*Latt%L1_p - Z_hop = Z_hop * exp(cmplx( 0.d0, -Chi(x_p, x1_p,B,pi),kind(0.d0))) + Z_hop = Z_hop * exp(cmplx( 0.d0, -Chi(x_p, x1_p,B),kind(0.d0))) x_p = Real(N1,Kind(0.d0))*Latt%L1_p x1_p = Xjp_p - Z_hop = Z_hop * exp(cmplx( 0.d0, -Chi(x_p, x1_p,B,pi),kind(0.d0))) + Z_hop = Z_hop * exp(cmplx( 0.d0, -Chi(x_p, x1_p,B),kind(0.d0))) endif Generic_hopping = Z_hop @@ -2137,12 +2141,13 @@ end function GENERIC_HOPPING !> Periodic boundary conditions for Landau gauge: c_{i+L} = e{-i Chi(L,i)} c_{i} !> !-------------------------------------------------------------------- - Real (Kind=kind(0.d0)) function Chi(L_p,X_p,B,pi) + Real (Kind=kind(0.d0)) function Chi(L_p,X_p,B) + Use Natural_Constants, only: twopi Implicit none - Real (Kind=Kind(0.d0)), Intent(In) :: L_p(2), X_p(2), B, pi + Real (Kind=Kind(0.d0)), Intent(In) :: L_p(2), X_p(2), B - Chi = - 2.d0 * pi *B * L_p(2) * X_p(1) + Chi = - twopi *B * L_p(2) * X_p(1) end function Chi end Module Predefined_Hoppings diff --git a/Prog/Predefined_Obs_mod.F90 b/Prog/Predefined_Obs_mod.F90 index f39fd922d..ef3923969 100644 --- a/Prog/Predefined_Obs_mod.F90 +++ b/Prog/Predefined_Obs_mod.F90 @@ -213,7 +213,7 @@ Subroutine Predefined_Obs_eq_Green_measure( Latt, Latt_unit, List, GR, GRC, N_S Type (Obser_Latt), Intent(inout) :: Obs ! Local - Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj,nf + Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj, nf Complex (Kind=Kind(0.d0)) :: Z If ( Size(List,1) .ne. Size(GR,1) .or. Size(List,2) .ne. 2 ) then @@ -343,7 +343,7 @@ Subroutine Predefined_Obs_tau_Green_measure( Latt, Latt_unit, List, NT, GT0,G0T, Type (Obser_Latt), Intent(inout) :: Obs ! Local - Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj,nf + Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj, nf Complex (Kind=Kind(0.d0)) :: Z @@ -406,7 +406,7 @@ Subroutine Predefined_Obs_tau_SpinSUN_measure( Latt, Latt_unit, List, NT, GT0,G0 Type (Obser_Latt), Intent(inout) :: Obs ! Local - Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj,nf + Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj Complex (Kind=Kind(0.d0)) :: Z If ( Obs%File_Latt .ne. "SpinZ" ) then @@ -467,7 +467,7 @@ Subroutine Predefined_Obs_tau_SpinMz_measure( Latt, Latt_unit, List, NT, GT0,G0T Type (Obser_Latt), Intent(inout) :: ObsZ, ObsXY, ObsXYZ ! Local - Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj,nf + Integer :: N_FL, I, I1, J, J1, no_I, no_J, imj Complex (Kind=Kind(0.d0)) :: ZZ, ZXY If ( Size(List,1) .ne. Size(GT0,1) .or. Size(List,2) .ne. 2 ) then @@ -615,7 +615,7 @@ Complex (Kind=Kind(0.d0)) function Predefined_Obs_V_Int(OP_Vint, GR, GRC, N_S ! Local Complex (Kind=Kind(0.d0)) :: Z, Z_tmp, ZC Integer :: N_FL, N, I, J, I1, J1, nf - Real (Kind=Kind(0.d0)) :: Zero=1.0D-16 + Real (Kind=Kind(0.d0)), parameter :: Zero=1.0D-16 If ( OP_Vint(1)%type .ne. 2 ) then Write(error_unit,*) 'Predefined_Obs_V_Int routine is defined fro tpye 2 vertices. ' diff --git a/Prog/Predefined_Trial_mod.F90 b/Prog/Predefined_Trial_mod.F90 index a3bf4698f..77d4d331e 100644 --- a/Prog/Predefined_Trial_mod.F90 +++ b/Prog/Predefined_Trial_mod.F90 @@ -114,6 +114,7 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, + Use Natural_Constants, only: pi, twopi Implicit none Character (len=64), Intent(IN) :: Lattice_type Integer, Intent(IN) :: Ndim, N_FL, N_part @@ -125,15 +126,16 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, Type(Operator), dimension(:,:), allocatable :: OP_tmp Type (Hopping_Matrix_type), allocatable :: Hopping_Matrix_tmp(:) - Real (Kind=Kind(0.d0)) :: Dtau, Ham_T, Ham_Chem, XB_X, XB_Y, Phi_X, Phi_Y, Dimer + Real (Kind=Kind(0.d0)) :: Dtau, Ham_T, Ham_Chem, Phi_X, Phi_Y Logical :: Checkerboard, Symm, Kekule_Trial Type (Lattice) :: Latt_Kekule Real (Kind=Kind(0.d0)) :: A1_p(2), A2_p(2), L1_p(2), L2_p(2), x_p(2),x1_p(2), hop(3), del_p(2) - Real (Kind=Kind(0.d0)) :: delta = 0.01, Ham_T1, Ham_T2, Ham_Tperp, dom, om + Real (Kind=Kind(0.d0)) :: delta = 1.d-2, Ham_T1, Ham_T2, Ham_Tperp, dom, om - Integer :: N, nf, I, I1, I2, nc, nc1, IK_u, I_u, J1, lp, J, N_Phi, den_file, Nom = 200 , nw - Logical :: Test=.false. , Bulk =.true. + Integer :: N, nf, I, I1, I2, nc, nc1, IK_u, I_u, J1, lp, N_Phi, den_file, Nom = 200 , nw + Logical, parameter :: Test=.false. + Logical :: Bulk =.true. Complex (Kind=Kind(0.d0)) :: Z_norm, Z Real (Kind=Kind(0.d0) ), allocatable :: Ham_T_vec(:), Ham_Tperp_vec(:), Ham_Chem_vec(:), Phi_X_vec(:), Phi_Y_vec(:),& @@ -302,7 +304,7 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, Enddo endif Case ("SQUARE") - Phi_X_vec = 0.01 + Phi_X_vec = 1.d-2 Call Set_Default_hopping_parameters_square(Hopping_Matrix_tmp,Ham_T_vec, Ham_Chem_vec, Phi_X_vec, Phi_Y_vec, & & Bulk, N_Phi_vec, N_FL, & & List, Invlist, Latt, Latt_unit ) @@ -317,7 +319,7 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, Case ("N_LEG_LADDER") Ham_T_vec = 1.d0 Ham_Tperp_vec = 1.d0 - Phi_X_vec = 0.01 + Phi_X_vec = 1.d-2 Call Set_Default_hopping_parameters_n_leg_ladder(Hopping_Matrix_tmp, Ham_T_vec, Ham_Tperp_vec, Ham_Chem_vec, Phi_X_vec, & & Phi_Y_vec, Bulk, N_Phi_vec, N_FL, & & List, Invlist, Latt, Latt_unit ) @@ -372,7 +374,7 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, Write(6,*) Op_tmp(1,1)%E(I) enddo Open(newunit=den_file, file="Den_H0", status="unknown") - delta = 2.d0 * Acos(-1.d0)/(Iscalar(Latt%L1_p,Latt%BZ1_p)) + delta = twopi/(Iscalar(Latt%L1_p,Latt%BZ1_p)) dom =(Op_tmp(1,1)%E(ndim) - Op_tmp(1,1)%E(1)) / dble(Nom) om = Op_tmp(1,1)%E(1) do nw = 1,Nom + 1 @@ -380,7 +382,7 @@ Subroutine Predefined_TrialWaveFunction(Lattice_type, Ndim, List,Invlist,Latt, Do I = 1,NDim Z = Z + 1.d0/cmplx(om - Op_tmp(1,1)%E(I),delta, kind=kind(0.d0)) enddo - Write(den_file,"(F14.7,2x,F14.7)") om, - Aimag(Z)/(dble(Ndim)*acos(-1.d0)) + Write(den_file,"(F14.7,2x,F14.7)") om, - Aimag(Z)/(dble(Ndim)*pi) om = om + dom enddo close(den_file) diff --git a/Prog/Set_random_mod.F90 b/Prog/Set_random_mod.F90 index cd00ff52d..8b37d4f8d 100644 --- a/Prog/Set_random_mod.F90 +++ b/Prog/Set_random_mod.F90 @@ -59,11 +59,12 @@ Subroutine Set_Random_number_Generator(File_seeds,Seed_in) Character (LEN=64), Intent(IN) :: File_seeds Integer, Intent(out) :: SEED_IN - Integer :: I, IERR + Integer :: IERR Integer, allocatable :: SEED_VEC(:) #ifdef MPI - INTEGER :: STATUS(MPI_STATUS_SIZE), irank_g, isize_g, igroup, ISIZE, IRANK + Integer :: I + INTEGER :: STATUS(MPI_STATUS_SIZE), ISIZE, IRANK CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) #endif diff --git a/Prog/WaveFunction_mod.F90 b/Prog/WaveFunction_mod.F90 index 501853dc5..e738770a8 100644 --- a/Prog/WaveFunction_mod.F90 +++ b/Prog/WaveFunction_mod.F90 @@ -80,7 +80,7 @@ Subroutine WF_overlap(WF_L, WF_R, Z_norm) ! Local - Integer :: N_Part, Ndim, n,ne + Integer :: N_Part, Ndim Complex (Kind=Kind(0.d0)), allocatable :: mat(:,:) Complex (Kind=Kind(0.d0)) :: alpha, beta diff --git a/Prog/Wrapgr_mod.F90 b/Prog/Wrapgr_mod.F90 index 5466219a1..d5e18e7ab 100644 --- a/Prog/Wrapgr_mod.F90 +++ b/Prog/Wrapgr_mod.F90 @@ -354,7 +354,7 @@ Subroutine Wrapgr_Random_update( GR,m,ntau, PHASE, N_Global_tau ) ! Space for local variables Integer :: n, Flip_length, nf, nf_eff, N_Type, ng_c, Flip_count - Real (Kind=Kind(0.d0)) :: T0_Proposal_ratio, T0_proposal,S0_ratio + Real (Kind=Kind(0.d0)) :: T0_Proposal_ratio, S0_ratio COMPLEX (Kind=Kind(0.d0)) :: Prev_Ratiotot, HS_Field, HS_New Logical :: Acc Character (Len=64) :: Mode @@ -444,7 +444,7 @@ subroutine Wrapgr_sort(Flip_length,Flip_list,Flip_value) ! Local integer :: swaps ! number of swaps made in one pass integer :: nc ! loop variable - integer :: temp, n ! temporary holder for making swap + integer :: temp ! temporary holder for making swap Complex (Kind=Kind(0.d0)) :: X diff --git a/Prog/cgr1_mod.F90 b/Prog/cgr1_mod.F90 index 2a16760a9..c811934cc 100644 --- a/Prog/cgr1_mod.F90 +++ b/Prog/cgr1_mod.F90 @@ -190,7 +190,10 @@ SUBROUTINE CGR(PHASE,NVAR, GRUP, udvr, udvl) COMPLEX (Kind=Kind(0.d0)), Dimension(:,:), Allocatable :: TPUP, RHS, TPUP_temp COMPLEX (Kind=Kind(0.d0)), Dimension(:) , Allocatable :: DUP INTEGER, Dimension(:), Allocatable :: IPVT, VISITED - COMPLEX (Kind=Kind(0.d0)) :: alpha, beta, Z, DLJ + COMPLEX (Kind=Kind(0.d0)) :: alpha, beta, Z +#if (defined(STAB3) || defined(STABLOG)) + COMPLEX (Kind=Kind(0.d0)) :: DLJ +#endif Integer :: I, J, N_size, info, LWORK, next, L Real (Kind=Kind(0.D0)) :: X, Xmax, sv diff --git a/Prog/main.F90 b/Prog/main.F90 index 3de7644e5..761d8a97f 100644 --- a/Prog/main.F90 +++ b/Prog/main.F90 @@ -155,10 +155,13 @@ Program Main Integer :: NBC, NSW Integer :: NTAU, NTAU1 - Character (len=64) :: file_seeds, file_dat, file_info, file_git_info + Character (len=64) :: file_seeds, file_info, file_git_info Integer :: Seed_in Complex (Kind=Kind(0.d0)) , allocatable, dimension(:,:) :: Initial_field +#ifdef HDF5 + Character (len=64) :: file_dat +#endif #ifdef HDF5 INTEGER(HID_T) :: file_id @@ -167,13 +170,17 @@ Program Main !General Integer :: NSTM, NT, NT1, NVAR - Integer :: Ierr, I,nf, nf_eff, nst, n, n1, N_op, NBin_eff + Integer :: I, nf, nf_eff, nst, n, n1, N_op, NBin_eff Integer :: tmp_Nt_sequential_start, tmp_Nt_sequential_end, tmp_N_Global_tau Logical :: Toggle, Toggle1 Complex (Kind=Kind(0.d0)) :: Phase, Z, Z1 Real (Kind=Kind(0.d0)) :: ZERO = 10D-8 Real (Kind=Kind(0.d0)) :: Mc_step_weight +#if defined(MPI) || defined(HDF5) + Integer :: Ierr +#endif + ! Storage for stabilization steps Integer, dimension(:), allocatable :: Stab_nt @@ -185,11 +192,10 @@ Program Main integer (kind=kind(0.d0)) :: count_bin_start, count_bin_end ! For MPI shared memory +#ifdef MPI character(64), parameter :: name="ALF_SHM_CHUNK_SIZE_GB" character(64) :: chunk_size_str Real (Kind=Kind(0.d0)) :: chunk_size_gb - -#ifdef MPI Integer :: Isize, Irank, Irank_g, Isize_g, color, key, igroup CALL MPI_INIT(ierr) diff --git a/Prog/observables_mod.F90 b/Prog/observables_mod.F90 index 7c693fe5c..879015186 100644 --- a/Prog/observables_mod.F90 +++ b/Prog/observables_mod.F90 @@ -143,11 +143,11 @@ Subroutine Obser_Latt_local_make(Obs, Nt, Filename, Latt, Latt_unit, Channel, dt !> \verbatim !> Name of file in which the bins will be written out. !> \endverbatim -!> @param [IN] Latt, Type(Lattice) +!> @param [INOUT] Latt, Type(Lattice) !> \verbatim !> Bravais lattice. Only gets linked, needs attribute target or pointer. !> \endverbatim -!> @param [IN] Latt_unit, Type(Unit_cell) +!> @param [INOUT] Latt_unit, Type(Unit_cell) !> \verbatim !> Unit cell. Only gets linked, needs attribute target or pointer. !> \endverbatim @@ -160,8 +160,12 @@ Subroutine Obser_Latt_local_make(Obs, Nt, Filename, Latt, Latt_unit, Channel, dt type(Obser_Latt_local), Intent(INOUT) :: Obs Integer, Intent(IN) :: Nt Character(len=64), Intent(IN) :: Filename - Type(Lattice), Intent(IN), target :: Latt - Type(Unit_cell), Intent(IN), target :: Latt_unit + ! INTENT(INOUT) required: flang rejects INTENT(IN) when a pointer to the dummy + ! argument is stored (Obs%Latt => Latt). Latt is not modified in this routine. + Type(Lattice), Intent(INOUT), target :: Latt + ! INTENT(INOUT) required: flang rejects INTENT(IN) when a pointer to the dummy + ! argument is stored (Obs%Latt_unit => Latt_unit). Latt_unit is not modified in this routine. + Type(Unit_cell), Intent(INOUT), target :: Latt_unit Character(len=*), Intent(IN) :: Channel Real(Kind=Kind(0.d0)), Intent(IN) :: dtau If (Nt > 1) then @@ -212,11 +216,11 @@ Subroutine Obser_Latt_make(Obs, Nt, Filename, Latt, Latt_unit, Channel, dtau) !> \verbatim !> Name of file in which the bins will be written out. !> \endverbatim -!> @param [IN] Latt, Type(Lattice) +!> @param [INOUT] Latt, Type(Lattice) !> \verbatim !> Bravais lattice. Only gets linked, needs attribute target or pointer. !> \endverbatim -!> @param [IN] Latt_unit, Type(Unit_cell) +!> @param [INOUT] Latt_unit, Type(Unit_cell) !> \verbatim !> Unit cell. Only gets linked, needs attribute target or pointer. !> \endverbatim @@ -233,8 +237,12 @@ Subroutine Obser_Latt_make(Obs, Nt, Filename, Latt, Latt_unit, Channel, dtau) type(Obser_Latt), Intent(INOUT) :: Obs Integer, Intent(IN) :: Nt Character(len=64), Intent(IN) :: Filename - Type(Lattice), Intent(IN), target :: Latt - Type(Unit_cell), Intent(IN), target :: Latt_unit + ! INTENT(INOUT) required: flang rejects INTENT(IN) when a pointer to the dummy + ! argument is stored (Obs%Latt => Latt). Latt is not modified in this routine. + Type(Lattice), Intent(INOUT), target :: Latt + ! INTENT(INOUT) required: flang rejects INTENT(IN) when a pointer to the dummy + ! argument is stored (Obs%Latt_unit => Latt_unit). Latt_unit is not modified in this routine. + Type(Unit_cell), Intent(INOUT), target :: Latt_unit Character(len=*), Intent(IN) :: Channel Real(Kind=Kind(0.d0)), Intent(IN) :: dtau Allocate (Obs%Obs_Latt(Latt%N, Nt, Latt_unit%Norb, Latt_unit%Norb)) @@ -361,10 +369,12 @@ Subroutine Print_bin_Latt(Obs, Group_Comm) ! Local Integer :: Ns, Nt, no, no1, I, Ntau Complex (Kind=Kind(0.d0)), allocatable, target :: Tmp(:,:,:,:) + Character (len=64) :: File_pr, File_suff +#if defined OBS_LEGACY Real (Kind=Kind(0.d0)) :: x_p(2) - Complex (Kind=Kind(0.d0)) :: Sign_bin - Character (len=64) :: File_pr, File_suff, File_aux, tmp_str + Character (len=64) :: File_aux, tmp_str logical :: File_exists +#endif #ifdef HDF5 Character (len=7), parameter :: File_h5 = "data.h5" Character (len=64) :: filename, groupname, obs_dsetname, bak_dsetname, sgn_dsetname @@ -377,7 +387,6 @@ Subroutine Print_bin_Latt(Obs, Group_Comm) #endif #ifdef MPI Complex (Kind=Kind(0.D0)), allocatable :: Tmp1(:) - Complex (Kind=Kind(0.d0)) :: Z Real (Kind=Kind(0.d0)) :: X Integer :: Ierr, Isize, Irank INTEGER :: irank_g, isize_g, igroup @@ -582,11 +591,14 @@ Subroutine Print_bin_Latt_Local(Obs, Group_Comm) Integer, Intent(In) :: Group_Comm ! Local - Integer :: Ns, Nt, no, no1, I, Ntau + Integer :: Ns, I, Ntau + Character (len=64) :: File_pr, File_suff +#if defined OBS_LEGACY + Integer :: Nt, no Real (Kind=Kind(0.d0)) :: x_r(2) - Complex (Kind=Kind(0.d0)) :: Sign_bin - Character (len=64) :: File_pr, File_suff, File_aux, tmp_str + Character (len=64) :: File_aux, tmp_str logical :: File_exists +#endif #ifdef HDF5 Character (len=7), parameter :: File_h5 = "data.h5" Character (len=64) :: filename, groupname, obs_dsetname, sgn_dsetname @@ -599,7 +611,6 @@ Subroutine Print_bin_Latt_Local(Obs, Group_Comm) #endif #ifdef MPI Complex (Kind=Kind(0.D0)), allocatable :: Tmp(:,:,:) - Complex (Kind=Kind(0.d0)) :: Z Real (Kind=Kind(0.d0)) :: X Integer :: Ierr, Isize, Irank INTEGER :: irank_g, isize_g, igroup @@ -773,8 +784,11 @@ Subroutine Print_bin_Vec(Obs,Group_Comm) ! Local Integer :: I - Character (len=64) :: File_pr, File_suff, File_aux + Character (len=64) :: File_pr +#if defined OBS_LEGACY + Character (len=64) :: File_aux logical :: File_exists +#endif #if defined HDF5 Character (len=7), parameter :: File_h5 = "data.h5" diff --git a/Prog/parse_ham_mod.py b/Prog/parse_ham_mod.py index c8c8988b6..c32f2d5c3 100644 --- a/Prog/parse_ham_mod.py +++ b/Prog/parse_ham_mod.py @@ -172,7 +172,6 @@ def create_read_write_par(filename, parameters, ham_name): #ifdef MPI Integer :: Isize, Irank, igroup, irank_g, isize_g - Integer :: STATUS(MPI_STATUS_SIZE) 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) diff --git a/Prog/tau_m_mod.F90 b/Prog/tau_m_mod.F90 index 344bcc46b..52db95442 100644 --- a/Prog/tau_m_mod.F90 +++ b/Prog/tau_m_mod.F90 @@ -71,8 +71,8 @@ SUBROUTINE TAU_M(udvst, GR, PHASE, NSTM, NWRAP, STAB_NT, LOBS_ST, LOBS_EN) Complex (Kind=Kind(0.d0)), Dimension(:,:), Allocatable :: HLP4, HLP5, HLP6 Complex (Kind=Kind(0.d0)) :: Z - Integer :: I, J, nf, nf_eff, NT, NT1, NTST, NST, N, N_type - Real (Kind=Kind(0.d0)) :: spin, Mc_step_Weight + Integer :: I, J, nf, nf_eff, NT, NT1, NTST, NST + Real (Kind=Kind(0.d0)) :: Mc_step_Weight Allocate( HLP4(Ndim,Ndim), HLP5(Ndim,Ndim), HLP6(Ndim,Ndim) ) Allocate( G00(Ndim,Ndim,N_FL), G0T(Ndim,Ndim,N_FL), GT0(Ndim,Ndim,N_FL), GTT(Ndim,Ndim,N_FL) ) diff --git a/Prog/tau_p_mod.F90 b/Prog/tau_p_mod.F90 index 0cd4e2244..92d8d8233 100644 --- a/Prog/tau_p_mod.F90 +++ b/Prog/tau_p_mod.F90 @@ -87,17 +87,13 @@ SUBROUTINE Tau_p(udvl, udvr, udvst, GR, PHASE, NSTM, STAB_NT, NST_IN, LOBS_ST, L ! Local. CLASS(UDV_State), Dimension(:), ALLOCATABLE :: udvr_local - Complex (Kind=Kind(0.d0)) :: DETZ, ZK, DET1(2) + Complex (Kind=Kind(0.d0)) :: DETZ Complex (Kind=Kind(0.d0)), Dimension(:,:,:), Allocatable :: GRUPB, GRUP Complex (Kind=Kind(0.d0)), Dimension(:,:,:), Allocatable :: G00UP, G0TUP, GT0UP, GTTUP Complex (Kind=Kind(0.d0)), Dimension(:,:,:), Allocatable :: G00UP_T, G0TUP_T, GT0UP_T, GTTUP_T - Complex (Kind=Kind(0.d0)), allocatable :: TEMP(:,:), TMPUP(:,:) + Complex (Kind=Kind(0.d0)), allocatable :: TEMP(:,:) - Real (Kind=kind(0.d0)) :: XMEAN_DYN, XMAX_DYN - - Integer :: NTAUIN, NTDM, LFAM, NFAM, N_Part, LQ , I, NCON, NF, nf_eff, NFLAG, NL, NT1, NT_ST, NT, NTAU, NTAU1, n, NCHECK - - Real (Kind=Kind(0.d0)) :: XMEAN, XMAX + Integer :: LQ , I, NF, nf_eff, NT1, NT_ST, NT, NTAU, NTAU1, NCHECK Real (Kind=Kind(0.d0)) :: Mc_step_weight LQ = ndim diff --git a/Prog/udv_state_mod.F90 b/Prog/udv_state_mod.F90 index 71c2e8d6c..7f45184dd 100644 --- a/Prog/udv_state_mod.F90 +++ b/Prog/udv_state_mod.F90 @@ -256,8 +256,6 @@ SUBROUTINE testscale_UDV_state(this) CLASS(UDV_State), INTENT(IN) :: this #if !defined(STABLOG) - real (Kind=Kind(this%D(1))) :: dummy_dp - ! Check if any scale is NaN if ( any(this%D /= this%D) ) then write(error_unit,*) @@ -268,16 +266,16 @@ SUBROUTINE testscale_UDV_state(this) ! ATTENTION, the test assumes a (mostly) sorted array D [real and positive numbers]! ! Check if largest scale is approaching the largest representable value - if ( dble(this%D(1)) > 0.1*huge(dummy_dp) .and. trigger_scale_warning) then + if ( dble(this%D(1)) > 0.1*huge(real(0.0, kind=kind(this%D(1)))) .and. trigger_scale_warning) then write(error_unit,*) write(error_unit,*) "Warning: Largest scale is approaching the largest representable value." write(error_unit,*) " Consider switching to LOG." trigger_scale_warning = .false. end if ! Check if myVariable is approaching the smallest representable value - if ( dble(this%D(this%n_part)) < 10.0*tiny(dummy_dp) .and. trigger_scale_warning) then + if ( dble(this%D(this%n_part)) < 10.0*tiny(real(0.0, kind=kind(this%D(1)))) .and. trigger_scale_warning) then write(error_unit,*) - write(error_unit,*) "Warning: Smallest scale is approaching the smalles representable value." + write(error_unit,*) "Warning: Smallest scale is approaching the smallest representable value." write(error_unit,*) " Consider switching to LOG." trigger_scale_warning = .false. end if diff --git a/Prog/upgrade_mod.F90 b/Prog/upgrade_mod.F90 index dad9d77f4..11dbd3299 100644 --- a/Prog/upgrade_mod.F90 +++ b/Prog/upgrade_mod.F90 @@ -129,7 +129,7 @@ Subroutine Upgrade2(GR,N_op,NT,PHASE,Hs_new, Prev_Ratiotot, S0_ratio, T0_proposa Integer :: n,m,nf, nf_eff, i, Op_dim, op_dim_nf Complex (Kind=Kind(0.d0)) :: Z, D_Mat, myexp, s1, s2 - Real (Kind=Kind(0.d0)) :: Weight, tmp_r + Real (Kind=Kind(0.d0)) :: Weight Complex (Kind=Kind(0.d0)) :: alpha, beta, g_loc Complex (Kind=Kind(0.d0)), Dimension(:, :), Allocatable :: Mat, Delta Complex (Kind=Kind(0.d0)), Dimension(:, :), Allocatable :: u, v diff --git a/configure.sh b/configure.sh index d501a1f76..3953ce79b 100755 --- a/configure.sh +++ b/configure.sh @@ -207,7 +207,10 @@ GNUOPTFLAGS="-cpp -O3 -ffree-line-length-none -ffast-math" # uncomment the next line if you want to use additional openmp parallelization GNUOPTFLAGS="${GNUOPTFLAGS} -fopenmp" # GNUDEVFLAGS="-Wconversion -Werror -fcheck=all -ffpe-trap=invalid,zero,overflow,underflow,denormal" -GNUDEVFLAGS="-Wconversion -Werror -Wno-error=cpp -fcheck=all -g -fbacktrace -fmax-errors=10" +GNUDEVFLAGS="-Wconversion -fcheck=all -g -fbacktrace -fmax-errors=10" +GNUDEVFLAGS="${GNUDEVFLAGS} -pedantic" +# GNUDEVFLAGS="${GNUDEVFLAGS} -Wall -Wno-error=unused-function -Wno-error=unused-variable -Wno-error=unused-dummy-argument -Wno-error=maybe-uninitialized" +GNUDEVFLAGS="${GNUDEVFLAGS} -Werror -Wno-error=cpp" GNUUSEFULFLAGS="-std=f2008" # default optimization flags for PGI compiler @@ -533,6 +536,7 @@ if [ -n "${ALF_FLAGS_EXT+x}" ]; then fi ALF_FLAGS_QRREF="${F90OPTFLAGS} ${ALF_FLAGS_EXT}" +ALF_FLAGS_QRREF="$(echo "$ALF_FLAGS_QRREF" | sed 's| -pedantic||')" # Modules need to know the programm configuration since entanglement needs MPI ALF_FLAGS_MODULES="${F90OPTFLAGS} ${PROGRAMMCONFIGURATION} ${ALF_FLAGS_EXT}" ALF_FLAGS_ANA="${F90USEFULFLAGS} ${F90OPTFLAGS} ${ALF_INC} ${ALF_FLAGS_EXT}" diff --git a/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90 b/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90 index 8437e33a0..6834fd7e5 100644 --- a/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90 +++ b/testsuite/Analysis.tests/1-maxent-kernel-overflow.F90 @@ -18,14 +18,12 @@ Program Test_MaxEnt_Kernel_Overflow Use MaxEnt_Wrapper_mod implicit none - real(Kind=Kind(0.d0)) :: beta, tau, om, pi + real(Kind=Kind(0.d0)) :: beta, tau, om real(Kind=Kind(0.d0)) :: val, ref, relerr real(Kind=Kind(0.d0)) :: tol real(Kind=Kind(0.d0)) :: taus(4) integer :: it, nerr - ! Use the same single-precision pi constant as maxent_wrapper_mod - pi = 3.1415927 beta = 10.d0 tol = 1.d-12 nerr = 0