diff --git a/Documentation/Doc.idx b/Documentation/Doc.idx index 0b177450a5b293aed654e9dabb91445a291f3a0e..3c1d526e8b6d36649b687f58da8ca432f1d91518 100644 --- a/Documentation/Doc.idx +++ b/Documentation/Doc.idx @@ -5,8 +5,8 @@ \indexentry{\texttt{Precision Green} }{27} \indexentry{\texttt{Precision Phase}}{27} \indexentry{\path{NWrap}}{27} -\indexentry{\texttt{Square} }{67} -\indexentry{\path{Bilayer_square}}{68} -\indexentry{\path{N_leg_ladder}}{68} -\indexentry{\path{Honeycomb}}{68} -\indexentry{\texttt{Bilayer\_honeycomb}}{68} +\indexentry{\texttt{Square} }{69} +\indexentry{\path{Bilayer_square}}{69} +\indexentry{\path{N_leg_ladder}}{70} +\indexentry{\path{Honeycomb}}{70} +\indexentry{\texttt{Bilayer\_honeycomb}}{70} diff --git a/Documentation/Hubbard_Plain.tex b/Documentation/Hubbard_Plain.tex index abca9a1792d72e675bce7598b2550b8b594d6234..830e40ff283d1ca57228c43e7f80934b5ae35e99 100644 --- a/Documentation/Hubbard_Plain.tex +++ b/Documentation/Hubbard_Plain.tex @@ -83,8 +83,9 @@ We can test the correct formatting of the parameters by calling: The main program will call the subroutine \texttt{Ham\_set} in the submodule \path{Hamiltonian_Hubbard_Plain_Vanilla_smod.F90} which specify the model. The routine \texttt{Ham\_set} will first read the parameter file \texttt{parameters} (see Sec.~\ref{sec:input}) \texttt{Call read\_parameters}; then set the lattice: \texttt{Call Ham\_latt}; set the hopping: \texttt{Call Ham\_hop}; - set the interaction: \texttt{call Ham\_V}; and if required, set the trial wave function: \texttt{call Ham\_trial}. - + set the interaction: \texttt{call Ham\_V}; and if required, set the trial wave function: \texttt{call Ham\_trial}. In the subroutine \texttt{Ham\_set} one + will equally have to specify if a symmetry relates different flavors. This functionality is described in Sec.~\ref{sec:flavor_sym} and one enables it by allocating tine array \texttt{Calc\_Fl}. + %------------------------------------------------------------ \subsection{The lattice: \texttt{Ham\_latt}} \label{U_PV_Ham_latt} %------------------------------------------------------------ @@ -167,16 +168,18 @@ Do nf = 1,N_FL if (nf == 2) X = -1.d0 Do i = 1,Ndim nc = nc + 1 - Op_V(i,nf)%P(1) = I + Op_V(i,nf)%P(1) = i Op_V(i,nf)%O(1,1) = cmplx(1.d0, 0.d0, kind(0.D0)) Op_V(i,nf)%g = X*SQRT(CMPLX(DTAU*ham_U/2.d0, 0.D0, kind(0.D0))) - Op_V(i,nf)%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) + Op_V(i,nf)%alpha = cmplx(-0.5d0, 0.d0, kind(0.D0)) Op_V(i,nf)%type = 2 Call Op_set( Op_V(i,nf) ) Enddo Enddo \end{lstlisting} -The code above makes it explicit that there is a sign difference between the coupling of the HS field in the two flavor sectors. +The code above makes it explicit that there is a sign difference between the coupling of the HS field in the two flavor sectors. Hence, +\texttt{Op\_V(i,nf)} encodes $ e^{ X \sqrt{\Delta \tau U/2} \left( \hat{c}^{\dagger}_{\ve{i},{\texttt{nf}}} \hat{c}^{\phantom\dagger}_{\ve{i},{\texttt{nf}}} + \alpha \right) } $ with $X=1 $ for $\texttt{nf}=1$ and $X=-1 $ for $\texttt{nf}=2$. Strictly speaking $X$ can be omitted. However, it is required when using the flavor symmetry option in the presence of particle-hole symmetry (see Sec.~\ref{sec:flavor_sym}). + %------------------------------------------------------------ \subsection{The trial wave function: \texttt{Ham\_Trial}} \label{U_PV_Ham_Trial} @@ -349,6 +352,53 @@ are then given by: The handling of time-displaced correlation functions is identical to that of equal-time correlations. + +%------------------------------------------------------------------------------------- +\subsection[Flavor symmetries]{ Flavor symmetries. \texttt{ weight\_reconstruction(weight), GR\_reconstruction(GR)}, and \texttt{GRT\_reconstruction(GT0, G0T)}} \label{sec:flavor_symm_vanilla_hubbard} +%------------------------------------------------------------------------------------- +At zero chemical potential, and for repulsive interactions, the plain-vanilla Hubbard model enjoys a partial particle-hole symmetry that +for each HS field configuration, maps one flavor one flavor (i.e. spin sector) onto the other: +\begin{equation} + \hat{P} z \hat{c}^{\dagger}_{\ve{i},\downarrow} \hat{P}^{-1} = z^{*} e^{i \, \ve{i} \cdot \ve{Q} }\hat{c}^{\phantom\dagger}_{\ve{i},\uparrow}. +\end{equation} +Here $ \ve{Q}$ is the antiferromagnetic wave vector. +Note that in the presence of an orbital magnetic field, or of twisted boundary conditions that couples symmetrically to the flavor (spin) +degree of freedom, +the anti-unitarity of the transformation is required. A consequence of this symmetry, and for a given HS field configuration, the following holds for the equal time Green function. + +\begin{align} +\begin{aligned} + \texttt{G00(x,y,}\uparrow\texttt{)} &= \langle \langle \hat{c}^{\phantom\dagger}_{\ve{x},\uparrow} (0) \hat{c}^{\dagger}_{\ve{y},\uparrow} (0) \rangle \rangle + \;=\; e^{i (\ve{y} - \ve{x}) \cdot \ve{Q}} \, \, \overline{ \langle \langle \hat{c}^{\dagger}_{\ve{x},\downarrow} (0) + \hat{c}^{\phantom\dagger}_{\ve{y},\uparrow} (0) \rangle } + \\ & = \delta_{\ve{x},\ve{y}} - e^{i (\ve{y} - \ve{x}) \cdot \ve{Q}}\, \, \overline{ \texttt{G00(y,x,}\downarrow\texttt{)} } +\end{aligned} +\end{align} + + +For the attractive Hubbard model $U<0$, the up and down sectors are related by time reversal symmetry. +\begin{equation} +\hat{T} z +\begin{pmatrix} + \hat{c}^{\phantom\dagger}_{\ve{i},\uparrow} \\ + \hat{c}^{\phantom\dagger}_{\ve{i},\downarrow} +\end{pmatrix} + \hat{T }^{-1} = z^{*} i \sigma_y +\begin{pmatrix} + \hat{c}^{\phantom\dagger}_{\ve{i},\uparrow} \\ + \hat{c}^{\phantom\dagger}_{\ve{i},\downarrow} +\end{pmatrix} +\end{equation} +Of course, we have assumed that the hopping remains invariant under time reversal. As a consequence of this symmetry: +\begin{align} +\begin{aligned} + \texttt{G00(x,y,}\uparrow\texttt{)} = \overline{ \texttt{G00(x,y}\downarrow\texttt{)} } +\end{aligned} +\end{align} + +The usage of the flavor symmetry is described in Sec.~\ref{sec:flavor_sym}. Only one flavor has to be computed and the routines +\texttt{GR\_reconstruction(GR)}, and \texttt{GRT\_reconstruction(GT0, G0T)} will reconstruct the equal time and time displaced Green functions respectively for one flavor given the other. Hence we gain a factor two in computing time. We note that since both symmetries are anti-unitary, the weights between the two sectors are related by a complex conjugation. This is specified in the routine \texttt{ weight\_reconstruction(weight)}. + %------------------------------------------------------------------------------------- \subsection{Numerical precision}\label{sec:prec_spin} %------------------------------------------------------------------------------------- diff --git a/Documentation/doc.pdf b/Documentation/doc.pdf index 040e54a033867ad4a383680e108b78bdd307d5c6..253e53c04d9720b7451059267cad6f9c676ecbcd 100644 Binary files a/Documentation/doc.pdf and b/Documentation/doc.pdf differ diff --git a/Documentation/implementation.tex b/Documentation/implementation.tex index 1e6664d4490fdd84957ce69ae7949d60bf33ceae..945ad2e7646dd57741a3df1f3dc556d5523fc2ef 100644 --- a/Documentation/implementation.tex +++ b/Documentation/implementation.tex @@ -472,6 +472,7 @@ To simplify the implementation of a new Hamiltonian, ALF comes with a set of pr \hl{\texttt{Projector}} & \path{logical} & Enable projector code \\ \texttt{Group\_Comm} & \path{int} & Group communicator for MPI \\ \hl{\texttt{Symm}} & \path{logical} & Symmetric trotter \\ + \texttt{Calc\_Fl} & \path{logical} & Explicitly calculate flavors (optional) \\ \bottomrule \toprule Private Variable & Type & Description \\ @@ -549,7 +550,13 @@ To simplify the implementation of a new Hamiltonian, ALF comes with a set of pr \texttt{Pr\_obs} & Writes the observables to disk by calling \texttt{Print\_bin} of the \texttt{Observables} module. - Usually, this doesn't have to be modified. & + Usually, this doesn't have to be modified. & \\ + \texttt{weight\_reconstruction} & + Reconstruction of the weight for `inactive' flavors. & \ref{sec:flavor_sym} \\ + \texttt{GR\_reconstruction} & + Reconstruction of the Greens function GR for the `inactive' flavors. & \ref{sec:flavor_sym} \\ + \texttt{GRT\_reconstruction} & + Reconstruction of the time-displaced Greens functions G0T and GT0. & \ref{sec:flavor_sym} \\\bottomrule \end{tabular} \caption{Typebound procedures bound to type \texttt{ham\_base}. To define a new model, at least \texttt{Ham\_Set} has to be overloaded in the Hamiltonian submodule. For measurements \texttt{Alloc\_obs}, \texttt{Obser} (and \texttt{ObserT} for time displaced observables) are necessary. @@ -614,4 +621,26 @@ we define the corresponding doubled array \texttt{Op\_V(M$_V$+M$_I$,N$_\mathrm \end{itemize} +\subsubsection{Flavor symmetries} +\label{sec:flavor_sym} +This code allows the use of time-reversal or particle-hole symmetry to accelerate the algorithm by only explicitly calculating a subset of flavors and reconstructing the complement by symmetry. +Here, a pair of flavors $(n_\mathrm{fl},\bar{n}_\mathrm{fl})$ is related by a unitary or anti-unitary symmetry, including particle-hole transformations, such that +\begin{equation} + \mathcal{U}^{-1} \texttt{Op\_V(i,n$_\mathrm{fl}$)} \mathcal{U} = \texttt{Op\_V(i,$\bar{n}_\mathrm{fl}$)}; \quad + \mathcal{U}^{-1} \texttt{Op\_T(i,n$_\mathrm{fl}$)} \mathcal{U} = \texttt{Op\_T(i,$\bar{n}_\mathrm{fl}$)} +\end{equation} +for any given $i$. Note that \texttt{Op\_V(i,n$_\mathrm{fl}$)} includes a constant shift $\alpha$ to absorb contributions from the commutator in case of particle-hole symmetries. +For example, the particle-hole symmetry requires a non-zero shift $\alpha=1/2$ in the $M_z$ decoupling of the Hubbard interaction to map the upspin to the downspin. +This acceleration is activated by allocating \texttt{Calc\_Fl} in \texttt{Ham\_set} and setting the `active' flavors $n_\mathrm{fl}$ to \texttt{.True.} and the symmetry-related flavors $\bar{n}_\mathrm{fl}$ to \texttt{.False.}. + +This symmetry allows to reconstruct one flavor, say $\bar{n}_\mathrm{fl}$, from the other, $n_\mathrm{fl}$. For unitary symmetries, the weight is given by $Z(\bar{n}_\mathrm{fl})=Z(n_\mathrm{fl})$, while anti-unitary symmetries lead to $Z(\bar{n}_\mathrm{fl})=Z(n_\mathrm{fl})^*$. This relation has to be provided by the user in \texttt{weight\_reconstruction}. Upto entry, all weights $Z(n_\mathrm{fl})$ have be explicilty calculated by ALF and the user has to fill all `inactive' flavors $\bar{n}_\mathrm{fl}$ of the array $Z$. +Similarly, the subroutines \texttt{GR\_reconstruction} and \texttt{GRT\_reconstruction} have to be override to provide the symmetry reconstruction of the `inactive' flavors $\bar{n}_\mathrm{fl}$ of the equal-time and time-displaced Greens-function out of $n_\mathrm{fl}$, respectively. + +Finally we note that if the projective algorithm is used then the trial wave function also has to satisfy the aforementioned symmetry. +In particular, +assume that the trial wave function corresponds to the ground state of a single particle Hamiltonian, $\hat{H}_T(n_\mathrm{fl}) $, then we +will require that +\begin{equation} +\mathcal{U}^{-1} \hat{H}_T(n_\mathrm{fl}) \mathcal{U} = \hat{H}_T(\bar{n}_\mathrm{fl}) +\end{equation} \ No newline at end of file diff --git a/Prog/Global_mod.F90 b/Prog/Global_mod.F90 index fe651e2af258f69e936bd3903e1049243abba2b7..068c637021b5880abd19c9badc164d3ea5ab8934 100644 --- a/Prog/Global_mod.F90 +++ b/Prog/Global_mod.F90 @@ -121,13 +121,13 @@ Module Global_mod !> Local variables. - Integer :: NST, NSTM, NF, 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, I_Partner, n_step, N_count, N_part Type (Fields) :: nsigma_old Real (Kind=Kind(0.d0)) :: T0_Proposal_ratio, Weight, Weight1 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)), 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) + Complex (Kind=Kind(0.d0)) :: Ratio(2), Ratio_p(2),Phase_array(N_FL) Logical :: TOGGLE, L_Test Integer, allocatable :: List_partner(:), List_masters(:) @@ -170,11 +170,15 @@ Module Global_mod ! Set old weight. storage = "Full" Call Compute_Fermion_Det(Phase_det_old,Det_Vec_old, udvl, udvst, Stab_nt, storage) - Phase_old = cmplx(1.d0,0.d0,kind=kind(0.d0)) - Do nf = 1,N_Fl - Phase_old = Phase_old*Phase_det_old(nf) + Phase_array=1.0d0 + Do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) + Phase_array(nf) = Phase_det_old(nf) + Call Op_phase(Phase_array(nf),OP_V,Nsigma,nf) Enddo - Call Op_phase(Phase_old,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase_old=product(Phase_array) + Phase_old=Phase_old**N_SUN endif !> Store old configuration nsigma_old%f = nsigma%f @@ -240,11 +244,15 @@ Module Global_mod storage = "Empty" Call Compute_Fermion_Det(Phase_det_new,Det_Vec_new, udvl, udvst, Stab_nt, storage) - Phase_new = cmplx(1.d0,0.d0,kind=kind(0.d0)) - Do nf = 1,N_Fl - Phase_new = Phase_new*Phase_det_new(nf) + Phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0)) + Do nf_eff = 1,N_Fl_eff + nf=Calc_FL_map(nf_eff) + Phase_array(nf) = Phase_det_new(nf) + Call Op_phase(Phase_array(nf),OP_V,Nsigma,nf) Enddo - Call Op_phase(Phase_new,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase_new=product(Phase_array) + Phase_new=Phase_new**N_SUN T0_Proposal_ratio = 1.d0 Ratiotot = Compute_Ratio_Global(Phase_Det_old, Phase_Det_new, & @@ -309,11 +317,12 @@ Module Global_mod if (Tempering_calc_det) then ! If move has been accepted, no use to recomute storage If (.not.TOGGLE) then - DO nf = 1,N_FL + DO nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvl(nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%reset('l') + CALL udvl(nf_eff)%reset('l') endif ENDDO DO NST = NSTM-1,1,-1 @@ -321,8 +330,8 @@ Module Global_mod NT = Stab_nt(NST ) !Write(6,*) NT1,NT, NST CALL WRAPUL(NT1,NT, udvl) - Do nf = 1,N_FL - udvst(NST, nf) = udvl(nf) + Do nf_eff = 1,N_FL_eff + udvst(NST, nf_eff) = udvl(nf_eff) ENDDO ENDDO NT1 = stab_nt(1) @@ -330,12 +339,16 @@ Module Global_mod Endif ! Compute the Green functions so as to provide correct starting point for the sequential updates. NVAR = 1 - Phase = cmplx(1.d0,0.d0,kind(0.d0)) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) - Phase = Phase*Z + Phase_array = cmplx(1.d0,0.d0,kind(0.d0)) + do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) + CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf_eff), udvl(nf_eff)) + call Op_phase(Z,OP_V,Nsigma,nf) + Phase_array(nf)=Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase=product(Phase_array) + Phase=Phase**N_SUN else ! Send >> Phase, GR, udvr, udvl, udvst << to new node ! First step: Each node sends to IRANK=0 its value nsigma_irank, @@ -367,15 +380,15 @@ Module Global_mod CALL MPI_Sendrecv_Replace(GR, n_GR, MPI_COMPLEX16, nsigma_old_irank, 0, & & nsigma_irank, 0, MPI_COMM_WORLD, STATUS, IERR) - do nf = 1,N_Fl - CALL udvr(nf)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) + do nf_eff = 1,N_Fl_eff + CALL udvr(nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) enddo - do nf = 1,N_Fl - CALL udvl(nf)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) + do nf_eff = 1,N_Fl_eff + CALL udvl(nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) enddo do NST = 1, NSTM - do nf = 1,N_Fl - CALL udvst(NST, nf)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) + do nf_eff = 1,N_Fl_eff + CALL udvst(NST, nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) enddo enddo endif @@ -449,14 +462,14 @@ Module Global_mod ! Local variables. - Integer :: NST, NSTM, NF, NT, NT1, NVAR,N, N1,N2, I, NC, N_part,j + Integer :: NST, NSTM, NF, NT, NT1, NVAR,N, N1,N2, I, NC, N_part,j, nf_eff Real (Kind=Kind(0.d0)) :: 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)), 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 - Complex (Kind=Kind(0.d0)) :: Ratio(2) + Complex (Kind=Kind(0.d0)) :: Ratio(2), Phase_array(N_Fl) Logical :: TOGGLE, L_Test Real (Kind=Kind(0.d0)) :: size_clust Real (Kind=Kind(0.d0)) :: ratio_2_test @@ -469,42 +482,53 @@ Module Global_mod NSTM = Size(udvst, 1) call nsigma_old%make(n1, n2) - Allocate ( Det_vec_old(NDIM,N_FL), Det_vec_new(NDIM,N_FL), Det_vec_test(NDIM,N_FL) ) - Allocate ( Phase_Det_new(N_FL), Phase_Det_old(N_FL) ) + Allocate ( Det_vec_old(NDIM,N_FL_eff), Det_vec_new(NDIM,N_FL_eff), Det_vec_test(NDIM,N_FL_eff) ) + Allocate ( Phase_Det_new(N_FL_eff), Phase_Det_old(N_FL_eff) ) L_test = .false. ! Set old weight. storage = "Full" Call Compute_Fermion_Det(Phase_det_old,Det_Vec_old, udvl, udvst, Stab_nt, storage) - Phase_old = cmplx(1.d0,0.d0,kind=kind(0.d0)) - Do nf = 1,N_Fl - Phase_old = Phase_old*Phase_det_old(nf) + Phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0)) + Do nf_eff = 1,N_Fl_eff + nf=Calc_FL_map(nf_eff) + Phase_array(nf) = Phase_det_old(nf_eff) + Call Op_phase(Phase_array(nf),OP_V,Nsigma,nf) Enddo - Call Op_phase(Phase_old,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase_old=product(Phase_array) + Phase_old=Phase_old**N_SUN If (L_test) then ! Testing - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvr(nf)%reset('r',WF_R(nf)%P) + CALL udvr(nf_eff)%reset('r',WF_R(nf)%P) else - CALL udvr(nf)%reset('r') + CALL udvr(nf_eff)%reset('r') endif Enddo NVAR = 1 - Phase = cmplx(1.d0,0.d0,kind(0.d0)) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) - Phase = Phase*Z + Phase_array = cmplx(1.d0,0.d0,kind(0.d0)) + do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) + CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf_eff), udvl(nf_eff)) + call Op_phase(Z,OP_V,Nsigma,nf) + Phase_array(nf)=Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) - Do Nf = 1,N_FL - Call DET_C_LU(GR(:,:,nf),Det_vec_test(:,nf),Ndim) - Z = Phase_det_old(nf) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase=product(Phase_array) + Phase=Phase**N_SUN + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) + Call DET_C_LU(GR(:,:,nf),Det_vec_test(:,nf_eff),Ndim) + !!!ATTENTION HOW DOES BELOW DEPEND ON NF BLOCK symm + Z = Phase_det_old(nf_eff) ratio_2_test=0.d0 DO I = 1,Ndim - Z = Z*Det_vec_test(I,nf)/ABS(Det_vec_test(I,nf)) - ratio_2_test=ratio_2_test+log(ABS(Det_vec_test(I,nf)))+Det_vec_old(I,nf) + Z = Z*Det_vec_test(I,nf_eff)/ABS(Det_vec_test(I,nf_eff)) + ratio_2_test=ratio_2_test+log(ABS(Det_vec_test(I,nf_eff)))+Det_vec_old(I,nf_eff) Enddo Z=Z*cmplx(exp(ratio_2_test),0.d0,kind(0.d0)) Write(6,*) 'Testing weight: ', Z @@ -527,11 +551,15 @@ Module Global_mod storage = "Empty" Call Compute_Fermion_Det(Phase_det_new,Det_Vec_new, udvl, udvst, Stab_nt, storage) - Phase_new = cmplx(1.d0,0.d0,kind=kind(0.d0)) - Do nf = 1,N_Fl - Phase_new = Phase_new*Phase_det_new(nf) + Phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0)) + Do nf_eff = 1,N_Fl_eff + nf=Calc_FL_map(nf_eff) + Phase_array(nf) = Phase_det_new(nf) + Call Op_phase(Phase_array(nf),OP_V,Nsigma,nf) Enddo - Call Op_phase(Phase_new,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase_new=product(Phase_array) + Phase_new=Phase_new**N_SUN Ratiotot = Compute_Ratio_Global(Phase_Det_old, Phase_Det_new, & & Det_vec_old, Det_vec_new, nsigma_old, T0_Proposal_ratio, Ratio) @@ -566,11 +594,12 @@ Module Global_mod If (NC > 0 ) then If (.not.TOGGLE) then - DO nf = 1,N_FL + DO nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvl(nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%reset('l') + CALL udvl(nf_eff)%reset('l') endif ENDDO DO NST = NSTM-1,1,-1 @@ -578,8 +607,8 @@ Module Global_mod NT = Stab_nt(NST ) !Write(6,*) NT1,NT, NST CALL WRAPUL(NT1,NT, udvl) - Do nf = 1,N_FL - udvst(NST, nf) = udvl(nf) + Do nf_eff = 1,N_FL_eff + udvst(NST, nf_eff) = udvl(nf_eff) ENDDO ENDDO NT1 = stab_nt(1) @@ -587,12 +616,16 @@ Module Global_mod Endif !Compute the Green functions so as to provide correct starting point for the sequential updates. NVAR = 1 - Phase = cmplx(1.d0,0.d0,kind(0.d0)) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) - Phase = Phase*Z + Phase_array = cmplx(1.d0,0.d0,kind(0.d0)) + do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) + CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf_eff), udvl(nf_eff)) + call Op_phase(Z,OP_V,Nsigma,nf) + Phase_array(nf)=Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase=product(Phase_array) + Phase=Phase**N_SUN endif @@ -645,29 +678,33 @@ Module Global_mod Complex (Kind=Kind(0.d0)), INTENT(out) :: Ratio(2) ! Local - Integer :: Nf, i, nt - Complex (Kind=Kind(0.d0)) :: Z, Z1 + Integer :: Nf, i, nt, nf_eff + Complex (Kind=Kind(0.d0)) :: Z, Z1, Ratio_1_array(N_FL), Ratio_2_array(N_FL) Real (Kind=Kind(0.d0)) :: X, Ratio_2 Ratio = cmplx(0.d0,0.d0,kind(0.d0)) - Ratio_2 = 0.d0 + Ratio_2_array = 0.d0 !X = 1.d0 - Do nf = 1,N_Fl + Do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) DO I = 1,Ndim !X= X*real(Det_vec_new(I,nf),kind(0.d0)) / Real(Det_vec_old(I,nf),kind(0.d0) ) - Ratio_2 = Ratio_2 + Det_vec_new(I,nf) - Det_vec_old(I,nf) + Ratio_2_array(nf) = Ratio_2_array(nf) + Det_vec_new(I,nf_eff) - Det_vec_old(I,nf_eff) enddo + !Z = Z**N_SUN + Ratio_2_array(nf) = real(N_SUN,kind(0.d0))*Ratio_2_array(nf) enddo !Z = cmplx(X,0.d0,kind(0.d0)) - Ratio(1) = cmplx(1.d0,0.d0,kind(0.d0)) - Do nf = 1,N_FL + Ratio_1_array = cmplx(1.d0,0.d0,kind(0.d0)) + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) !Z = Z*Phase_Det_new(nf)/Phase_Det_old(nf) - Ratio(1) = Ratio(1) * Phase_Det_new(nf)/Phase_Det_old(nf) + Ratio_1_array(nf_eff) = Ratio_1_array(nf) * Phase_Det_new(nf_eff)/Phase_Det_old(nf_eff) + !Z = Z**N_SUN + Ratio_1_array(nf) = Ratio_1_array(nf)**N_SUN enddo - !Z = Z**N_SUN - Ratio(1) = Ratio(1)**N_SUN - Ratio_2 = real(N_SUN,kind(0.d0))*Ratio_2 + Ratio(1)=1.d0 Do I = 1,Size(Op_V,1) X = 0.d0 Do nt = 1,Ltrot @@ -675,16 +712,22 @@ Module Global_mod Ratio(1) = Ratio(1) * cmplx( nsigma%Gama(i,nt)/nsigma_old%Gama(i,nt),0.d0,kind(0.d0) ) !You could put this in Ratio_2 X = X + nsigma%Phi(i,nt) - nsigma_old%Phi(i,nt) Enddo - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) !Z = Z * exp(cmplx( X*Real(N_SUN,Kind(0.d0)), 0.d0,kind(0.d0)) * Op_V(i,nf)%g * Op_V(i,nf)%alpha ) - Ratio(1) = Ratio(1) * exp(cmplx( X*Real(N_SUN,Kind(0.d0)), 0.d0,kind(0.d0)) * Op_V(i,nf)%g * Op_V(i,nf)%alpha ) + Ratio_1_array(nf) = Ratio_1_array(nf) * exp(cmplx( X*Real(N_SUN,Kind(0.d0)), 0.d0,kind(0.d0)) * Op_V(i,nf)%g * Op_V(i,nf)%alpha ) Enddo Enddo + if (reconstruction_needed) then + call ham%weight_reconstruction(Ratio_1_array) + call ham%weight_reconstruction(Ratio_2_array) + endif + Ratio(1)=Ratio(1)*product(Ratio_1_array) !Z = Z * cmplx( ham%Delta_S0_global(Nsigma_old),0.d0,kind(0.d0) ) !Z = Z * cmplx( T0_Proposal_ratio, 0.d0,kind(0.d0)) - Ratio_2 = Ratio_2 + log(ham%Delta_S0_global(Nsigma_old)) + log(T0_Proposal_ratio) + Ratio(2) = sum(Ratio_2_array) + Ratio(2) = Ratio(2) + log(ham%Delta_S0_global(Nsigma_old)) + log(T0_Proposal_ratio) - Ratio(2) = Ratio_2 Compute_Ratio_Global = Ratio(1)*exp(Ratio(2)) @@ -746,13 +789,14 @@ Module Global_mod !> Local variables - Integer :: N_size, NCON, I, J, N_part, info, NSTM, N, nf, nst, nt, nt1 + Integer :: N_size, NCON, I, J, N_part, info, NSTM, N, nf, nst, nt, nt1, nf_eff Integer, allocatable :: ipiv(:) COMPLEX (Kind=Kind(0.d0)) :: alpha,beta, Z, Z1 TYPE(UDV_State) :: udvlocal COMPLEX (Kind=Kind(0.d0)), Dimension(:,:), Allocatable :: TP!, U, V COMPLEX (Kind=Kind(0.d0)), Dimension(:), Allocatable :: D + !! TODO adapt to flavor symmetry if(udvl(1)%side .ne. "L" .and. udvl(1)%side .ne. "l" ) then write(error_unit,*) "Compute_Fermion_Det: calling wrong decompose" @@ -761,11 +805,12 @@ Module Global_mod NSTM = Size(udvst, 1) If (storage == "Empty" ) then - DO nf = 1,N_FL + DO nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvl(nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%reset('l') + CALL udvl(nf_eff)%reset('l') endif ENDDO DO NST = NSTM-1,1,-1 @@ -773,8 +818,8 @@ Module Global_mod NT = Stab_nt(NST ) !Write(6,*) 'Call wrapul', NT1,NT, udvl(1)%d(1), udvl(1)%N_part, udvl(1)%Ndim CALL WRAPUL(NT1,NT, udvl) - Do nf = 1,N_FL - udvst(NST, nf) = udvl(nf) + Do nf_eff = 1,N_FL_eff + udvst(NST, nf_eff) = udvl(nf_eff) ENDDO ENDDO NT1 = stab_nt(1) @@ -786,27 +831,28 @@ Module Global_mod N_size=udvl(1)%ndim Det_vec = 0.d0 Allocate (TP(N_part,N_part), ipiv(N_part)) - do nf=1,N_FL + do nf_eff=1,N_FL_eff + nf=Calc_Fl_map(nf_eff) N_part=udvst(1,nf)%N_part do i=1,NSTM-1 do n=1,N_part #if !defined(STABLOG) - Det_Vec(n,nf)=Det_Vec(n,nf)+log(dble(udvst(i,nf)%D(n))) + Det_Vec(n,nf_eff)=Det_Vec(n,nf_eff)+log(dble(udvst(i,nf)%D(n))) #else - Det_Vec(n,nf)=Det_Vec(n,nf)+udvst(i,nf)%L(n) + Det_Vec(n,nf_eff)=Det_Vec(n,nf_eff)+udvst(i,nf)%L(n) #endif enddo enddo do n=1,N_part #if !defined(STABLOG) - Det_Vec(n,nf)=Det_Vec(n,nf)+log(dble(udvl(nf)%D(n))) + Det_Vec(n,nf_eff)=Det_Vec(n,nf_eff)+log(dble(udvl(nf_eff)%D(n))) #else - Det_Vec(n,nf)=Det_Vec(n,nf)+udvl(nf)%L(n) + Det_Vec(n,nf_eff)=Det_Vec(n,nf_eff)+udvl(nf_eff)%L(n) #endif enddo alpha=1.d0 beta=0.d0 - CALL ZGEMM('C','N',N_part,N_part,N_size,alpha,udvl(nf)%U(1,1),N_size,WF_R(nf)%P(1,1),N_size,beta,TP(1,1),N_part) + CALL ZGEMM('C','N',N_part,N_part,N_size,alpha,udvl(nf_eff)%U(1,1),N_size,WF_R(nf)%P(1,1),N_size,beta,TP(1,1),N_part) ! ZGETRF computes an LU factorization of a general M-by-N matrix A ! using partial pivoting with row interchanges. call ZGETRF(N_part, N_part, TP(1,1), N_part, ipiv, info) @@ -817,7 +863,7 @@ Module Global_mod endif Z = Z * TP(J,J) enddo - Phase_det(nf) = Z/abs(Z) + Phase_det(nf_eff) = Z/abs(Z) Det_vec(1,nf) = Det_vec(1,nf)+log(abs(Z)) enddo Deallocate(TP,ipiv) @@ -831,28 +877,29 @@ Module Global_mod beta = cmplx(0.d0,0.d0,kind(0.d0)) Allocate (TP(N_Size,N_Size),D(N_size)) CALL udvlocal%alloc(N_size) - Do nf = 1,N_FL - TP = udvl(nf)%U !udvl stores U^dag instead of U !CT(udvl%U) + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) + TP = udvl(nf_eff)%U !udvl stores U^dag instead of U !CT(udvl%U) #if !defined(STABLOG) #if !defined(STAB3) DO J = 1,N_size - TP(:,J) = TP(:,J) + udvl(nf)%V(:,J)*udvl(nf)%D(J) + TP(:,J) = TP(:,J) + udvl(nf_eff)%V(:,J)*udvl(nf_eff)%D(J) ENDDO #else DO J = 1,N_size - if ( dble(udvl(nf)%D(J)) <= 1.d0 ) then - TP(:,J) = TP(:,J) + udvl(nf)%V(:,J)*udvl(nf)%D(J) + if ( dble(udvl(nf_eff)%D(J)) <= 1.d0 ) then + TP(:,J) = TP(:,J) + udvl(nf_eff)%V(:,J)*udvl(nf_eff)%D(J) else - TP(:,J) = TP(:,J)/udvl(nf)%D(J) + udvl(nf)%V(:,J) + TP(:,J) = TP(:,J)/udvl(nf_eff)%D(J) + udvl(nf_eff)%V(:,J) endif ENDDO #endif #else DO J = 1,N_size - if ( udvl(nf)%L(J) <= 0.d0 ) then - TP(:,J) = TP(:,J) + udvl(nf)%V(:,J)*cmplx(exp(udvl(nf)%L(J)),0.d0,kind(0.d0)) + if ( udvl(nf_eff)%L(J) <= 0.d0 ) then + TP(:,J) = TP(:,J) + udvl(nf_eff)%V(:,J)*cmplx(exp(udvl(nf_eff)%L(J)),0.d0,kind(0.d0)) else - TP(:,J) = TP(:,J)*cmplx(exp(-udvl(nf)%L(J)),0.d0,kind(0.d0)) + udvl(nf)%V(:,J) + TP(:,J) = TP(:,J)*cmplx(exp(-udvl(nf_eff)%L(J)),0.d0,kind(0.d0)) + udvl(nf_eff)%V(:,J) endif ENDDO #endif @@ -861,33 +908,33 @@ Module Global_mod Call UDV_WRAP(TP,udvlocal%U, D, udvlocal%V, NCON ) Z = DET_C(udvlocal%V, N_size) ! Det destroys its argument ! Call MMULT(TP, udvl%U, udvlocal%U) - CALL ZGEMM('C', 'N', N_size, N_size, N_size, alpha, udvl(nf)%U(1,1), N_size, udvlocal%U(1,1), N_size, beta, TP, N_size) + CALL ZGEMM('C', 'N', N_size, N_size, N_size, alpha, udvl(nf_eff)%U(1,1), N_size, udvlocal%U(1,1), N_size, beta, TP, N_size) Z1 = Det_C(TP, N_size) - Phase_det(nf) = Z*Z1/ABS(Z*Z1) + Phase_det(nf_eff) = Z*Z1/ABS(Z*Z1) #if !defined(STABLOG) #if !defined(STAB3) - Det_vec(:,nf) = log(real(D(:))) - Det_vec(1,nf) = log(real(D(1))) + log(ABS(Z*Z1)) + Det_vec(:,nf_eff) = log(real(D(:))) + Det_vec(1,nf_eff) = log(real(D(1))) + log(ABS(Z*Z1)) #else - Det_vec(1,nf) = log(real(D(1))*ABS(Z*Z1)) - if (dble(udvl(nf)%D(1)) > 1.d0) Det_vec(1,nf)=Det_Vec(1,nf)+log(dble(udvl(nf)%D(1))) + Det_vec(1,nf_eff) = log(real(D(1))*ABS(Z*Z1)) + if (dble(udvl(nf_eff)%D(1)) > 1.d0) Det_vec(1,nf_eff)=Det_Vec(1,nf_eff)+log(dble(udvl(nf_eff)%D(1))) Do J=2,Ndim - if (dble(udvl(nf)%D(J))<=1.d0) then - Det_vec(J,nf) = log(real(D(J))) + if (dble(udvl(nf_eff)%D(J))<=1.d0) then + Det_vec(J,nf_eff) = log(real(D(J))) else - Det_vec(J,nf) = log(real(D(J)))+log(dble(udvl(nf)%D(J))) + Det_vec(J,nf_eff) = log(real(D(J)))+log(dble(udvl(nf_eff)%D(J))) endif enddo #endif #else - Det_vec(:,nf) = log(real(D(:))) - Det_vec(1,nf) = log(real(D(1))*ABS(Z*Z1)) - if (udvl(nf)%L(1) > 0.d0) Det_vec(1,nf)=Det_Vec(1,nf)+udvl(nf)%L(1) + Det_vec(:,nf_eff) = log(real(D(:))) + Det_vec(1,nf_eff) = log(real(D(1))*ABS(Z*Z1)) + if (udvl(nf_eff)%L(1) > 0.d0) Det_vec(1,nf_eff)=Det_Vec(1,nf_eff)+udvl(nf_eff)%L(1) Do J=2,Ndim - if (udvl(nf)%L(J)<=0.d0) then - Det_vec(J,nf) = log(real(D(J))) + if (udvl(nf_eff)%L(J)<=0.d0) then + Det_vec(J,nf_eff) = log(real(D(J))) else - Det_vec(J,nf) = log(real(D(J)))+udvl(nf)%L(J) + Det_vec(J,nf_eff) = log(real(D(J)))+udvl(nf_eff)%L(J) endif enddo #endif diff --git a/Prog/Hamiltonian_main_mod.F90 b/Prog/Hamiltonian_main_mod.F90 index 21100c64eec31c0b681a9c9790df891d57e3769e..baed7cdf286e174b37f3f4f97dac91d5847a0af2 100644 --- a/Prog/Hamiltonian_main_mod.F90 +++ b/Prog/Hamiltonian_main_mod.F90 @@ -152,6 +152,9 @@ procedure, nopass :: Delta_S0_global => Delta_S0_global_base procedure, nopass :: S0 => S0_base procedure, nopass :: Ham_Langevin_HMC_S0 => Ham_Langevin_HMC_S0_base + procedure, nopass :: weight_reconstruction => weight_reconstruction_base + procedure, nopass :: GR_reconstruction => GR_reconstruction_base + procedure, nopass :: GRT_reconstruction => GRT_reconstruction_base #ifdef HDF5 procedure, nopass :: write_parameters_hdf5 => write_parameters_hdf5_base #endif @@ -163,15 +166,18 @@ Type (Operator), dimension(:,:), allocatable, public :: Op_T Type (WaveFunction), dimension(:), allocatable, public :: WF_L Type (WaveFunction), dimension(:), allocatable, public :: WF_R + Logical , dimension(:), allocatable, public :: Calc_Fl + Integer , dimension(:), allocatable, public :: Calc_Fl_map Type (Fields), public :: nsigma Integer , public :: Ndim - Integer , public :: N_FL + Integer , public :: N_FL, N_FL_eff Integer , public :: N_SUN Integer , public :: Ltrot Integer , public :: Thtrot Logical , public :: Projector Integer , public :: Group_Comm Logical , public :: Symm + Logical , public :: reconstruction_needed !> Privat Observables @@ -201,7 +207,7 @@ subroutine Alloc_Ham() Implicit none - Integer :: ierr + Integer :: ierr, I Character (len=64) :: ham_name NAMELIST /VAR_HAM_NAME/ ham_name @@ -607,16 +613,78 @@ Forces_0 = 0.d0 end Subroutine Ham_Langevin_HMC_S0_base - - -#ifdef HDF5 - subroutine write_parameters_hdf5_base(filename) + +!-------------------------------------------------------------------- +!> @brief +!> Reconstructs dependent flavors of the configuration's weight. +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!-------------------------------------------------------------------- + subroutine weight_reconstruction_base(weight) implicit none + complex (Kind=Kind(0.d0)), Intent(inout) :: weight(:) + write(error_unit, *) 'weight_reconstruction not defined!' + error stop 1 + end subroutine weight_reconstruction_base + + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Reconstructs dependent flavors of equal time Greens function +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!> @param [INOUT] Gr Complex(:,:,:) +!> \verbatim +!> Green function: Gr(I,J,nf) = on time slice ntau +!> \endverbatim +!------------------------------------------------------------------- + subroutine GR_reconstruction_base(GR) + + Implicit none - Character (len=64), intent(in) :: filename + Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) - end subroutine write_parameters_hdf5_base -#endif + write(error_unit, *) "Warning: GR_reconstruction not implemented." + error stop 1 + end Subroutine GR_reconstruction_base +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Reconstructs dependent flavors of time displaced Greens function G0T and GT0 +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!> @param [INOUT] GT0, G0T, Complex(:,:,:) +!> \verbatim +!> Green functions: +!> GT0(I,J,nf) = +!> G0T(I,J,nf) = +!> \endverbatim +!------------------------------------------------------------------- + Subroutine GRT_reconstruction_base(GT0, G0T) + Implicit none + + Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GT0(Ndim,Ndim,N_FL), G0T(Ndim,Ndim,N_FL) + + write(error_unit, *) "Warning: GRT_reconstruction not implemented." + error stop 1 + end Subroutine GRT_reconstruction_base + + +#ifdef HDF5 + subroutine write_parameters_hdf5_base(filename) + implicit none + + Character (len=64), intent(in) :: filename + + end subroutine write_parameters_hdf5_base +#endif + + end Module Hamiltonian_main diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 index 7e0429b7f41b5f2a2142e6b41b7e13862605808c..4cae85437196e5cf95b06e7e06e9dbbc18732ed2 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 @@ -139,6 +139,9 @@ procedure, nopass :: Alloc_obs procedure, nopass :: Obser procedure, nopass :: ObserT + procedure, nopass :: weight_reconstruction + procedure, nopass :: GR_reconstruction + procedure, nopass :: GRT_reconstruction #ifdef HDF5 procedure, nopass :: write_parameters_hdf5 #endif @@ -165,6 +168,7 @@ Type (Lattice), target :: Latt Type (Unit_cell), target :: Latt_unit + INTEGER :: nf_calc, nf_reconst contains @@ -263,6 +267,12 @@ Write(unit_info,*) 't : ', Ham_T Write(unit_info,*) 'Ham_U : ', Ham_U Write(unit_info,*) 'Ham_chem : ', Ham_chem + If ( Ham_U >=0.d0 .and. Ham_chem == 0.d0 ) then + Write(unit_info,*) 'Assuming particle hole symmetry' + endif + If ( Ham_U < 0.d0 ) then + Write(unit_info,*) 'Assuming time reversal symmetry' + endif if (Projector) then Do nf = 1,N_FL Write(unit_info,*) 'Degen of right trial wave function: ', WF_R(nf)%Degen @@ -274,7 +284,18 @@ Endif #endif - end Subroutine Ham_Set + ! Particle-hole symmetry for repulsive U only if chemical potential vanishes + ! Time reversal symmetry for attractive U + If ( (Ham_U >= 0.d0 .and. Ham_chem == 0.d0) .or. Ham_U < 0.d0 ) then + allocate(Calc_Fl(N_FL)) + nf_calc=2 + nf_reconst=1 + Calc_Fl(nf_calc)=.True. + Calc_Fl(nf_reconst)=.False. + endif + + + End Subroutine Ham_Set !-------------------------------------------------------------------- !> @author @@ -428,7 +449,7 @@ Op_V(i,nf)%P(1) = I Op_V(i,nf)%O(1,1) = cmplx(1.d0, 0.d0, kind(0.D0)) Op_V(i,nf)%g = X*SQRT(CMPLX(DTAU*ham_U/2.d0, 0.D0, kind(0.D0))) - Op_V(i,nf)%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) + Op_V(i,nf)%alpha = cmplx(-0.5d0, 0.d0, kind(0.D0)) Op_V(i,nf)%type = 2 Call Op_set( Op_V(i,nf) ) Enddo @@ -725,4 +746,104 @@ end Subroutine OBSERT +!-------------------------------------------------------------------- +!> @brief +!> Reconstructs dependent flavors of the configuration's weight. +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!-------------------------------------------------------------------- + subroutine weight_reconstruction(weight) + implicit none + complex (Kind=Kind(0.d0)), Intent(inout) :: weight(:) + + weight(nf_reconst) = conjg(Weight(nf_calc)) + + end subroutine weight_reconstruction + + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Reconstructs dependent flavors of equal time Greens function +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!> @param [INOUT] Gr Complex(:,:,:) +!> \verbatim +!> Green function: Gr(I,J,nf) = on time slice ntau +!> \endverbatim +!------------------------------------------------------------------- + subroutine GR_reconstruction(GR) + + Implicit none + + Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: GR(Ndim,Ndim,N_FL) + Integer :: I,J,imj + real (kind=kind(0.d0)) :: X, ZZ + + If (Ham_U >= 0.d0) then + Do J = 1,Ndim + Do I = 1,Ndim + X=-1.0 + imj = latt%imj(I,J) + if (mod(Latt%list(imj,1)+Latt%list(imj,2),2)==0) X=1.d0 + ! write(*,*) Latt%list(I,:),Latt%list(J,:),mod(Latt%list(imj,1)+Latt%list(imj,2),2), X + ZZ=0.d0 + if (I==J) ZZ=1.d0 + GR(I,J,nf_reconst) = ZZ - X*GR(J,I,nf_calc) + Enddo + Enddo + else + Do J = 1,Ndim + Do I = 1,Ndim + GR(I,J,nf_reconst) = Conjg(GR(I,J,nf_calc)) + Enddo + Enddo + Endif + end Subroutine GR_reconstruction + + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Reconstructs dependent flavors of time displaced Greens function G0T and GT0 +!> @details +!> This has to be overloaded in the Hamiltonian submodule. +!> @param [INOUT] GT0, G0T, Complex(:,:,:) +!> \verbatim +!> Green functions: +!> GT0(I,J,nf) = +!> G0T(I,J,nf) = +!> \endverbatim +!------------------------------------------------------------------- + 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)) :: X, ZZ + + If (Ham_U >= 0.d0) then + Do J = 1,Latt%N + Do I = 1,Latt%N + X=-1.0 + imj = latt%imj(I,J) + if (mod(Latt%list(imj,1)+Latt%list(imj,2),2)==0) X=1.d0 + G0T(I,J,nf_reconst) = -X*conjg(GT0(J,I,nf_calc)) + GT0(I,J,nf_reconst) = -X*conjg(G0T(J,I,nf_calc)) + enddo + enddo + else + Do J = 1,Latt%N + Do I = 1,Latt%N + G0T(I,J,nf_reconst) = conjg(G0T(I,J,nf_calc)) + GT0(I,J,nf_reconst) = conjg(GT0(I,J,nf_calc)) + enddo + enddo + endif + end Subroutine GRT_reconstruction + end submodule ham_Hubbard_Plain_Vanilla_smod diff --git a/Prog/Hop_mod.F90 b/Prog/Hop_mod.F90 index b2427b2213dad15699a88c3063406d848a33f04f..91b8899f1a87d6d60728234bde3afb2768faadb0 100644 --- a/Prog/Hop_mod.F90 +++ b/Prog/Hop_mod.F90 @@ -270,11 +270,12 @@ COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), Intent(Out):: Out COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), Intent(IN):: In - Integer :: nf, nc + Integer :: nf, nc, nf_eff class(ContainerElementBase), pointer :: dummy Out = In - Do nf = 1, size(In,3) + Do nf_eff = 1, N_FL_eff !size(In,3) + nf=Calc_Fl_map(nf_eff) do nc = Ncheck,1,-1 dummy => ExpOpT_vec(nf)%at(nc) call dummy%adjointaction(Out(:, :, nf)) diff --git a/Prog/Langevin_HMC_mod.F90 b/Prog/Langevin_HMC_mod.F90 index b5dc5ae7d3796232c8a43175e4d2cd2d5afb1cc6..c6f1bc4680c84cda674ecadff8ec1ff81b4f7b3d 100644 --- a/Prog/Langevin_HMC_mod.F90 +++ b/Prog/Langevin_HMC_mod.F90 @@ -30,6 +30,8 @@ ! to the ALF project or to mark your material in a reasonable way as different from the original version. +! TODO ATTENTION: this module still has to be updated for flavor symmetries!!! + Module Langevin_HMC_mod Use Hamiltonian_main @@ -103,8 +105,8 @@ !Local - Integer :: NSTM, n, nf, NST, NTAU, nt, nt1, Ntau1, NVAR, N_Type, I, J - Complex (Kind=Kind(0.d0)) :: Z, Z1 + Integer :: NSTM, n, nf, nf_eff, NST, NTAU, nt, nt1, Ntau1, NVAR, N_Type, I, J + Complex (Kind=Kind(0.d0)) :: Z, Z1, Phase_array(N_FL) Real (Kind=Kind(0.d0)) :: spin NSTM = Size(Stab_nt,1) - 1 @@ -113,11 +115,11 @@ !Enddo Langevin_HMC%Forces = cmplx(0.d0,0.d0,Kind(0.d0)) - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff if (Projector) then - CALL udvr(nf)%reset('r',WF_R(nf)%P) + CALL udvr(nf_eff)%reset('r',WF_R(nf_eff)%P) else - CALL udvr(nf)%reset('r') + CALL udvr(nf_eff)%reset('r') endif Enddo NST = 1 @@ -129,18 +131,22 @@ If (NTAU1 == Stab_nt(NST) ) then NT1 = Stab_nt(NST-1) CALL WRAPUR(NT1, NTAU1, udvr) - Z = cmplx(1.d0, 0.d0, kind(0.D0)) - Do nf = 1, N_FL + Phase_array = cmplx(1.d0, 0.d0, kind(0.D0)) + Do nf_eff = 1, N_FL_eff + nf=Calc_FL_map(nf_eff) ! Read from storage left propagation from LTROT to NTAU1 - udvl(nf) = udvst(NST, nf) + udvl(nf_eff) = udvst(NST, nf_eff) NVAR = 1 IF (NTAU1 .GT. LTROT/2) NVAR = 2 TEST(:,:) = GR(:,:,nf) - CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf), UDVL(nf)) - Z = Z*Z1 + CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf_eff), UDVL(nf_eff)) Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) + call Op_phase(Z1,OP_V,Nsigma,nf) + Phase_array(nf)=Z1 ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Z=product(Phase_array) + Z=Z**N_SUN Call Control_PrecisionP(Z,Phase) Phase = Z NST = NST + 1 @@ -149,8 +155,10 @@ IF (NTAU1 .GE. LOBS_ST .AND. NTAU1 .LE. LOBS_EN .and. Calc_Obser_eq ) THEN If (Symm) then Call Hop_mod_Symm(GR_Tilde,GR) + If (reconstruction_needed) Call ham%GR_reconstruction( GR_Tilde ) CALL ham%Obser( GR_Tilde, PHASE, Ntau1,Langevin_HMC%Delta_t_running ) else + If (reconstruction_needed) Call ham%GR_reconstruction( GR ) CALL ham%Obser( GR, PHASE, Ntau1, Langevin_HMC%Delta_t_running ) endif endif @@ -179,36 +187,43 @@ !Local - Complex (Kind=Kind(0.d0)) :: Z, Z1 - Integer :: nf, I, J, n, N_type + Complex (Kind=Kind(0.d0)) :: Z(N_FL), Z1 + Integer :: nf, I, J, n, N_type, nf_eff Real(Kind=Kind(0.d0)) :: spin - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_FL_map(nf_eff) CALL HOP_MOD_mmthr (GR(:,:,nf), nf ) CALL HOP_MOD_mmthl_m1(GR(:,:,nf), nf ) Enddo Do n = 1, size(OP_V,1) this%Forces(n,nt1) = cmplx(0.d0,0.d0,Kind(0.d0)) - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_FL_map(nf_eff) spin = nsigma%f(n,nt1) ! Phi(nsigma(n,ntau1),Op_V(n,nf)%type) N_type = 1 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) enddo + !TODO how does flavor symmetry effect section below? I feel like skipping some flavors is incorrect! if (OP_V(n,1)%type == 3 ) then - Do nf = 1, N_Fl - Z = cmplx(0.d0,0.d0,Kind(0.d0)) + Z = cmplx(0.d0,0.d0,Kind(0.d0)) + Do nf_eff = 1, N_Fl_eff + nf=Calc_FL_map(nf_eff) do I = 1,size(OP_V(n,nf)%P,1) do J = 1,size(OP_V(n,nf)%P,1) Z1 = cmplx(0.d0,0.d0,Kind(0.d0)) if ( I == J ) Z1 = cmplx(1.d0,0.d0,Kind(0.d0)) - Z = Z + Op_V(n,nf)%O(I,J) * ( Z1 - Gr(Op_V(n,nf)%P(J),Op_V(n,nf)%P(I), nf) ) + Z(nf) = Z(nf) + Op_V(n,nf)%O(I,J) * ( Z1 - Gr(Op_V(n,nf)%P(J),Op_V(n,nf)%P(I), nf) ) Enddo Enddo + Enddo + if (reconstruction_needed) call ham%weight_reconstruction(Z) + Do nf = 1, N_Fl this%Forces(n,nt1) = this%Forces(n,nt1) - & - & Op_V(n,nf)%g * Z * cmplx(real(N_SUN,Kind(0.d0)), 0.d0, Kind(0.d0)) + & Op_V(n,nf)%g * Z(nf) * cmplx(real(N_SUN,Kind(0.d0)), 0.d0, Kind(0.d0)) Enddo endif enddo @@ -235,18 +250,19 @@ Integer, intent(in), dimension(:), allocatable :: Stab_nt ! Local - Integer :: NSTM, nf, nt, nt1, NST, NVAR - Complex (Kind=Kind(0.d0)) :: Z + Integer :: NSTM, nf, nt, nt1, NST, NVAR, nf_eff + Complex (Kind=Kind(0.d0)) :: Z, Phase_array(N_FL) NSTM = Size(Stab_nt,1) - 1 - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_FL_map(nf_eff) if (Projector) then - CALL udvl(nf)%reset('l',WF_L(nf)%P) - CALL udvst(NSTM, nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%reset('l',WF_L(nf)%P) + CALL udvst(NSTM, nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%reset('l') - CALL udvst(NSTM, nf)%reset('l') + CALL udvl(nf_eff)%reset('l') + CALL udvst(NSTM, nf_eff)%reset('l') endif ENDDO @@ -256,28 +272,33 @@ NT = Stab_nt(NST ) !Write(6,*)'Hi', NT1,NT, NST CALL WRAPUL(NT1, NT, UDVL) - Do nf = 1,N_FL - UDVST(NST, nf) = UDVL(nf) + Do nf_eff = 1,N_FL_eff + UDVST(NST, nf_eff) = UDVL(nf_eff) ENDDO ENDDO NT1 = stab_nt(1) CALL WRAPUL(NT1, 0, UDVL) - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_FL_map(nf_eff) if (Projector) then - CALL udvr(nf)%reset('r',WF_R(nf)%P) + CALL udvr(nf_eff)%reset('r',WF_R(nf)%P) else - CALL udvr(nf)%reset('r') + CALL udvr(nf_eff)%reset('r') endif ENDDO NVAR = 1 - Phase = cmplx(1.d0, 0.d0, kind(0.D0)) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), UDVR(nf), UDVL(nf)) - Phase = Phase*Z + Phase_array = cmplx(1.d0, 0.d0, kind(0.D0)) + do nf_eff = 1,N_Fl_eff + nf=Calc_FL_map(nf_eff) + CALL CGR(Z, NVAR, GR(:,:,nf), UDVR(nf_eff), UDVL(nf_eff)) + call Op_phase(Z,OP_V,Nsigma,nf) + Phase_array(nf)=Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase=product(Phase_array) + Phase=Phase**N_SUN end Subroutine Langevin_HMC_Reset_storage diff --git a/Prog/Operator_mod.F90 b/Prog/Operator_mod.F90 index 6f27dc379675732186859224198640d5335d06fb..c94bfc425767f5fdcfb6f0713712fd4a42333a39 100644 --- a/Prog/Operator_mod.F90 +++ b/Prog/Operator_mod.F90 @@ -147,30 +147,27 @@ Contains !> @param[in] Nsigma !> Type(Fields) !> * Fields -!> @param[in] N_SUN +!> @param[in] nf !> Integer -!> * Number of colors +!> * flavor index !-------------------------------------------------------------------- - Subroutine Op_phase(Phase,OP_V,Nsigma,N_SUN) + Subroutine Op_phase(Phase,OP_V,Nsigma,nf) Implicit none Complex (Kind=Kind(0.d0)), Intent(Inout) :: Phase - Integer, Intent(IN) :: N_SUN + Integer, Intent(IN) :: nf Type (Fields), Intent(IN) :: Nsigma Type (Operator), dimension(:,:), Intent(In) :: Op_V Real (Kind=Kind(0.d0)) :: angle - Integer :: n, nf, nt + Integer :: n, nt - do nf = 1,Size(Op_V,2) - do n = 1,size(Op_V,1) - do nt = 1,size(nsigma%f,2) - angle = Aimag( Op_V(n,nf)%g * Op_V(n,nf)%alpha ) * nsigma%Phi(n,nt) - Phase = Phase*CMPLX(cos(angle),sin(angle), Kind(0.D0)) - enddo + do n = 1,size(Op_V,1) + do nt = 1,size(nsigma%f,2) + angle = Aimag( Op_V(n,nf)%g * Op_V(n,nf)%alpha ) * nsigma%Phi(n,nt) + Phase = Phase*CMPLX(cos(angle),sin(angle), Kind(0.D0)) enddo enddo - Phase = Phase**N_SUN end Subroutine Op_phase diff --git a/Prog/Predefined_Int_mod.F90 b/Prog/Predefined_Int_mod.F90 index b30b4f4505c9a5d854eede6b169db78eaf5d8997..09527316c4f28188af168e9034c12940dfcca39d 100644 --- a/Prog/Predefined_Int_mod.F90 +++ b/Prog/Predefined_Int_mod.F90 @@ -113,7 +113,7 @@ ! !> @brief !> Sets Hubbard U Mz -!> - U/2 ( n_{i,up} - n_{i,do} ) ^2 +!> - U/2 ( (n_{i,up} -1/2) - (n_{i,do} - 1/2) ) ^2 !> Here N_FL = 2 and the routine sets both the up and down operators. !> The routine uses a descrete HS transformation with four fields !> (see documentation) @@ -132,12 +132,14 @@ Op_up%O(1,1) = cmplx(1.d0 ,0.d0, kind(0.D0)) Op_up%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) Op_up%g = SQRT(CMPLX(DTAU*U/2.d0, 0.D0, kind(0.D0))) + Op_up%alpha = cmplx(-0.5d0,0.d0, kind(0.D0)) Op_up%type = 2 Op_do%P(1) = I Op_do%O(1,1) = cmplx(1.d0 ,0.d0, kind(0.D0)) Op_do%alpha = cmplx(0.d0, 0.d0, kind(0.D0)) Op_do%g = -SQRT(CMPLX(DTAU*U/2.d0, 0.D0, kind(0.D0))) + Op_do%alpha = cmplx(-0.5d0,0.d0, kind(0.D0)) Op_do%type = 2 Call Op_set( Op_up ) diff --git a/Prog/Wrapgr_mod.F90 b/Prog/Wrapgr_mod.F90 index 2d6b01f1ce3619ee7a57e9f703bceef14f1a5311..b4b34ac9a73d65ea482094d05e6cea9e2757d3d1 100644 --- a/Prog/Wrapgr_mod.F90 +++ b/Prog/Wrapgr_mod.F90 @@ -99,7 +99,7 @@ Contains INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau !Local - Integer :: nf, N_Type, NTAU1,n, m + Integer :: nf, nf_eff, N_Type, NTAU1,n, m Complex (Kind=Kind(0.d0)) :: Prev_Ratiotot Real (Kind=Kind(0.d0)) :: T0_proposal, T0_Proposal_ratio, S0_ratio, spin, HS_new Character (Len=64) :: Mode @@ -107,12 +107,14 @@ Contains ! Wrap up, upgrade ntau1. with B^{1}(tau1) NTAU1 = NTAU + 1 - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) CALL HOP_MOD_mmthr (GR(:,:,nf), nf ) CALL HOP_MOD_mmthl_m1(GR(:,:,nf), nf ) Enddo Do n = Nt_sequential_start,Nt_sequential_end - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) spin = nsigma%f(n,ntau1) ! Phi(nsigma(n,ntau1),Op_V(n,nf)%type) N_type = 1 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) @@ -137,7 +139,8 @@ Contains toggle1 = .false. Call Control_upgrade_eff(toggle1) Endif - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) enddo @@ -177,7 +180,7 @@ Contains INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau ! Local - Integer :: nf, N_Type, n, m + Integer :: nf, nf_eff, N_Type, n, m Complex (Kind=Kind(0.d0)) :: Prev_Ratiotot Real (Kind=Kind(0.d0)) :: T0_proposal, T0_Proposal_ratio, S0_ratio, spin, HS_new Character (Len=64) :: Mode @@ -196,7 +199,8 @@ Contains N_type = 2 nf = 1 spin = nsigma%f(n,ntau) - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type) enddo !Write(6,*) 'Upgrade : ', ntau,n @@ -225,11 +229,13 @@ Contains nf = 1 spin = nsigma%f(n,ntau) ! Phi(nsigma(n,ntau),Op_V(n,nf)%type) N_type = 1 - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type ) enddo enddo - DO nf = 1,N_FL + DO nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call Hop_mod_mmthl (GR(:,:,nf), nf) Call Hop_mod_mmthr_m1(GR(:,:,nf), nf) enddo @@ -262,7 +268,7 @@ Contains Integer, INTENT(IN) :: m,m1, ntau !Local - Integer :: n, nf, N_Type + Integer :: n, nf, nf_eff, N_Type Real (Kind=Kind(0.d0)) :: spin If (m == m1) then @@ -270,12 +276,14 @@ Contains elseif ( m1 > m ) then !Write(6,*) "Wrapup from ", m + 1, "to", m1, " on tau=", ntau Do n = m+1,m1 - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) spin = nsigma%f(n,ntau) N_type = 1 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) enddo - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) enddo @@ -286,14 +294,16 @@ Contains N_type = 2 nf = 1 spin = nsigma%f(n,ntau) - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type) enddo !Write(6,*) 'Upgrade : ', ntau,n nf = 1 spin = nsigma%f(n,ntau) N_type = 1 - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), spin, Ndim, N_Type ) enddo enddo @@ -343,7 +353,7 @@ Contains ! Space for local variables - Integer :: n, Flip_length, nf, N_Type, ng_c, Flip_count + 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, HS_new, spin COMPLEX (Kind=Kind(0.d0)) :: Prev_Ratiotot Logical :: Acc @@ -374,7 +384,8 @@ Contains Call Wrapgr_PlaceGR(GR,m, n-1, ntau) !Write(6,*) "Back from PlaceGR", m, n-1,ntau If ( Flip_count == 1 .and. Flip_length > 1 ) GR_st = Gr - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) spin = nsigma%f(n,ntau) N_type = 1 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) @@ -385,7 +396,8 @@ Contains HS_new = Flip_value(Flip_count) Call Upgrade2(GR,n,ntau,PHASE, HS_new , & & Prev_Ratiotot, S0_ratio, T0_Proposal_ratio, Acc, mode ) - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),spin,Ndim,N_Type) enddo @@ -397,7 +409,8 @@ Contains & Prev_Ratiotot, S0_ratio, T0_Proposal_ratio, Acc, mode ) !Write(6,*) "Back from up mode final", n,ntau !Write(6,*) "Acceptance", Acc - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),Spin,Ndim,N_Type) enddo @@ -472,7 +485,7 @@ Contains Integer :: ntau ! Local - Integer :: m, m1, N_size, nf, nth + Integer :: m, m1, N_size, nf, nf_eff, nth Real (Kind=kind(0.d0)) :: Xmean, Xmax !Input is the Green function on time slice tau @@ -490,7 +503,8 @@ Contains Write(6,*) m, N_size Xmax = 0.d0 - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) Call COMPARE(GR_st(:,:,nf),GR(:,:,nf),XMAX,XMEAN) Enddo Write(6,*) 'Compare Global_tau_mod_Test ', Xmax diff --git a/Prog/main.F90 b/Prog/main.F90 index 1d6ea6aa59b90752cba19bbcc61c5353c6c76aaa..94e5804df3686c39d004e2aa914f8db122f8973e 100644 --- a/Prog/main.F90 +++ b/Prog/main.F90 @@ -146,6 +146,7 @@ Program Main COMPLEX (Kind=Kind(0.d0)), Dimension(:,:) , Allocatable :: TEST COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), Allocatable :: GR, GR_Tilde CLASS(UDV_State), DIMENSION(:), ALLOCATABLE :: udvl, udvr + COMPLEX (Kind=Kind(0.d0)), Dimension(:) , Allocatable :: Phase_array Integer :: Nwrap, NSweep, NBin, NBin_eff,Ltau, NSTM, NT, NT1, NVAR, LOBS_EN, LOBS_ST, NBC, NSW @@ -184,7 +185,7 @@ Program Main ! General - Integer :: Ierr, I,nf, nst, n, N_op + Integer :: Ierr, I,nf, nf_eff, nst, n, N_op Complex (Kind=Kind(0.d0)) :: Z_ONE = cmplx(1.d0, 0.d0, kind(0.D0)), Phase, Z, Z1 Real (Kind=Kind(0.d0)) :: ZERO = 10D-8, X, X1 Real (Kind=Kind(0.d0)) :: Mc_step_weight @@ -317,12 +318,36 @@ Program Main Call Fields_init() Call Alloc_Ham() Call ham%Ham_set() + + ! Test if user has initialized Calc_FL array + If ( .not. allocated(Calc_Fl)) then + allocate(Calc_Fl(N_FL)) + Calc_Fl=.True. + endif + ! Count number of flavors to be calculated + N_FL_eff=0 + Do I=1,N_Fl + if (Calc_Fl(I)) N_FL_eff=N_FL_eff+1 + Enddo + reconstruction_needed=.false. + If (N_FL_eff /= N_FL) reconstruction_needed=.true. + !initialize the flavor map + allocate(Calc_Fl_map(N_FL_eff),Phase_array(N_FL)) + N_FL_eff=0 + Do I=1,N_Fl + if (Calc_Fl(I)) then + N_FL_eff=N_FL_eff+1 + Calc_Fl_map(N_FL_eff)=I + endif + Enddo + if(Projector) then if (.not. allocated(WF_R) .or. .not. allocated(WF_L)) then write(error_unit,*) "Projector is selected but there are no trial wave functions!" error stop 1 endif - do nf=1,N_fl + do nf_eff=1,N_fl_eff + nf=Calc_Fl_map(nf_eff) if (.not. allocated(WF_R(nf)%P) .or. .not. allocated(WF_L(nf)%P)) then write(error_unit,*) "Projector is selected but there are no trial wave functions!" error stop 1 @@ -513,19 +538,20 @@ Program Main !Call Test_Hamiltonian Allocate ( Test(Ndim,Ndim), GR(NDIM,NDIM,N_FL), GR_Tilde(NDIM,NDIM,N_FL) ) - ALLOCATE(udvl(N_FL), udvr(N_FL), udvst(NSTM, N_FL)) - do nf = 1, N_FL + ALLOCATE(udvl(N_FL_eff), udvr(N_FL_eff), udvst(NSTM, N_FL_eff)) + do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) do n = 1, NSTM - CALL udvst(n,nf)%alloc(ndim) + CALL udvst(n,nf_eff)%alloc(ndim) ENDDO if (Projector) then - CALL udvl(nf)%init(ndim,'l',WF_L(nf)%P) - CALL udvr(nf)%init(ndim,'r',WF_R(nf)%P) - CALL udvst(NSTM, nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%init(ndim,'l',WF_L(nf)%P) + CALL udvr(nf_eff)%init(ndim,'r',WF_R(nf)%P) + CALL udvst(NSTM, nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%init(ndim,'l') - CALL udvr(nf)%init(ndim,'r') - CALL udvst(NSTM, nf)%reset('l') + CALL udvl(nf_eff)%init(ndim,'l') + CALL udvr(nf_eff)%init(ndim,'r') + CALL udvst(NSTM, nf_eff)%reset('l') endif enddo @@ -534,8 +560,8 @@ Program Main NT = Stab_nt(NST ) !Write(6,*)'Hi', NT1,NT, NST CALL WRAPUL(NT1, NT, UDVL) - Do nf = 1,N_FL - UDVST(NST, nf) = UDVL(nf) + Do nf_eff = 1,N_FL_eff + UDVST(NST, nf_eff) = UDVL(nf_eff) ENDDO ENDDO NT1 = stab_nt(1) @@ -544,12 +570,16 @@ Program Main NVAR = 1 - Phase = cmplx(1.d0, 0.d0, kind(0.D0)) - do nf = 1,N_Fl - CALL CGR(Z, NVAR, GR(:,:,nf), UDVR(nf), UDVL(nf)) - Phase = Phase*Z + Phase_array = cmplx(1.d0, 0.d0, kind(0.D0)) + do nf_eff = 1,N_Fl_eff + nf=Calc_Fl_map(nf_eff) + CALL CGR(Z, NVAR, GR(:,:,nf), UDVR(nf_eff), UDVL(nf_eff)) + call Op_phase(Z,OP_V,Nsigma,nf) + Phase_array(nf)=Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Phase=product(Phase_array) + Phase=Phase**N_SUN #ifdef MPI !WRITE(6,*) 'Phase is: ', Irank, PHASE, GR(1,1,1) #else @@ -603,11 +633,12 @@ Program Main If (Sequential) then ! Propagation from 1 to Ltrot ! Set the right storage to 1 - do nf = 1,N_FL + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvr(nf)%reset('r',WF_R(nf)%P) + CALL udvr(nf_eff)%reset('r',WF_R(nf)%P) else - CALL udvr(nf)%reset('r') + CALL udvr(nf_eff)%reset('r') endif Enddo @@ -619,20 +650,24 @@ Program Main If (NTAU1 == Stab_nt(NST) ) then NT1 = Stab_nt(NST-1) CALL WRAPUR(NT1, NTAU1, udvr) - Z = cmplx(1.d0, 0.d0, kind(0.D0)) - Do nf = 1, N_FL + Phase_array = cmplx(1.d0, 0.d0, kind(0.D0)) + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) ! Read from storage left propagation from LTROT to NTAU1 - udvl(nf) = udvst(NST, nf) + udvl(nf_eff) = udvst(NST, nf_eff) ! Write in storage right prop from 1 to NTAU1 - udvst(NST, nf) = udvr(nf) + udvst(NST, nf_eff) = udvr(nf_eff) NVAR = 1 IF (NTAU1 .GT. LTROT/2) NVAR = 2 TEST(:,:) = GR(:,:,nf) - CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf), UDVL(nf)) - Z = Z*Z1 + CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf_eff), UDVL(nf_eff)) Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) + call Op_phase(Z1,OP_V,Nsigma,nf) + Phase_array(nf)=Z1 ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Z=product(Phase_array) + Z=Z**N_SUN Call Control_PrecisionP(Z,Phase) Phase = Z NST = NST + 1 @@ -646,18 +681,23 @@ Program Main Mc_step_weight = 1.d0 If (Symm) then Call Hop_mod_Symm(GR_Tilde,GR) + !reconstruction of NOT calculated block!!! + If (reconstruction_needed) Call ham%GR_reconstruction( GR_Tilde ) CALL ham%Obser( GR_Tilde, PHASE, Ntau1, Mc_step_weight ) else + !reconstruction of NOT calculated block!!! + If (reconstruction_needed) Call ham%GR_reconstruction( GR ) CALL ham%Obser( GR, PHASE, Ntau1, Mc_step_weight ) endif ENDIF ENDDO - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) if (Projector) then - CALL udvl(nf)%reset('l',WF_L(nf)%P) + CALL udvl(nf_eff)%reset('l',WF_L(nf)%P) else - CALL udvl(nf)%reset('l') + CALL udvl(nf_eff)%reset('l') endif ENDDO @@ -671,8 +711,12 @@ Program Main Mc_step_weight = 1.d0 If (Symm) then Call Hop_mod_Symm(GR_Tilde,GR) + !reconstruction of NOT calculated block!!! + If (reconstruction_needed) Call ham%GR_reconstruction( GR_Tilde ) CALL ham%Obser( GR_Tilde, PHASE, Ntau1, Mc_step_weight ) else + !reconstruction of NOT calculated block!!! + If (reconstruction_needed) Call ham%GR_reconstruction( GR ) CALL ham%Obser( GR, PHASE, Ntau1,Mc_step_weight ) endif ENDIF @@ -681,20 +725,24 @@ Program Main !Write(6,*) 'Wrapul : ', NT1, NTAU1 CALL WRAPUL(NT1, NTAU1, udvl) !Write(6,*) 'Write UL, read UR ', NTAU1, NST - Z = cmplx(1.d0, 0.d0, kind(0.D0)) - do nf = 1,N_FL + Phase_array = cmplx(1.d0, 0.d0, kind(0.D0)) + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) ! Read from store the right prop. from 1 to LTROT/NWRAP-1 - udvr(nf) = udvst(NST, nf) + udvr(nf_eff) = udvst(NST, nf_eff) ! WRITE in store the left prop. from LTROT/NWRAP-1 to 1 - udvst(NST, nf) = udvl(nf) + udvst(NST, nf_eff) = udvl(nf_eff) NVAR = 1 IF (NTAU1 .GT. LTROT/2) NVAR = 2 TEST(:,:) = GR(:,:,nf) - CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf), UDVL(nf)) - Z = Z*Z1 + CALL CGR(Z1, NVAR, GR(:,:,nf), UDVR(nf_eff), UDVL(nf_eff)) Call Control_PrecisionG(GR(:,:,nf),Test,Ndim) + call Op_phase(Z1,OP_V,Nsigma,nf) + Phase_array(nf)=Z1 ENDDO - call Op_phase(Z,OP_V,Nsigma,N_SUN) + if (reconstruction_needed) call ham%weight_reconstruction(Phase_array) + Z=product(Phase_array) + Z=Z**N_SUN Call Control_PrecisionP(Z,Phase) Phase = Z IF( LTAU == 1 .and. Projector .and. Stab_nt(NST)<=THTROT+1 .and. THTROT+1 0) then - Allocate ( Mat(Op_dim,Op_Dim), Delta(Op_dim,N_FL), u(Ndim,Op_dim), v(Ndim,Op_dim) ) + Allocate ( Mat(Op_dim,Op_Dim), Delta(Op_dim,N_FL_eff), u(Ndim,Op_dim), v(Ndim,Op_dim) ) Allocate ( y_v(Ndim,Op_dim), xp_v(Ndim,Op_dim), x_v(Ndim,Op_dim) ) endif @@ -155,14 +156,15 @@ module upgrade_mod nf = 1 nsigma_new%f(1,1) = Hs_new !real(ns_new,kind=kind(0.d0)) nsigma_new%t(1) = Op_V(n_op,nf)%Type - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) !Z1 = Op_V(n_op,nf)%g * ( Phi_st(ns_new,Op_V(n_op,nf)%type) - Phi_st(ns_old,Op_V(n_op,nf)%type)) Z1 = Op_V(n_op,nf)%g * ( nsigma_new%Phi(1,1) - nsigma%Phi(n_op,nt) ) op_dim_nf = Op_V(n_op,nf)%N_non_zero Do m = 1,op_dim_nf myexp = exp( Z1* Op_V(n_op,nf)%E(m) ) Z = myexp - 1.d0 - Delta(m,nf) = Z + Delta(m,nf_eff) = Z do n = 1,op_dim_nf Mat(n,m) = - Z * GR( Op_V(n_op,nf)%P(n), Op_V(n_op,nf)%P(m),nf ) Enddo @@ -187,6 +189,9 @@ module upgrade_mod Ratio(nf) = D_Mat * exp( Z1*Op_V(n_op,nf)%alpha ) Enddo + !call reconstruct weight subroutine to fill the non-calculated blocks + if (reconstruction_needed) call ham%weight_reconstruction(Ratio) + Ratiotot = Product(Ratio) nf = 1 !Ratiotot = (Ratiotot**dble(N_SUN)) * Gama_st(ns_new, Op_V(n_op,nf)%type)/Gama_st(ns_old, Op_V(n_op,nf)%type) @@ -215,7 +220,8 @@ module upgrade_mod If (mode == "Final" ) Phase = Phase * Ratiotot/sqrt(Ratiotot*conjg(Ratiotot)) !Write(6,*) 'Accepted : ', Ratiotot - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) ! Setup u(i,n), v(n,i) op_dim_nf = Op_V(n_op,nf)%N_non_zero if (op_dim_nf > 0) then @@ -223,7 +229,7 @@ module upgrade_mod call zlaset('N', Ndim, op_dim_nf, beta, beta, u, size(u, 1)) call zlaset('N', Ndim, op_dim_nf, beta, beta, v, size(v, 1)) do n = 1,op_dim_nf - u( Op_V(n_op,nf)%P(n), n) = Delta(n,nf) + u( Op_V(n_op,nf)%P(n), n) = Delta(n,nf_eff) do i = 1,Ndim v(i,n) = - GR( Op_V(n_op,nf)%P(n), i, nf ) enddo diff --git a/Prog/wrapul_mod.F90 b/Prog/wrapul_mod.F90 index 18160491832538efddede1afe4ac46fa98d26f6f..bd467d5e200de540faf3dee70e6bc04bd9250bdb 100644 --- a/Prog/wrapul_mod.F90 +++ b/Prog/wrapul_mod.F90 @@ -69,7 +69,7 @@ module wrapul_mod ! Working space. COMPLEX (Kind=Kind(0.d0)) :: U1(Ndim,Ndim), V1(Ndim,Ndim), TMP(Ndim,Ndim), TMP1(Ndim,Ndim) COMPLEX (Kind=Kind(0.d0)) :: Z_ONE, beta - Integer :: NT, NCON, n, nf + Integer :: NT, NCON, n, nf, nf_eff Real (Kind=Kind(0.d0)) :: X @@ -77,7 +77,8 @@ module wrapul_mod Z_ONE = cmplx(1.d0, 0.d0, kind(0.D0)) beta = 0.D0 - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) CALL INITD(TMP,Z_ONE) DO NT = NTAU1, NTAU+1 , -1 Do n = Size(Op_V,1),1,-1 @@ -89,16 +90,16 @@ module wrapul_mod ENDDO !Carry out U,D,V decomposition. - CALL ZGEMM('C', 'N', Ndim, UDVL(nf)%N_part, Ndim, Z_ONE, TMP, Ndim, udvl(nf)%U(1, 1), Ndim, beta, TMP1, Ndim) - if( ALLOCATED(UDVL(nf)%V) ) then - DO n = 1,UDVL(nf)%N_part - TMP1(:, n) = TMP1(:, n) * udvl(nf)%D(n) + CALL ZGEMM('C', 'N', Ndim, UDVL(nf_eff)%N_part, Ndim, Z_ONE, TMP, Ndim, udvl(nf_eff)%U(1, 1), Ndim, beta, TMP1, Ndim) + if( ALLOCATED(UDVL(nf_eff)%V) ) then + DO n = 1,UDVL(nf_eff)%N_part + TMP1(:, n) = TMP1(:, n) * udvl(nf_eff)%D(n) ENDDO - CALL UDV_WRAP_Pivot(TMP1,udvl(nf)%U,udvl(nf)%D,V1,NCON,Ndim,Ndim) - CALL ZGEMM('N', 'C', Ndim, Ndim, Ndim, Z_ONE, udvl(nf)%V(1,1), Ndim, V1, Ndim, beta, TMP1, Ndim) - udvl(nf)%V = TMP1 + CALL UDV_WRAP_Pivot(TMP1,udvl(nf_eff)%U,udvl(nf_eff)%D,V1,NCON,Ndim,Ndim) + CALL ZGEMM('N', 'C', Ndim, Ndim, Ndim, Z_ONE, udvl(nf_eff)%V(1,1), Ndim, V1, Ndim, beta, TMP1, Ndim) + udvl(nf_eff)%V = TMP1 else - CALL UDV_WRAP_Pivot(TMP1(:,1:UDVL(nf)%N_part),udvl(nf)%U,udvl(nf)%D,V1,NCON,Ndim,UDVL(nf)%N_part) + CALL UDV_WRAP_Pivot(TMP1(:,1:UDVL(nf_eff)%N_part),udvl(nf_eff)%U,udvl(nf_eff)%D,V1,NCON,Ndim,UDVL(nf_eff)%N_part) endif ENDDO @@ -110,18 +111,19 @@ module wrapul_mod CLASS(UDV_State), intent(inout), allocatable, dimension(:) :: UDVL Integer, intent(in) :: NTAU1, NTAU - Integer :: NT, n, nf + Integer :: NT, n, nf, nf_eff - Do nf = 1, N_FL + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) DO NT = NTAU1, NTAU+1 , -1 Do n = Size(Op_V,1),1,-1 - Call Op_mmultR(udvl(nf)%U,Op_V(n,nf),nsigma%f(n,nt),'c') + Call Op_mmultR(udvl(nf_eff)%U,Op_V(n,nf),nsigma%f(n,nt),'c') enddo - Call Hop_mod_mmthlc (udvl(nf)%U,nf) + Call Hop_mod_mmthlc (udvl(nf_eff)%U,nf) ENDDO !Carry out U,D,V decomposition. - CALL UDVL(nf)%decompose + CALL UDVL(nf_eff)%decompose Enddo #endif END SUBROUTINE WRAPUL diff --git a/Prog/wrapur_mod.F90 b/Prog/wrapur_mod.F90 index 4e2e1ab838d3458b8ccda694ddadbb963a39bf92..dec454e214bb616565312f7773dec9e779b6ce4a 100644 --- a/Prog/wrapur_mod.F90 +++ b/Prog/wrapur_mod.F90 @@ -66,13 +66,14 @@ module wrapur_mod ! Working space. Complex (Kind=Kind(0.d0)) :: Z_ONE COMPLEX (Kind=Kind(0.d0)) :: V1(Ndim,Ndim), TMP(Ndim,Ndim), TMP1(Ndim,Ndim) - Integer ::NCON, NT, I, J, n, nf + Integer ::NCON, NT, I, J, n, nf, nf_eff Real (Kind=Kind(0.d0)) :: X NCON = 0 ! Test for UDV :::: 0: Off, 1: On. Z_ONE = cmplx(1.d0, 0.d0, kind(0.D0)) - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) CALL INITD(TMP,Z_ONE) DO NT = NTAU + 1, NTAU1 !CALL MMULT(TMP1,Exp_T(:,:,nf) ,TMP) @@ -83,17 +84,17 @@ module wrapur_mod Call Op_mmultR(Tmp,Op_V(n,nf),nsigma%f(n,nt),'n') ENDDO ENDDO - CALL MMULT(TMP1,TMP, udvr(nf)%U) - if(allocated(udvr(nf)%V)) then + CALL MMULT(TMP1,TMP, udvr(nf_eff)%U) + if(allocated(udvr(nf_eff)%V)) then DO J = 1,NDim DO I = 1,NDim - TMP1(I,J) = TMP1(I,J)*udvr(nf)%D(J) + TMP1(I,J) = TMP1(I,J)*udvr(nf_eff)%D(J) ENDDO ENDDO - TMP = udvr(nf)%V + TMP = udvr(nf_eff)%V endif - CALL UDV_WRAP_Pivot(TMP1(:,1:UDVR(nf)%N_part), udvr(nf)%U, udvr(nf)%D, V1,NCON,Ndim,UDVR(nf)%N_part) - if(allocated(udvr(nf)%V)) CALL MMULT(udvr(nf)%V, V1, TMP) + CALL UDV_WRAP_Pivot(TMP1(:,1:UDVR(nf_eff)%N_part), udvr(nf_eff)%U, udvr(nf_eff)%D, V1,NCON,Ndim,UDVR(nf_eff)%N_part) + if(allocated(udvr(nf_eff)%V)) CALL MMULT(udvr(nf_eff)%V, V1, TMP) ENDDO #else Implicit None @@ -104,17 +105,18 @@ module wrapur_mod ! Working space. - Integer :: NT, n, nf + Integer :: NT, n, nf, nf_eff - Do nf = 1,N_FL + Do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) DO NT = NTAU + 1, NTAU1 - Call Hop_mod_mmthR(UDVR(nf)%U,nf) + Call Hop_mod_mmthR(UDVR(nf_eff)%U,nf) Do n = 1,Size(Op_V,1) - Call Op_mmultR(UDVR(nf)%U,Op_V(n,nf),nsigma%f(n,nt),'n') + Call Op_mmultR(UDVR(nf_eff)%U,Op_V(n,nf),nsigma%f(n,nt),'n') ENDDO ENDDO - CALL UDVR(nf)%decompose + CALL UDVR(nf_eff)%decompose ENDDO #endif diff --git a/testsuite/Prog.tests/9-Op-Phase.F90 b/testsuite/Prog.tests/9-Op-Phase.F90 index a02f8f2b6589c0b8c084dacc9961e39d71d247b0..b1fa57947e2ad9b46fd8f9a8e01cb841499c02b0 100644 --- a/testsuite/Prog.tests/9-Op-Phase.F90 +++ b/testsuite/Prog.tests/9-Op-Phase.F90 @@ -34,7 +34,10 @@ Program TESTOP_PHASE Phasenew = 1.D0 Phaseold = 1.D0 - Call Op_Phase(Phasenew, Op, NSigma, N_SUN) + Do nf = 1, 3 + Call Op_Phase(Phasenew, Op, NSigma, nf) + enddo + Phasenew=Phasenew**N_SUN ! check against old version do nf = 1, Size(Op,2)