epsi
/
msspec_python3
3
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

Add Renormalization option in the python interface

This commit is contained in:
Sylvain Tricot 2019-11-27 14:09:55 +01:00
parent 13a43b0b70
commit 4a5f6f1161
9 changed files with 354 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Jenkinsfile vendored
View File

@ -13,7 +13,7 @@ pipeline {
sh '/bin/bash ./CI/CI.bash -t ci_venv'
}
}
stage('Creatin a setup file and test installation...') {
stage('Creating a setup file and test installation...') {
steps {
sh '/bin/bash ./CI/CI.bash -p ci_venv'
sh '/bin/bash ./package/MsSpec*.setup --accept -- -p ./ci_venv -y'

8
src/msspec/calcio.py
View File

@ -637,6 +637,12 @@ class SpecIO(object):
line = fillstr(line, str(p.get_parameter('calc_ispher')), 29, 'left')
line = fillstr(line, str(p.get_parameter('calc_igr')), 39, 'left')
content += line
line = create_line("I_REN,N_REN,REN_R,REN_I")
line = fillstr(line, str(p.get_parameter('calc_iren')), 9, 'left')
line = fillstr(line, str(p.get_parameter('calc_nren')), 19, 'left')
line = fillstr(line, str(p.get_parameter('calc_renr')), 29, 'decimal')
line = fillstr(line, str(p.get_parameter('calc_reni')), 39, 'decimal')
content += line
line = create_line("ISFLIP,IR_DIA,ITRTL,I_TEST")
line = fillstr(line, str(p.get_parameter('calc_isflip')), 9, 'left')
line = fillstr(line, str(p.get_parameter('calc_irdia')), 19, 'left')
@ -856,7 +862,7 @@ class SpecIO(object):
'N_CL_L_M': 0,
'NE_M': self.phagenio.ne,
'NL_M': self.phagenio.nlmax + 1,
'LI_M': get_li(self.parameters.extra_level),
'LI_M': get_li(self.parameters.extra_level) + 1,
'NEMET_M': 1,
'NO_ST_M': self.parameters.calc_no,
}

15
src/msspec/calculator.py
View File

@ -165,6 +165,19 @@ class _MSCALCULATOR(Calculator):
self.calculation_parameters = CalculationParameters(
self.global_parameters, self.phagen_parameters, self.spec_parameters)
# initialize all parameters with defaults values
LOGGER.info("Set default values =========================================")
for p in (list(self.global_parameters) +
list(self.muffintin_parameters) +
list(self.tmatrix_parameters) +
list(self.source_parameters) +
list(self.detector_parameters) +
list(self.scan_parameters) +
list(self.calculation_parameters) +
list(self.spectroscopy_parameters)):
p.set(p.default)
LOGGER.info("End of default values ======================================")
# updated global parameters with provided keywords
self.global_parameters.spectroscopy = spectroscopy
self.global_parameters.algorithm = algorithm
@ -344,7 +357,7 @@ class _MSCALCULATOR(Calculator):
'N_CL_L_M' : 1,
'NE_M' : self.phagenio.ne,
'NL_M' : self.phagenio.nlmax + 1,
'LI_M' : get_li(self.spec_parameters.extra_level),
'LI_M' : get_li(self.spec_parameters.extra_level) + 1,
'NEMET_M' : 1,
'NO_ST_M' : self.spec_parameters.calc_no,
'NDIF_M' : 10,

90
src/msspec/parameters.py
View File

@ -1,4 +1,5 @@
# coding: utf-8
# vim: set et ts=4 sw=4 sts nu fdm=indent:
"""
Module parameters
=================
@ -573,6 +574,12 @@ class SpecParameters(BaseParameters):
Parameter('calc_ispher', types=int, limits=[0, 1], default=1,
fmt='d'),
Parameter('calc_igr', types=int, limits=[0, 2], default=0, fmt='d'),
Parameter('calc_iren', types=int, limits=[0, 4], default=1, fmt='d'),
Parameter('calc_nren', types=int, limits=[1, None], default=1, fmt='d'),
Parameter('calc_renr', types=float, limits=[None, None], default=1., fmt='.3f'),
Parameter('calc_reni', types=float, limits=[None, None], default=0., fmt='.3f'),
Parameter('calc_isflip', types=int, limits=[0, 1], default=0,
fmt='d'),
Parameter('calc_irdia', types=int, limits=[0, 1], default=0,
@ -1233,14 +1240,14 @@ class ScanParameters(BaseParameters):
spectro = self.global_parameters.spectroscopy
scantype = self.get_parameter('type').value
if scantype == 'scatf':
comment = 'with the scattering factor scan type.'
if spectro == 'EXAFS':
comment = 'in EXAFS spetroscopy.'
if spectro in ('EXAFS',) or scantype in ('scatf',):
msg = 'Setting the theta angle is not possible %s' % comment
LOGGER.error('Incompatible options!')
raise ValueError(msg)
#if scantype == 'scatf':
# comment = 'with the scattering factor scan type.'
#if spectro == 'EXAFS':
# comment = 'in EXAFS spetroscopy.'
#if spectro in ('EXAFS',) or scantype in ('scatf',):
# msg = 'Setting the theta angle is not possible %s' % comment
# LOGGER.error('Incompatible options!')
# raise ValueError(msg)
# p._value = np.array(p.value, dtype=np.float).flatten()
arr = np.array(p.value, dtype=np.float).flatten()
@ -1268,14 +1275,14 @@ class ScanParameters(BaseParameters):
spectro = self.global_parameters.spectroscopy
scantype = self.get_parameter('type').value
if scantype == 'scatf':
comment = 'with scattering factor scan type.'
if spectro == 'EXAFS':
comment = 'in EXAFS spetroscopy.'
if spectro in ('EXAFS',) or scantype in ('scatf',):
msg = 'Setting the phi angle is not possible %s' % comment
LOGGER.error('Incompatible options')
raise ValueError(msg)
#if scantype == 'scatf':
# comment = 'with scattering factor scan type.'
#if spectro == 'EXAFS':
# comment = 'in EXAFS spetroscopy.'
#if spectro in ('EXAFS',) or scantype in ('scatf',):
# msg = 'Setting the phi angle is not possible %s' % comment
# LOGGER.error('Incompatible options')
# raise ValueError(msg)
arr = np.array(p.value, dtype=np.float).flatten()
@ -1365,6 +1372,18 @@ class CalculationParameters(BaseParameters):
doc=textwrap.dedent("""
The scattering order. Only meaningful for the 'expansion' algorithm.
Its value is limited to 10.""")),
Parameter('renormalization_mode', allowed_values=(None, 'Sigma_n', 'G_n'),
types=(type(None), str), default=None,
doc=textwrap.dedent("""
Enable the calculation of the coefficients for the renormalization of
the multiple scattering series.
Only meaningful for the 'expansion' algorithm. You can choose to
renormalize in terms of the Sigma_n or G_n matrices.""")),
Parameter('renormalization_omega', types=(int,float,complex),
default=1.+0j,
doc=textwrap.dedent("""
The :math:`\\omega` coefficient used to initialize the
renormalization alogorithm.""")),
Parameter('RA_cutoff_damping', types=int, limits=(0, None),
default=0, doc=textwrap.dedent("""
The Rehr-Albers truncation order. If > 0, the *RA_cutoff* is
@ -1489,6 +1508,22 @@ class CalculationParameters(BaseParameters):
self.spec_parameters.calc_ndif = p.value
LOGGER.info('Scattering order set to %s', p.value)
def bind_renormalization_mode(self, p):
if p.value is None:
self.spec_parameters.calc_iren = 0
else:
if p.value.lower() == 'Sigma_n'.lower():
self.spec_parameters.calc_iren = 1
elif p.value.lower() == 'G_n'.lower():
self.spec_parameters.calc_iren = 2
LOGGER.info(f"Renormalization activated with \'{p.value}\' method")
def bind_renormalization_omega(self, p):
omega = complex(p.value)
self.spec_parameters.calc_renr = omega.real
self.spec_parameters.calc_reni = omega.imag
LOGGER.info(f"Renormalization omega set to \'{p.value}\'")
def bind_RA_cutoff_damping(self, p):
self.spec_parameters.calc_ino = p.value
LOGGER.info('Rehr-Albers cutoff damping set to %s', p.value)
@ -1502,7 +1537,8 @@ class CalculationParameters(BaseParameters):
LOGGER.info('Type of basis functions: \'%s\'', p.value)
def bind_spin_flip(self, p):
if self.global_parameters.spinpol is False:
isflip = int(p.value)
if self.global_parameters.spinpol is False and isflip:
err_msg = (
"'{}' is ignored since the 'spinpol' global parameter is set "
"to False. Enable spin polarization in the constructor of "
@ -1510,26 +1546,24 @@ class CalculationParameters(BaseParameters):
).format(p.name)
LOGGER.error("Incompatible options!")
raise ValueError(err_msg)
isflip = int(p.value)
self.spec_parameters.calc_isflip = isflip
LOGGER.info('Spin-flip set to %s', p.value)
def bind_integrals(self, p):
if self.global_parameters.spinpol is False:
err_msg = (
"'{}' is ignored since the 'spinpol' global parameter is set "
"to False. Enable spin polarization in the constructor of "
"your Calculator if you want to use this option."
).format(p.name)
LOGGER.error("Incompatible options!")
raise ValueError(err_msg)
irdia = 0 if p.value == 'all' else 1
self.spec_parameters.calc_irdia = irdia
LOGGER.info('Radial integrals taken into account: %s', p.value)
if self.global_parameters.spinpol is False:
LOGGER.warning(
f"'{p.name}' is ignored since the 'spinpol' global parameter is set "
"to False. Enable spin polarization in the constructor of "
"your Calculator if you want to use this option.")
def bind_path_filtering(self, p):
ifwd = ipw = ilength = 0
ipp = 1
if p.value != None:
if ('plane_wave_spin_averaged' in p.value and 'plane_wave_normal' in
p.value):
err_msg = (
@ -1538,8 +1572,6 @@ class CalculationParameters(BaseParameters):
"with other filters if needed).")
LOGGER.error('Incompatible options!')
raise ValueError(err_msg)
if p.value != None:
if 'forward_scattering' in p.value:
ifwd = 1
if 'backward_scattering' in p.value:
@ -1702,7 +1734,7 @@ class PEDParameters(BaseParameters):
None: 0,
'single': 1,
'both': 2}
self.spec_parameters.ped_so = somap[p.value]
self.spec_parameters.ped_iso = somap[p.value]
class EIGParameters(BaseParameters):

2
src/msspec/spec/fortran/Makefile
View File

@ -1,6 +1,6 @@
COMP=gfortran
objects_src := dim_mod.f modules.f allocation.f spec.f
objects_src := dim_mod.f modules.f renormalization.f allocation.f spec.f
objects := $(patsubst %.f,%.o, $(objects_src))
OPTS := -g -Wall -Wextra -Warray-temporaries -Wconversion -fbacktrace -ffree-line-length-0 -fcheck=all -ffpe-trap=zero,overflow,underflow -finit-real=nan

4
src/msspec/spec/fortran/allocation.f
View File

@ -10,6 +10,7 @@
USE BASES_MOD
USE CLUSLIM_MOD
USE COOR_MOD
USE C_RENORM_MOD
USE DEBWAL_MOD
USE INDAT_MOD
USE INIT_A_MOD
@ -27,6 +28,7 @@
USE PARCAL_A_MOD
USE RELADS_MOD
USE RELAX_MOD
USE RENORM_MOD
USE RESEAU_MOD
USE SPIN_MOD
USE TESTS_MOD
@ -113,6 +115,7 @@
CALL ALLOC_BASES()
CALL ALLOC_CLUSLIM()
CALL ALLOC_COOR()
CALL ALLOC_C_RENORM()
CALL ALLOC_DEBWAL()
CALL ALLOC_INDAT()
CALL ALLOC_INIT_A()
@ -130,6 +133,7 @@
CALL ALLOC_PARCAL_A()
CALL ALLOC_RELADS()
CALL ALLOC_RELAX()
CALL ALLOC_RENORM()
CALL ALLOC_RESEAU()
CALL ALLOC_SPIN()
CALL ALLOC_TESTS()

24
src/msspec/spec/fortran/modules.f
View File

@ -1416,4 +1416,26 @@ C=======================================================================
ALLOCATE(JPON(NPATH_M,NDIF_M))
END SUBROUTINE ALLOC_PRINTP
END MODULE PRINTP_MOD
C=======================================================================
MODULE RENORM_MOD
IMPLICIT NONE
INTEGER :: I_REN, N_REN
REAL :: REN_R, REN_I
CONTAINS
SUBROUTINE ALLOC_RENORM()
USE DIM_MOD
END SUBROUTINE ALLOC_RENORM
END MODULE RENORM_MOD
C=======================================================================
MODULE C_RENORM_MOD
IMPLICIT NONE
COMPLEX, ALLOCATABLE, DIMENSION(:) :: C_REN
CONTAINS
SUBROUTINE ALLOC_C_RENORM()
USE DIM_MOD
IF (ALLOCATED(C_REN)) THEN
DEALLOCATE(C_REN)
ENDIF
ALLOCATE(C_REN(0:NDIF_M))
END SUBROUTINE ALLOC_C_RENORM
END MODULE C_RENORM_MOD

210
src/msspec/spec/fortran/renormalization.f Normal file
View File

@ -0,0 +1,210 @@
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
USE DIM_MOD
USE C_RENORM_MOD
USE RENORM_MOD
C
REAL C(0:NDIF_M,0:NDIF_M)
C
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
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
c=0.0
CCCC 2019.06.09 Aika
C(0,0)=1.
C(1,0)=1.
C(1,1)=1.
DO N=2,NDIF
C(N,0)=1.
C(N,N)=1.
DO K=1,N-1
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
IF(I_REN.EQ.1) THEN
C
C.....(g_n,G_n) renormalization
C
REN2=REN**N_REN ! g_n = omega^n
C
ELSEIF(I_REN.EQ.2) THEN
C
C.....(s_{n},Sigma_n) renormalization
C
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
REN2=-(ONCE-REN)**(N_REN+1) ! zeta_n
C
ENDIF
C
C AT & MTD 2019.04.17 - summation over j ?
DO K=0,NDIF
c_ren(k)=zeroc
DO J=K,NDIF
C_REN(K)=C_REN(K)+c(j,k)*(ONEC-REN2)**(J-K)
ENDDO
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
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
END

25
src/msspec/spec/fortran/spec.f
View File

@ -4976,6 +4976,7 @@ C
USE RA_MOD
USE RELADS_MOD
USE RELAX_MOD
USE RENORM_MOD
USE RESEAU_MOD
USE SPECTRUM_MOD
USE SPIN_MOD
@ -5601,6 +5602,7 @@ C
READ(ICOM,1) RIEN
C
READ(ICOM,21) NO,NDIF,ISPHER,I_GR
READ(ICOM,50) I_REN,N_REN,REN_R,REN_I
C
IF(ISPHER.EQ.0) THEN
IDWSPH=0
@ -6211,6 +6213,7 @@ C
IF(SPECTRO.NE.'EIG') THEN
C
WRITE(IUO1,121) NO,NDIF,ISPHER,I_GR
WRITE(IUO1,150) I_REN,N_REN,REN_R,REN_I
C
IF(SPECTRO.EQ.'XAS') NDIF=NDIF+1
C
@ -6316,6 +6319,14 @@ C
ENDIF
ENDIF
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
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
@ -6667,6 +6678,7 @@ C
47 FORMAT(5X,I5,6X,I4,6X,I4,8X,F6.3)
48 FORMAT(9X,I1,9X,I1,9X,I1,9X,I1)
49 FORMAT(8X,I2,6X,F7.2,3X,F7.2)
50 FORMAT(9X,I1,9X,I1,6X,F8.3,2X,F8.3)
C
C
C.................... Write FORMAT ....................
@ -6733,6 +6745,7 @@ C
147 FORMAT(8X,I5,6X,I4,6X,I4,8X,F6.3,5X,'N_MAX,N_ITER,N_TABLE,SHIFT')
148 FORMAT(12X,I1,9X,I1,9X,I1,9X,I1,9X,'I_XN,I_VA,I_GN,I_WN')
149 FORMAT(11X,I2,6X,F7.2,3X,F7.2,16X,'L,ALPHA,BETA')
150 FORMAT(12X,I1,9X,I1,6X,F8.3,2X,F8.3,5X,'I_REN,N_REN,REN_R,REN_I')
C
201 FORMAT(///,21X,10A4,////)
203 FORMAT('**************************************************',
@ -10575,6 +10588,7 @@ C
USE DIM_MOD
C
USE APPROX_MOD
USE C_RENORM_MOD
USE EXPFAC_MOD
USE EXTREM_MOD
USE INIT_L_MOD
@ -10585,6 +10599,7 @@ C
USE PATH_MOD
USE PRINTP_MOD
USE RA_MOD
USE RENORM_MOD
USE ROT_MOD
USE SCATMAT_MOD , F => F21
USE TESTS_MOD
@ -10672,6 +10687,13 @@ C
COEF=COEF*CEXDW(JSC)
ENDDO
C
C Renormalization of the path
C
IF(I_REN.GE.1) THEN
COEF=COEF*C_REN(JORDP)
write(354,*) JORDP,C_REN(JORDP)
ENDIF
C
C Call of the subroutines used for the R-A termination matrix
C This termination matrix is now merged into PATHOP
C
@ -11245,6 +11267,7 @@ C
USE AMPLI_MOD
USE APPROX_MOD
USE COOR_MOD , NTCLU => NATCLU, NTP => NATYP
USE C_RENORM_MOD
USE DEBWAL_MOD
USE DIRECT_MOD , RTHETA => RTHEXT
USE EXTREM_MOD
@ -11261,6 +11284,7 @@ C
USE PARCAL_MOD
USE PATH_MOD
USE PRINTP_MOD
USE RENORM_MOD
USE RESEAU_MOD
USE SPIN_MOD
USE TESTPA_MOD
@ -12135,6 +12159,7 @@ C
ELSE
R2=TLT(LF,1,1,JE)
ENDIF
IF(I_REN.GE.1) R2=R2*C_REN(0)
DO MF=-LF,LF
MR=2+MF-MI
LMR=LRR+MR