From d863c98d75ab2b043f8033e33af91266fa14941a Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Thu, 27 Feb 2025 15:47:32 +0100 Subject: [PATCH] Fix bug in Pi_1 renormalization. The 'L_n' scheme is still not working yet but 'G_n', 'Sigma_n' and 'Pi_1' are working. --- src/msspec/parameters.py | 14 +- .../spec/fortran/common_sub/read_data.f | 11 +- .../fortran/renormalization/renormalization.f | 321 ++++++++---------- 3 files changed, 164 insertions(+), 182 deletions(-) diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index f4fa39a..8f1e8b3 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -19,8 +19,8 @@ # along with this msspec. If not, see . # # Source file : src/msspec/parameters.py -# Last modified: Tue, 15 Feb 2022 15:37:28 +0100 -# Committed by : Sylvain Tricot +# Last modified: Thu, 27 Feb 2025 15:47:32 +0100 +# Committed by : Sylvain Tricot """ @@ -400,7 +400,7 @@ class SpecParameters(BaseParameters): fmt='d'), Parameter('calctype_ipol', types=int, limits=(-1, 2), default=0, fmt='d'), - Parameter('calctype_iamp', types=int, limits=(0, 1), default=1, + Parameter('calctype_iamp', types=int, limits=(0, 1), default=0, fmt='d'), Parameter('ped_li', types=str, default='1s'), @@ -1446,14 +1446,14 @@ class CalculationParameters(BaseParameters): The scattering order. Only meaningful for the 'expansion' algorithm. Its value is limited to 10."""), Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n', - 'Z_n', 'Pi_1', 'Lowdin'), + 'Z_n', 'Pi_1', 'L_n'), types=(type(None), str), default=None, doc=""" Enable the calculation of the coefficients for the renormalization of the multiple scattering series. You can choose to renormalize in terms of the :math:`G_n`, the - :math:`\\Sigma_n`, the :math:`Z_n`, the :math:`\\Pi_1` or the Lowdin - :math:`K^2` matrices"""), + :math:`\\Sigma_n`, the :math:`Z_n`, the Löwdin :math:`\\Pi_1` or + :math:`L_n` matrices"""), Parameter('renormalization_omega', types=(int,float,complex), default=1.+0j, doc=""" @@ -1589,7 +1589,7 @@ class CalculationParameters(BaseParameters): 'Sigma_n': 2, 'Z_n' : 3, 'Pi_1' : 4, - 'Lowdin' : 5} + 'L_n' : 5} # Check that the method is neither 'Z_n' nor 'K^2' for other # 'spetroscopy' than EIG try: diff --git a/src/msspec/spec/fortran/common_sub/read_data.f b/src/msspec/spec/fortran/common_sub/read_data.f index 0145a27..4354ebd 100644 --- a/src/msspec/spec/fortran/common_sub/read_data.f +++ b/src/msspec/spec/fortran/common_sub/read_data.f @@ -1404,11 +1404,12 @@ C C C Computing the renormalization coefficients C - IF(I_REN.LE.4) THEN - CALL COEF_RENORM(NDIF) - ELSEIF(I_REN.EQ.5) THEN - CALL COEF_LOEWDIN(NDIF) - ENDIF +C IF(I_REN.LE.4) THEN +C CALL COEF_RENORM(NDIF) +C ELSEIF(I_REN.EQ.5) THEN +C CALL COEF_LOEWDIN(NDIF) +C ENDIF + CALL COEF_RENORM(NDIF) C C Storage of the logarithm of the Gamma function GLD(N+1,N_INT) C for integer (N_INT=1) and semi-integer (N_INT=2) values : diff --git a/src/msspec/spec/fortran/renormalization/renormalization.f b/src/msspec/spec/fortran/renormalization/renormalization.f index ef0dde6..1609d1c 100644 --- a/src/msspec/spec/fortran/renormalization/renormalization.f +++ b/src/msspec/spec/fortran/renormalization/renormalization.f @@ -1,53 +1,66 @@ SUBROUTINE COEF_RENORM(NDIF) -C -C This subroutine computes the coefficients for the renormalization -C of the multiple scattering series. These coefficients are -C expressed as C_REN(K) where K is the multiple scattering order. -C REN2 is the value of the mixing (or renormalization) parameter. -C -C NDIF is the scattering order at which the series is truncated, -C so that K varies from 0 to NDIF. -C -C COMMON /RENORM/: -C -C I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN) -C = 2 : renormalization in terms of the Sigma_n matrices -C = 3 : renormalization in terms of the Z_n matrices -C = 4 : renormalization in terms of the Pi_1 matrix -C -C N_REN = n -C -C REN = REN_R+IC*REN_I : omega -C -C Last modified : 11 Apr 2019 -C +! +! This subroutine computes the coefficients for the renormalization +! of the multiple scattering series. These coefficients are +! expressed as C_REN(K) where K is the multiple scattering order. +! REN2 is the value of the mixing (or renormalization) parameter. +! +! NDIF is the scattering order at which the series is truncated, +! so that K varies from 0 to NDIF. +! +! COMMON /RENORM/: +! +! I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN) +! = 2 : renormalization in terms of the Sigma_n matrices +! = 3 : renormalization in terms of the Z_n matrices +! = 4 : Löwdin renormalization in terms of the Pi_1 matrices +! = 5 : Löwdin renormalization in terms of the L_n matrices +! +! N_REN = renormalization order n +! +! REN = REN_R + i * REN_I : omega +! +! +! Reference: A. Takatsu, S. Tricot, P. Schieffer, K. Dunseath, +! M. Terao-Dunseath, K. Hatada and D. Sébilleau, +! Phys. Chem. Chem. Phys., 2022, 24, 5658 +! +! +! Authors : D. Sébilleau, A. Takatsu, M. Terao-Dunseath, K. Dunseath +! +! +! Last modified (DS,ST): 26 Feb 2025 +! USE DIM_MOD USE C_RENORM_MOD USE RENORM_MOD -C +! + INTEGER M_MIN, M_MAX +! REAL C(0:NDIF_M,0:NDIF_M) -C +! + COMPLEX X(0:NDIF_M,0:NDIF_M),SUM_L,POWER COMPLEX REN,REN2,COEF1,COEF2,ZEROC,ONEC,IC COMPLEX Y1(0:NDIF_M,0:NDIF_M) -C -C +! +! ZEROC=(0.,0.) ONEC=(1.,0.) IC=(0.,1.) -C +! REN=REN_R+IC*REN_I ! omega -C -C Initialisation of renormalization coefficients -C +! +! Initialisation of renormalization coefficients +! DO J=0,NDIF C_REN(J)=ZEROC ENDDO -C -C Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)! -C (K) -CCCC 2019.06.09 Aika +! +! Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)! +! (K) +!CCC 2019.06.09 Aika c=0.0 -CCCC 2019.06.09 Aika +!CCC 2019.06.09 Aika C(0,0)=1. C(1,0)=1. C(1,1)=1. @@ -58,153 +71,121 @@ CCCC 2019.06.09 Aika C(N,K)=C(N-1,K)+C(N-1,K-1) ENDDO ENDDO -C +! IF(I_REN.LE.3) THEN -C -C Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n) -C +! +! Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n) +! IF(I_REN.EQ.1) THEN -C -C.....(g_n,G_n) renormalization -C +! +!.....(g_n,G_n) renormalization +! REN2=REN**N_REN ! g_n = omega^n -C +! ELSEIF(I_REN.EQ.2) THEN -C -C.....(s_{n},Sigma_n) renormalization -C +! +!.....(s_{n},Sigma_n) renormalization +! REN2=(ONEC-REN**(N_REN+1))/(FLOAT(N_REN+1)*(ONEC-REN)) ! s_n -C +! ELSEIF(I_REN.EQ.3) THEN -C -C.....(zeta_{n},Z_n) renormalization -C -C 2019.04.29 -C REN2=(REN-ONEC)**(N_REN+1) ! zeta_n -C 2019.06.09 -C REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n +! +!.....(zeta_{n},Z_n) renormalization +! +! 2019.04.29 +! REN2=(REN-ONEC)**(N_REN+1) ! zeta_n +! 2019.06.09 +! REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n REN2=-(ONCE-REN)**(N_REN+1) ! zeta_n -C +! ENDIF -C -C AT & MTD 2019.04.17 - summation over j ? +! +! AT & MTD 2019.04.17 - summation over j ? DO K=0,NDIF - c_ren(k)=zeroc + C_REN(K)=ZEROC DO J=K,NDIF - C_REN(K)=C_REN(K)+c(j,k)*(ONEC-REN2)**(J-K) + C_REN(K)=C_REN(K)+C(J,K)*(ONEC-REN2)**(J-K) ENDDO - c_ren(k)=c_ren(k)*ren2**(k+1) + C_REN(K)=C_REN(K)*REN2**(K+1) ENDDO -C -C DO K=0,NDIF -C COEF1=REN2**(K+1) -C DO J=K,NDIF -C COEF2=(ONEC-REN2)**(J-K) -C C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K) -C ENDDO -C ENDDO -C +! +! DO K=0,NDIF +! COEF1=REN2**(K+1) +! DO J=K,NDIF +! COEF2=(ONEC-REN2)**(J-K) +! C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K) +! ENDDO +! ENDDO +! ELSEIF(I_REN.EQ.4) THEN -C -C Loewdin (Pi_1) renormalization for n = 1 -C -C Notation: Y1(M,K) : [Y_1^m]_k -C 2019.06.09 -C Notation: Y1(K,M) : [Y_1^k]_m -C - COEF1=ONEC-REN ! (1 - omega) - DO K=0,NDIF - Y1(K,K)=COEF1**K - IF(K.LT.NDIF) THEN - DO M=K+1,NDIF - COEF2=(REN**(M-K))*(COEF1**(2*K-M)) -C 2019.04.19 AT & MTD -C Y1(M,K)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1)) - Y1(K,M)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1)) - ENDDO - ENDIF - ENDDO -C - DO K=0,NDIF - IN=INT(K/2) - C_REN(K)=ZEROC - DO M=IN,K - C_REN(K)=C_REN(K)+Y1(M,K) - ENDDO - ENDDO -C - ENDIF -C - END -C -C======================================================================= -C - SUBROUTINE COEF_LOEWDIN(NDIF) -C -C This subroutine computes the coefficients for the Loewdin expansion -C of the multiple scattering series. These coefficients are -C expressed as C_LOW(K) where K is the multiple scattering order. -C REN is the value of the mixing (or renormalization) parameter. -C NDIF is the scattering order at which the series is truncated, -C so that K varies from 0 to NDIF. -C -C Corresponds to parameter I_REN = 5 -C -C Notation: X(K,N) = X_n(omega,k) -C -C -C Last modified : 11 Apr 2019 -C - USE DIM_MOD - USE C_RENORM_MOD, C_LOW => C_REN - USE RENORM_MOD -C - COMPLEX X(0:NDIF_M,0:NDIF_M),SUM_L,POWER - COMPLEX REN,ZEROC,ONEC,IC -C -C - ZEROC=(0.,0.) - ONEC=(1.,0.) - IC=(0.,1.) -C - REN=REN_R+IC*REN_I ! omega -C -C Initialisation of renormalization coefficients -C - DO J=0,NDIF - C_LOW(J)=ZEROC - ENDDO -C -C Computing the X(N,K) coefficients, with K <= N -C - POWER=ONEC/REN - DO N=0,NDIF - POWER=POWER*REN ! omega^n - IF(N.EQ.0) THEN - X(N,0)=ONEC - ELSE - X(N,0)=ZEROC - ENDIF - DO K=1,NDIF - IF(K.GT.N) THEN - X(N,K)=ZEROC - ELSEIF(K.EQ.N) THEN - X(N,K)=POWER*X(N-1,K-1) - ELSE - X(N,K)=X(N-1,K)*(REN-POWER) + POWER*X(N-1,K-1) - ENDIF - ENDDO - ENDDO -C -C Calculation of L_n(omega,NDIF) -C - DO N=0,NDIF - SUM_L=ZEROC - DO K=N,NDIF - SUM_L=SUM_L+X(K,N) - ENDDO - C_LOW(N)=SUM_L - ENDDO +! +! Loewdin (Pi_1) renormalization for n = 1 +! +! Notation: Y1(K,M) : [Y_1^k]_m +! +! with k : scattering order +! m : summation index +! + COEF1 = ONEC - REN ! (1 - omega) +! + Y1(0,0) = ONEC ! +! + DO K = 1, NDIF ! + M_MAX = MIN(K,NDIF) ! + M_MIN = INT(K / 2) ! + DO M = M_MIN, M_MAX ! + COEF2 = (REN**(K-M)) * (COEF1**(2*M-K)) ! + Y1(K,M) = COEF2 * ( C(M,K-M) + COEF1 * C(M,K-M-1) ) ! + END DO ! + END DO +! ! + C_REN(0) = ONEC ! + C_REN(1) = ONEC ! + DO K = 2, NDIF ! + C_REN(K) = ZEROC ! + DO M = M_MIN, M_MAX ! + C_REN(K) = C_REN(K) + Y1(M,K) ! + END DO ! + END DO ! + ELSE IF(I_REN.EQ.5) THEN +! +! Loewdin L_n(omega,NDIF) renormalization +! +! Notation: X(K,N) = X_n(omega,k) +! +! +! Computing the X(N,K) coefficients, with K <= N +! + POWER = ONEC / REN ! + DO N = 0, NDIF ! + POWER = POWER * REN ! omega^n + IF(N == 0) THEN ! + X(N,0) = ONEC ! + ELSE ! + X(N,0) = ZEROC ! + END IF ! + DO K = 1, NDIF ! + IF(K > N) THEN ! + X(N,K) = ZEROC ! + ELSE IF(K == N) THEN ! + X(N,K) = POWER * X(N-1,K-1) ! + ELSE ! + X(N,K) = X(N-1,K) * (REN - POWER) + POWER * X(N-1,K-1)! + END IF ! + END DO ! + END DO ! +! +! Calculation of L_n(omega,NDIF) +! + DO N = 0, NDIF ! + SUM_L = ZEROC ! + DO K = N, NDIF ! + SUM_L = SUM_L + X(K,N) ! + END DO ! + C_REN(N) = SUM_L ! + END DO ! +! + END IF ! +! END -