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.
This commit is contained in:
Sylvain Tricot 2025-02-27 15:47:32 +01:00
parent 5c5c57be2b
commit d863c98d75
3 changed files with 164 additions and 182 deletions

View File

@ -19,8 +19,8 @@
# along with this msspec. If not, see <http://www.gnu.org/licenses/>. # along with this msspec. If not, see <http://www.gnu.org/licenses/>.
# #
# Source file : src/msspec/parameters.py # Source file : src/msspec/parameters.py
# Last modified: Tue, 15 Feb 2022 15:37:28 +0100 # Last modified: Thu, 27 Feb 2025 15:47:32 +0100
# Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes1.fr> # Committed by : Sylvain Tricot <sylvain.tricot@univ-rennes.fr>
""" """
@ -400,7 +400,7 @@ class SpecParameters(BaseParameters):
fmt='d'), fmt='d'),
Parameter('calctype_ipol', types=int, limits=(-1, 2), default=0, Parameter('calctype_ipol', types=int, limits=(-1, 2), default=0,
fmt='d'), fmt='d'),
Parameter('calctype_iamp', types=int, limits=(0, 1), default=1, Parameter('calctype_iamp', types=int, limits=(0, 1), default=0,
fmt='d'), fmt='d'),
Parameter('ped_li', types=str, default='1s'), Parameter('ped_li', types=str, default='1s'),
@ -1446,14 +1446,14 @@ class CalculationParameters(BaseParameters):
The scattering order. Only meaningful for the 'expansion' algorithm. The scattering order. Only meaningful for the 'expansion' algorithm.
Its value is limited to 10."""), Its value is limited to 10."""),
Parameter('renormalization_mode', allowed_values=(None, 'G_n', 'Sigma_n', 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, types=(type(None), str), default=None,
doc=""" doc="""
Enable the calculation of the coefficients for the renormalization of Enable the calculation of the coefficients for the renormalization of
the multiple scattering series. the multiple scattering series.
You can choose to renormalize in terms of the :math:`G_n`, the 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:`\\Sigma_n`, the :math:`Z_n`, the Löwdin :math:`\\Pi_1` or
:math:`K^2` matrices"""), :math:`L_n` matrices"""),
Parameter('renormalization_omega', types=(int,float,complex), Parameter('renormalization_omega', types=(int,float,complex),
default=1.+0j, default=1.+0j,
doc=""" doc="""
@ -1589,7 +1589,7 @@ class CalculationParameters(BaseParameters):
'Sigma_n': 2, 'Sigma_n': 2,
'Z_n' : 3, 'Z_n' : 3,
'Pi_1' : 4, 'Pi_1' : 4,
'Lowdin' : 5} 'L_n' : 5}
# Check that the method is neither 'Z_n' nor 'K^2' for other # Check that the method is neither 'Z_n' nor 'K^2' for other
# 'spetroscopy' than EIG # 'spetroscopy' than EIG
try: try:

View File

@ -1404,11 +1404,12 @@ C
C C
C Computing the renormalization coefficients C Computing the renormalization coefficients
C C
IF(I_REN.LE.4) THEN C IF(I_REN.LE.4) THEN
CALL COEF_RENORM(NDIF) C CALL COEF_RENORM(NDIF)
ELSEIF(I_REN.EQ.5) THEN C ELSEIF(I_REN.EQ.5) THEN
CALL COEF_LOEWDIN(NDIF) C CALL COEF_LOEWDIN(NDIF)
ENDIF C ENDIF
CALL COEF_RENORM(NDIF)
C C
C Storage of the logarithm of the Gamma function GLD(N+1,N_INT) 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 : C for integer (N_INT=1) and semi-integer (N_INT=2) values :

View File

@ -1,53 +1,66 @@
SUBROUTINE COEF_RENORM(NDIF) SUBROUTINE COEF_RENORM(NDIF)
C !
C This subroutine computes the coefficients for the renormalization ! This subroutine computes the coefficients for the renormalization
C of the multiple scattering series. These coefficients are ! of the multiple scattering series. These coefficients are
C expressed as C_REN(K) where K is the multiple scattering order. ! expressed as C_REN(K) where K is the multiple scattering order.
C REN2 is the value of the mixing (or renormalization) parameter. ! REN2 is the value of the mixing (or renormalization) parameter.
C !
C NDIF is the scattering order at which the series is truncated, ! NDIF is the scattering order at which the series is truncated,
C so that K varies from 0 to NDIF. ! so that K varies from 0 to NDIF.
C !
C COMMON /RENORM/: ! COMMON /RENORM/:
C !
C I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN) ! I_REN = 1 : renormalization in terms of G_n matrices (n : N_REN)
C = 2 : renormalization in terms of the Sigma_n matrices ! = 2 : renormalization in terms of the Sigma_n matrices
C = 3 : renormalization in terms of the Z_n matrices ! = 3 : renormalization in terms of the Z_n matrices
C = 4 : renormalization in terms of the Pi_1 matrix ! = 4 : Löwdin renormalization in terms of the Pi_1 matrices
C ! = 5 : Löwdin renormalization in terms of the L_n matrices
C N_REN = n !
C ! N_REN = renormalization order n
C REN = REN_R+IC*REN_I : omega !
C ! REN = REN_R + i * REN_I : omega
C Last modified : 11 Apr 2019 !
C !
! 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 DIM_MOD
USE C_RENORM_MOD USE C_RENORM_MOD
USE RENORM_MOD USE RENORM_MOD
C !
INTEGER M_MIN, M_MAX
!
REAL C(0:NDIF_M,0:NDIF_M) 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 REN,REN2,COEF1,COEF2,ZEROC,ONEC,IC
COMPLEX Y1(0:NDIF_M,0:NDIF_M) COMPLEX Y1(0:NDIF_M,0:NDIF_M)
C !
C !
ZEROC=(0.,0.) ZEROC=(0.,0.)
ONEC=(1.,0.) ONEC=(1.,0.)
IC=(0.,1.) IC=(0.,1.)
C !
REN=REN_R+IC*REN_I ! omega REN=REN_R+IC*REN_I ! omega
C !
C Initialisation of renormalization coefficients ! Initialisation of renormalization coefficients
C !
DO J=0,NDIF DO J=0,NDIF
C_REN(J)=ZEROC C_REN(J)=ZEROC
ENDDO ENDDO
C !
C Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)! ! Computing the binomial coefficients C(N,K) = (N) = N! / K! (N-K)!
C (K) ! (K)
CCCC 2019.06.09 Aika !CCC 2019.06.09 Aika
c=0.0 c=0.0
CCCC 2019.06.09 Aika !CCC 2019.06.09 Aika
C(0,0)=1. C(0,0)=1.
C(1,0)=1. C(1,0)=1.
C(1,1)=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) C(N,K)=C(N-1,K)+C(N-1,K-1)
ENDDO ENDDO
ENDDO ENDDO
C !
IF(I_REN.LE.3) THEN IF(I_REN.LE.3) THEN
C !
C Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n) ! Computing the modified renormalization parameter REN2 (g_n,s_n,zeta_n)
C !
IF(I_REN.EQ.1) THEN IF(I_REN.EQ.1) THEN
C !
C.....(g_n,G_n) renormalization !.....(g_n,G_n) renormalization
C !
REN2=REN**N_REN ! g_n = omega^n REN2=REN**N_REN ! g_n = omega^n
C !
ELSEIF(I_REN.EQ.2) THEN ELSEIF(I_REN.EQ.2) THEN
C !
C.....(s_{n},Sigma_n) renormalization !.....(s_{n},Sigma_n) renormalization
C !
REN2=(ONEC-REN**(N_REN+1))/(FLOAT(N_REN+1)*(ONEC-REN)) ! s_n REN2=(ONEC-REN**(N_REN+1))/(FLOAT(N_REN+1)*(ONEC-REN)) ! s_n
C !
ELSEIF(I_REN.EQ.3) THEN ELSEIF(I_REN.EQ.3) THEN
C !
C.....(zeta_{n},Z_n) renormalization !.....(zeta_{n},Z_n) renormalization
C !
C 2019.04.29 ! 2019.04.29
C REN2=(REN-ONEC)**(N_REN+1) ! zeta_n ! REN2=(REN-ONEC)**(N_REN+1) ! zeta_n
C 2019.06.09 ! 2019.06.09
C REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n ! REN2=-(REN-ONEC)**(N_REN+1) ! zeta_n
REN2=-(ONCE-REN)**(N_REN+1) ! zeta_n REN2=-(ONCE-REN)**(N_REN+1) ! zeta_n
C !
ENDIF ENDIF
C !
C AT & MTD 2019.04.17 - summation over j ? ! AT & MTD 2019.04.17 - summation over j ?
DO K=0,NDIF DO K=0,NDIF
c_ren(k)=zeroc C_REN(K)=ZEROC
DO J=K,NDIF 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 ENDDO
c_ren(k)=c_ren(k)*ren2**(k+1) C_REN(K)=C_REN(K)*REN2**(K+1)
ENDDO ENDDO
C !
C DO K=0,NDIF ! DO K=0,NDIF
C COEF1=REN2**(K+1) ! COEF1=REN2**(K+1)
C DO J=K,NDIF ! DO J=K,NDIF
C COEF2=(ONEC-REN2)**(J-K) ! COEF2=(ONEC-REN2)**(J-K)
C C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K) ! C_REN(J)=C_REN(J)+COEF1*COEF2*C(J,K)
C ENDDO ! ENDDO
C ENDDO ! ENDDO
C !
ELSEIF(I_REN.EQ.4) THEN ELSEIF(I_REN.EQ.4) THEN
C !
C Loewdin (Pi_1) renormalization for n = 1 ! Loewdin (Pi_1) renormalization for n = 1
C !
C Notation: Y1(M,K) : [Y_1^m]_k ! Notation: Y1(K,M) : [Y_1^k]_m
C 2019.06.09 !
C Notation: Y1(K,M) : [Y_1^k]_m ! with k : scattering order
C ! m : summation index
COEF1=ONEC-REN ! (1 - omega) !
DO K=0,NDIF COEF1 = ONEC - REN ! (1 - omega)
Y1(K,K)=COEF1**K !
IF(K.LT.NDIF) THEN Y1(0,0) = ONEC !
DO M=K+1,NDIF !
COEF2=(REN**(M-K))*(COEF1**(2*K-M)) DO K = 1, NDIF !
C 2019.04.19 AT & MTD M_MAX = MIN(K,NDIF) !
C Y1(M,K)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1)) M_MIN = INT(K / 2) !
Y1(K,M)=COEF2*(C(K,M-K)+COEF1*C(K,M-K-1)) DO M = M_MIN, M_MAX !
ENDDO COEF2 = (REN**(K-M)) * (COEF1**(2*M-K)) !
ENDIF Y1(K,M) = COEF2 * ( C(M,K-M) + COEF1 * C(M,K-M-1) ) !
ENDDO END DO !
C END DO
DO K=0,NDIF ! !
IN=INT(K/2) C_REN(0) = ONEC !
C_REN(K)=ZEROC C_REN(1) = ONEC !
DO M=IN,K DO K = 2, NDIF !
C_REN(K)=C_REN(K)+Y1(M,K) C_REN(K) = ZEROC !
ENDDO DO M = M_MIN, M_MAX !
ENDDO C_REN(K) = C_REN(K) + Y1(M,K) !
C END DO !
ENDIF END DO !
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
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 END