diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 96fa276..d848c56 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -1536,7 +1536,7 @@ class CalculationParameters(BaseParameters): #assert( p.value in (None, 'Sigma_n', 'G_n') ) assert( p.value in p.allowed_values) elif (self.global_parameters.spectroscopy == 'EIG' - and self.global_parameters.algorithm == 'power'): + and self.global_parameters.algorithm == 'inversion'): assert( p.value in p.allowed_values) else: assert( p.value is None ) diff --git a/src/msspec/spec/fortran/eig/mi/new.f.hidden b/src/msspec/spec/fortran/eig/mi/new.f.hidden deleted file mode 100644 index 1c93f94..0000000 --- a/src/msspec/spec/fortran/eig/mi/new.f.hidden +++ /dev/null @@ -1,3972 +0,0 @@ -C -C -C ************************************************************ -C * ******************************************************** * -C * * * * -C * * MULTIPLE-SCATTERING SPIN-INDEPENDENT * * -C * * EIGENVALUE CALCULATION CODE * * -C * * * * -C * ******************************************************** * -C ************************************************************ -C -C -C -C -C Written by D. Sebilleau, Groupe Theorie, -C Departement Materiaux-Nanosciences, -C Institut de Physique de Rennes, -C UMR CNRS-Universite 6251, -C Universite de Rennes-1, -C 35042 Rennes-Cedex, -C France -C -C Contributions : M. Gavaza, H.-F. Zhao, K. Hatada -C -C----------------------------------------------------------------------- -C -C As a general rule in this code, although there might be a few -C exceptions (...), a variable whose name starts with a 'I' is a -C switch, with a 'J' is a loop index and with a 'N' is a number. -C -C The main subroutines are : -C -C * PHDDIF : computes the photoelectron diffraction -C formula -C -C * XASDIF : computes the EXAFS or XANES formula -C depending on the energy -C -C * AEDDIF : computes the Auger electron diffraction -C formula -C -C * FINDPATHS : generates the multiple scattering -C paths the electron will follow -C -C * PATHOP : calculates the contribution of a given -C path to the scattering path operator -C -C * MATDIF : computes the Rehr-Albers scattering -C matrices -C -C A subroutine called NAME_A is the Auger equivalent of subroutine -C NAME. The essentail difference between NAME and NAME_A is that -C they do not contain the same arrays. -C -C Always remember, when changing the input data file, to keep the -C format. The rule here is that the last digit of any integer or -C character data must correspond to the tab (+) while for real data, -C the tab precedes the point. -C -C Do not forget, before submitting a calculation, to check the -C consistency of the input data with the corresponding maximal -C values in the include file. -C -C----------------------------------------------------------------------- -C -C Please report any bug or problem to me at : -C -C didier.sebilleau@univ-rennes1.fr -C -C -C -C Last modified : 10 Jan 2016 -C -C======================================================================= -C - PROGRAM MAIN -C -C This routine reads the various input files and calls the subroutine -C performing the requested calculation -C -C INCLUDE 'spec.inc' -C - USE ADSORB_MOD - USE APPROX_MOD - USE ATOMS_MOD - USE AUGER_MOD - USE BASES_MOD - USE CLUSLIM_MOD - USE COOR_MOD - USE DEBWAL_MOD - USE INDAT_MOD - USE INIT_A_MOD - USE INIT_L_MOD - USE INIT_J_MOD - USE INIT_M_MOD - USE INFILES_MOD - USE INUNITS_MOD - USE LIMAMA_MOD - USE LPMOY_MOD - USE MASSAT_MOD - USE MILLER_MOD - USE OUTUNITS_MOD - USE PARCAL_MOD - USE PARCAL_A_MOD - USE RELADS_MOD - USE RELAX_MOD - USE RESEAU_MOD - USE SPIN_MOD - USE TESTS_MOD - USE TRANS_MOD - USE TL_AED_MOD - USE TYPCAL_MOD - USE TYPCAL_A_MOD - USE TYPEM_MOD - USE TYPEXP_MOD - USE VALIN_MOD - USE XMRHO_MOD -C - DIMENSION VEC(3,3),VB1(3),VB2(3),VB3(3),VBS(3) - DIMENSION ROT(3,3),EMET(3) - DIMENSION VAL2(NATCLU_M) - DIMENSION IRE(NATCLU_M,2) - DIMENSION REL(NATCLU_M),RHOT(NATM) - DIMENSION ATOME(3,NATCLU_M),COORD(3,NATCLU_M) - DIMENSION NTYP(NATCLU_M),NATYP_OLD(NATM) - DIMENSION LMAX_TMP(NATM,NE_M),DIST12(NATCLU_M,NATCLU_M) - DIMENSION IBWD_TMP(NATP_M),RTHFWD_TMP(NATP_M),RTHBWD_TMP(NATP_M) - DIMENSION UJ2_TMP(NATM),RHOT_TMP(NATM),XMT_TMP(NATM) -C - COMPLEX TLSTAR,RHOR(NE_M,NATM,0:18,2,NSPIN2_M) - COMPLEX TLSTAR_A - COMPLEX RHOR_A(0:NT_M,NATM,0:40,2,NSPIN2_M),RAD_D,RAD_E - COMPLEX RHOR1STAR,RHOR2STAR -C -C -C - CHARACTER RIEN - CHARACTER*1 B - CHARACTER*2 R -C -C -C -C -C -C - CHARACTER*30 TUNIT,DUMMY -C - DATA PI,BOHR,SMALL/3.141593,0.529177,0.001/ - DATA INV /1/ -C - READ(*,776) NFICHLEC - READ(*,776) ICOM - DO JF=1,NFICHLEC - READ(*,777) INDATA(JF) - ENDDO -C -C.......... Loop on the data files .......... -C - DO JFICH=1,NFICHLEC - OPEN(UNIT=ICOM, FILE=INDATA(JFICH), STATUS='OLD') - CALL READ_DATA(ICOM,NFICHLEC,JFICH,ITRTL,*2,*1,*55,*74,*99,*504,* - &520,*540,*550,*570,*580,*590,*630) -C -C.......... Atomic case index .......... -C - I_AT=0 - IF((SPECTRO.EQ.'PHD').AND.(I_TEST.EQ.2)) I_AT=1 - IF((SPECTRO.EQ.'AED').AND.(I_TEST_A.EQ.2)) I_AT=1 - IF((SPECTRO.EQ.'XAS').AND.(I_TEST.EQ.2)) I_AT=1 - IF(SPECTRO.EQ.'APC') THEN - IF((I_TEST.EQ.2).AND.(I_TEST_A.EQ.2)) I_AT=1 - ENDIF -C - IF(IBAS.EQ.1) THEN - IF(ITEST.EQ.0) THEN - NEQ=(2*NIV+1)**3 - ELSE - NEQ=(2*NIV+3)**3 - ENDIF - IF(NEQ*NATP_M.GT.NATCLU_M) GOTO 518 - ENDIF -C - IF(SPECTRO.EQ.'APC') THEN - N_EL=2 - ELSE - N_EL=1 - ENDIF - IF((INTERACT.EQ.'COULOMB').OR.(INTERACT.EQ.'DIPCOUL')) THEN - IF(I_MULT.EQ.0) THEN - LE_MIN=ABS(LI_C-ABS(LI_I-LI_A)) - LE_MAX=LI_C+LI_A+LI_I - ELSE - LE_MIN=ABS(LI_C-L_MUL) - LE_MAX=LI_C+L_MUL - ENDIF - ENDIF - IF(SPECTRO.EQ.'EIG') THEN - LE_MIN=1 - LE_MAX=1 - ENDIF -C -C.......... Test of the dimensions against the input values .......... -C - IF(NO.GT.NO_ST_M) GOTO 600 - IF(LE_MAX.GT.LI_M) GOTO 620 -C - OPEN(UNIT=IUI2, FILE=INFILE2, STATUS='OLD') - OPEN(UNIT=IUI3, FILE=INFILE3, STATUS='OLD') - IF(INTERACT.EQ.'DIPCOUL') THEN - OPEN(UNIT=IUI7, FILE=INFILE7, STATUS='OLD') - OPEN(UNIT=IUI8, FILE=INFILE8, STATUS='OLD') - ENDIF -C -C.......... Reading of the TL and radial matrix elements files .......... -C.......... (dipolar excitation case) .......... -C - IF((INTERACT.NE.'COULOMB')) THEN - IF(SPECTRO.EQ.'APC') WRITE(IUO1,418) - READ(IUI2,3) NAT1,NE1,ITL,IPOTC,LMAX_MODE - IF(ISPIN.EQ.0) THEN - IF(NAT1.EQ.1) THEN - WRITE(IUO1,561) - ELSE - WRITE(IUO1,560) NAT1 - ENDIF - ENDIF - IF((ITL.EQ.1).AND.(ISPIN.EQ.1)) THEN - READ(IUI2,530) E_MIN,E_MAX,DE - ENDIF - IF((ISPIN.EQ.0).AND.(ITL.EQ.0)) THEN - NLG=INT(NAT1-0.0001)/4 +1 - DO NN=1,NLG - NRL=4*NN - JD=4*(NN-1)+1 - IF(NN.EQ.NLG) NRL=NAT1 - READ(IUI2,555) (LMAX(JAT,1),JAT=JD,NRL) - WRITE(IUO1,556) (LMAX(JAT,1),JAT=JD,NRL) - ENDDO -C -C Temporary storage of LMAX. Waiting for a version of PHAGEN -C with LMAX dependent on the energy -C - DO JE=1,NE - DO JAT=1,NAT1 - LMAX(JAT,JE)=LMAX(JAT,1) - ENDDO - ENDDO -C - NL1=1 - DO JAT=1,NAT1 - NL1=MAX0(NL1,LMAX(JAT,1)+1) - ENDDO - IF(NL1.GT.NL_M) GOTO 184 - ENDIF - IF(ITL.EQ.0) READ(IUI3,101) NATR,NER - IF(ISPIN.EQ.1) THEN - READ(IUI3,106) L_IN,NATR,NER - IF(LI.NE.L_IN) GOTO 606 - ENDIF - NAT2=NAT+NATA - IF((NAT1.NE.NAT2).OR.(NE1.NE.NE)) GOTO 180 - IF((ITL.EQ.0).AND.((NATR.NE.NAT2).OR.(NER.NE.NE))) GOTO 182 -C -C.......... DL generated by MUFPOT and RHOR given .......... -C.......... by S. M. Goldberg, C. S. Fadley .......... -C.......... and S. Kono, J. Electron Spectr. .......... -C.......... Relat. Phenom. 21, 285 (1981) .......... -C - IF(ITL.EQ.0) THEN - DO JAT=1,NAT2 - IF((INITL.NE.0).AND.(IFTHET.NE.1)) THEN - READ(IUI3,102) RIEN - READ(IUI3,102) RIEN - READ(IUI3,102) RIEN - ENDIF - DO JE=1,NE - IF((IFTHET.EQ.1).OR.(INITL.EQ.0)) GOTO 121 - READ(IUI3,103) ENERGIE - READ(IUI3,102) RIEN - READ(IUI3,102) RIEN - READ(IUI3,102) RIEN - 121 CONTINUE - DO L=0,LMAX(JAT,JE) - READ(IUI2,7) VK(JE),TL(L,1,JAT,JE) - TL(L,1,JAT,JE)=CSIN(TL(L,1,JAT,JE))*CEXP((0., - & 1.)*TL(L,1,JAT,JE)) - ENDDO - IF((IFTHET.EQ.1).OR.(INITL.EQ.0)) GOTO 5 - DO LL=1,18 - READ(IUI3,104) RH1,RH2,DEF1,DEF2 - RHOR(JE,JAT,LL,1,1)=CMPLX(RH1) - RHOR(JE,JAT,LL,2,1)=CMPLX(RH2) - DLT(JE,JAT,LL,1)=CMPLX(DEF1) - DLT(JE,JAT,LL,2)=CMPLX(DEF2) - ENDDO - 5 CONTINUE - ENDDO - ENDDO - ELSE -C -C.......... TL and RHOR calculated by PHAGEN .......... -C - DO JE=1,NE - NLG=INT(NAT2-0.0001)/4 +1 - IF(NE.GT.1) WRITE(IUO1,563) JE - DO NN=1,NLG - NRL=4*NN - JD=4*(NN-1)+1 - IF(NN.EQ.NLG) NRL=NAT2 - READ(IUI2,555) (LMAX(JAT,JE),JAT=JD,NRL) - WRITE(IUO1,556) (LMAX(JAT,JE),JAT=JD,NRL) - ENDDO - NL1=1 - DO JAT=1,NAT2 - NL1=MAX0(NL1,LMAX(JAT,1)+1) - ENDDO - IF(NL1.GT.NL_M) GOTO 184 - DO JAT=1,NAT2 - READ(IUI2,*) DUMMY - DO L=0,LMAX(JAT,JE) - IF(LMAX_MODE.EQ.0) THEN - READ(IUI2,9) VK(JE),TLSTAR - ELSE - READ(IUI2,9) VK(JE),TLSTAR - ENDIF - TL(L,1,JAT,JE)=CONJG(TLSTAR) - VK(JE)=CONJG(VK(JE)) - ENDDO - ENDDO -C - IF((IFTHET.EQ.1).OR.(INITL.EQ.0)) GOTO 333 - IF(JE.EQ.1) THEN - DO JDUM=1,7 - READ(IUI3,102) RIEN - ENDDO - ENDIF - DO JEMET=1,NEMET - JM=IEMET(JEMET) - READ(IUI3,105) RHOR1STAR,RHOR2STAR - RHOR(JE,JM,NNL,1,1)=CONJG(RHOR1STAR) - RHOR(JE,JM,NNL,2,1)=CONJG(RHOR2STAR) - ENDDO - 333 VK(JE)=VK(JE)*A - VK2(JE)=CABS(VK(JE)*VK(JE)) - ENDDO - ENDIF -C - CLOSE(IUI2) - CLOSE(IUI3) -C -C.......... Suppression of possible zeros in the TL array .......... -C.......... (in case of the use of matrix inversion and .......... -C.......... for energy variations) .......... -C - IF((ISPIN.EQ.0).AND.(ITL.EQ.1).AND.(LMAX_MODE.NE.0)) THEN - CALL SUP_ZEROS(TL,LMAX,NE,NAT2,IUO1,ITRTL) - ENDIF - ENDIF -C -C.......... Reading of the TL and radial matrix elements files .......... -C.......... (Coulomb excitation case) .......... -C - IF((INTERACT.EQ.'COULOMB').OR.(INTERACT.EQ.'DIPCOUL')) THEN - IERR=0 - IF(INTERACT.EQ.'COULOMB') THEN - IRD1=IUI2 - IRD2=IUI3 - ELSEIF(INTERACT.EQ.'DIPCOUL') THEN - IRD1=IUI7 - IRD2=IUI8 - ENDIF - IF(SPECTRO.EQ.'APC') WRITE(IUO1,419) - READ(IRD1,3) NAT1_A,NE1_A,ITL_A,IPOTC_A,LMAX_MODE_A - IF(ISPIN.EQ.0) THEN - IF(NAT1_A.EQ.1) THEN - WRITE(IUO1,561) - ELSE - WRITE(IUO1,560) NAT1_A - ENDIF - ENDIF - IF((ITL_A.EQ.1).AND.(ISPIN.EQ.1)) THEN - READ(IRD1,530) E_MIN_A,E_MAX_A,DE_A - ENDIF - IF(ITL_A.EQ.1) THEN - READ(IRD2,107) LI_C2,LI_I2,LI_A2 - READ(IRD2,117) LE_MIN1,N_CHANNEL - LE_MAX1=LE_MIN1+N_CHANNEL-1 - IF(I_TEST_A.NE.1) THEN - IF((LE_MIN.NE.LE_MIN1).OR.(LE_MAX.NE.LE_MAX1)) GOTO - & 610 - ELSE - LI_C2=0 - LI_I2=1 - LI_A2=0 - LE_MIN1=1 - N_CHANNEL=1 - ENDIF - ENDIF - IF((ISPIN.EQ.0).AND.(ITL_A.EQ.0)) THEN - NLG=INT(NAT1_A-0.0001)/4 +1 - DO NN=1,NLG - NRL=4*NN - JD=4*(NN-1)+1 - IF(NN.EQ.NLG) NRL=NAT1_A - READ(IRD1,555) (LMAX_A(JAT,1),JAT=JD,NRL) - WRITE(IUO1,556) (LMAX_A(JAT,1),JAT=JD,NRL) - ENDDO -C -C Temporary storage of LMAX_A. Waiting for a version of PHAGEN -C with LMAX_A dependent on the energy -C - DO JE=1,NE1_A - DO JAT=1,NAT1_A - LMAX_A(JAT,JE)=LMAX_A(JAT,1) - ENDDO - ENDDO -C - NL1_A=1 - DO JAT=1,NAT1_A - NL1_A=MAX0(NL1_A,LMAX_A(JAT,1)+1) - ENDDO - IF(NL1_A.GT.NL_M) GOTO 184 - ENDIF - IF(ITL_A.EQ.0) READ(IRD2,101) NATR_A,NER_A - IF(ISPIN.EQ.1) THEN - READ(IRD2,106) L_IN_A,NATR_A,NER_A - IF(LI_C.NE.L_IN_A) GOTO 606 - ENDIF - NAT2_A=NAT+NATA - NAT2=NAT2_A - IF((NAT1_A.NE.NAT2_A).OR.(NE1_A.NE.NE_A)) GOTO 180 - IF((ITL_A.EQ.0).AND.((NATR_A.NE.NAT2_A).OR.(NER_A.NE.NE))) - & GOTO 182 -C -C.......... DL generated by MUFPOT and RHOR given .......... -C.......... by S. M. Goldberg, C. S. Fadley .......... -C.......... and S. Kono, J. Electron Spectr. .......... -C.......... Relat. Phenom. 21, 285 (1981) .......... -C - IF(ITL_A.EQ.0) THEN - CONTINUE - ELSE -C -C.......... TL_A and RHOR_A calculated by PHAGEN .......... -C - DO JE=1,NE_A - NLG=INT(NAT2_A-0.0001)/4 +1 - IF(NE_A.GT.1) WRITE(IUO1,563) JE - DO NN=1,NLG - NRL=4*NN - JD=4*(NN-1)+1 - IF(NN.EQ.NLG) NRL=NAT2_A - READ(IRD1,555) (LMAX_A(JAT,JE),JAT=JD,NRL) - WRITE(IUO1,556) (LMAX_A(JAT,JE),JAT=JD,NRL) - ENDDO - DO JAT=1,NAT2_A - READ(IRD1,*) DUMMY - DO L=0,LMAX_A(JAT,JE) - IF(LMAX_MODE_A.EQ.0) THEN - READ(IRD1,9) VK_A(JE),TLSTAR - ELSE - READ(IRD1,7) VK_A(JE),TLSTAR - ENDIF - TL_A(L,1,JAT,JE)=CONJG(TLSTAR) - VK_A(JE)=CONJG(VK_A(JE)) - ENDDO - ENDDO -C - IF(IFTHET_A.EQ.1) GOTO 331 - DO LE=LE_MIN,LE_MAX - DO JEMET=1,NEMET - JM=IEMET(JEMET) - READ(IRD2,109) L_E,LB_MIN,LB_MAX - IF(I_TEST_A.EQ.1) THEN - L_E=1 - LB_MIN=0 - LB_MAX=1 - ENDIF - IF(LE.NE.L_E) IERR=1 - L_BOUNDS(L_E,1)=LB_MIN - L_BOUNDS(L_E,2)=LB_MAX - DO LB=LB_MIN,LB_MAX - READ(IRD2,108) L_A,RAD_D,RAD_E - RHOR_A(LE,JM,L_A,1,1)=RAD_D - RHOR_A(LE,JM,L_A,2,1)=RAD_E - IF(I_TEST_A.EQ.1) THEN - IF(LB.EQ.LB_MIN) THEN - RHOR_A(LE,JM,L_A,1,1)=(0.0,0.0) - RHOR_A(LE,JM,L_A,2,1)=(1.0,0.0) - ELSEIF(LB.EQ.LB_MAX) THEN - RHOR_A(LE,JM,L_A,1,1)=(1.0,0.0) - RHOR_A(LE,JM,L_A,2,1)=(0.0,0.0) - ENDIF - ENDIF - ENDDO - ENDDO - ENDDO - 331 VK_A(JE)=VK_A(JE)*A - VK2_A(JE)=CABS(VK_A(JE)*VK_A(JE)) - ENDDO - ENDIF -C - CLOSE(IRD1) - CLOSE(IRD2) -C -C.......... Suppression of possible zeros in the TL array .......... -C.......... (in case of the use of matrix inversion and .......... -C.......... for energy variations) .......... -C - IF((ISPIN.EQ.0).AND.(ITL_A.EQ.1).AND.(LMAX_MODE_A.NE.0)) THEN - CALL SUP_ZEROS(TL_A,LMAX_A,NE_A,NAT2_A,IUO1,ITRTL) - ENDIF - IF(SPECTRO.EQ.'APC') WRITE(IUO1,420) -C - ENDIF -C -C.......... Check of the consistency of the two TL and radial .......... -C.......... matrix elements for APECS .......... -C - IF(SPECTRO.EQ.'APC') THEN -C - I_TL_FILE=0 - I_RD_FILE=0 -C - IF(NAT1.NE.NAT1_A) I_TL_FILE=1 - IF(NE1.NE.NE1_A) I_TL_FILE=1 - IF(ITL.NE.ITL_A) I_TL_FILE=1 - IF(IPOTC.NE.IPOTC_A) I_TL_FILE=1 -C - IF(LI_C.NE.LI_C2) I_RD_FILE=1 - IF(LI_I.NE.LI_I2) I_RD_FILE=1 - IF(LI_A.NE.LI_A2) I_RD_FILE=1 -C - IF(I_TL_FILE.EQ.1) GOTO 608 - IF(I_RD_FILE.EQ.1) GOTO 610 - IF(IERR.EQ.1) GOTO 610 -C - ENDIF -C -C.......... Calculation of the scattering factor (only) .......... -C - IF((IFTHET.EQ.0).AND.(IFTHET_A.EQ.0)) GO TO 8 - IF(IFTHET.EQ.1) THEN - CALL PLOTFD(A,LMAX,ITL,NL1,NAT2,NE) - ELSEIF(IFTHET_A.EQ.1) THEN -c CALL PLOTFD_A(A,LMAX_A,ITL_A,NL1_A,NAT2_A,NE_A) - ENDIF - WRITE(IUO1,57) - STOP -C - 8 IF(IBAS.EQ.0) THEN -C -C............... Reading of an external cluster ............... -C -C -C Cluster originating from CLUSTER_NEW.F : IPHA=0 -C Cluster originating from PHAGEN_NEW.F : IPHA=1 (atomic units), IPHA=2 (angstroems) -C Other cluster : the first line must be text; then -C free format : Atomic number,X,Y,Z,number -C of the corresponding prototypical atom ; -C All atoms corresponding to the same -C prototypical atom must follow each other. -C Moreover, the blocks of equivalent atoms -C must be ordered by increasing number of -C prototypical atom. -C - VALZ_MIN=1000.0 - VALZ_MAX=-1000.0 -C - OPEN(UNIT=IUI4, FILE=INFILE4, STATUS='OLD') - READ(IUI4,778,ERR=892) IPHA - GOTO 893 - 892 IPHA=3 - IF(UNIT.EQ.'ANG') THEN - CUNIT=1./A - TUNIT='ANGSTROEMS' - ELSEIF(UNIT.EQ.'LPU') THEN - CUNIT=1. - TUNIT='UNITS OF THE LATTICE PARAMETER' - ELSEIF(UNIT.EQ.'ATU') THEN - CUNIT=BOHR/A - TUNIT='ATOMIC UNITS' - ELSE - GOTO 890 - ENDIF - 893 NATCLU=0 - DO JAT=1,NAT2 - NATYP(JAT)=0 - ENDDO - IF(IPHA.EQ.0) THEN - CUNIT=1. - TUNIT='UNITS OF THE LATTICE PARAMETER' - ELSEIF(IPHA.EQ.1) THEN - CUNIT=BOHR/A - TUNIT='ATOMIC UNITS' - IEMET(1)=1 - ELSEIF(IPHA.EQ.2) THEN - CUNIT=1./A - TUNIT='ANGSTROEMS' - IEMET(1)=1 - ENDIF - IF(IPRINT.EQ.2) THEN - IF(I_AT.NE.1) THEN - WRITE(IUO1,558) IUI4,TUNIT - IF(IPHA.EQ.3) WRITE(IUO1,549) - ENDIF - ENDIF - JATM=0 - DO JLINE=1,10000 - IF(IPHA.EQ.0) THEN - READ(IUI4,125,END=780) R,NN,X,Y,Z,JAT - ELSEIF(IPHA.EQ.1) THEN - READ(IUI4,779,END=780) R,NN,X,Y,Z,JAT - ELSEIF(IPHA.EQ.2) THEN - READ(IUI4,779,END=780) R,NN,X,Y,Z,JAT - ELSEIF(IPHA.EQ.3) THEN - READ(IUI4,*,END=780) NN,X,Y,Z,JAT - ENDIF - JATM=MAX0(JAT,JATM) - NATCLU=NATCLU+1 - IF(IPHA.NE.3) THEN - CHEM(JAT)=R - ELSE - CHEM(JAT)='XX' - ENDIF - NZAT(JAT)=NN - NATYP(JAT)=NATYP(JAT)+1 - COORD(1,NATCLU)=X*CUNIT - COORD(2,NATCLU)=Y*CUNIT - COORD(3,NATCLU)=Z*CUNIT - VALZ(NATCLU)=Z*CUNIT -c IF((IPRINT.GE.2).AND.(I_AT.EQ.0)) THEN - IF(IPRINT.GE.2) THEN - WRITE(IUO1,557) NATCLU,COORD(1,NATCLU),COORD(2, - & NATCLU),COORD(3,NATCLU),JAT,NATYP(JAT),CHEM(JAT) - ENDIF - ENDDO - 780 NBZ=NATCLU - IF(JATM.NE.NAT) GOTO 514 - CLOSE(IUI4) -C - IF(NATCLU.GT.NATCLU_M) GOTO 510 - DO JA1=1,NATCLU - DO JA2=1,NATCLU - DIST12(JA1,JA2)=SQRT((COORD(1,JA1)-COORD(1,JA2))**2+( - & COORD(2,JA1)-COORD(2,JA2))**2+(COORD(3,JA1)-COORD(3,JA2))** - & 2) - IF((JA2.GT.JA1).AND.(DIST12(JA1,JA2).LT.0.001)) GOTO - & 895 - ENDDO - ENDDO -C - D_UP=VALZ_MAX-VALZ(1) - D_DO=VALZ(1)-VALZ_MIN - IF((D_DO.LE.D_UP).AND.(I_GR.EQ.2)) THEN - I_INV=1 - ELSE - I_INV=0 - ENDIF - ELSE -C -C............... Construction of an internal cluster ............... -C - CALL BASE - CALL ROTBAS(ROT) - IF(IVG0.EQ.2) THEN - NMAX=NIV+1 - ELSE - NMAX=(2*NIV+1)**3 - ENDIF - IF((IPRINT.EQ.2).AND.(IVG0.LE.1)) THEN - WRITE(IUO1,37) - WRITE(IUO1,38) NIV - DO NUM=1,NMAX - CALL NUMAT(NUM,NIV,IA,IB,IC) - WRITE(IUO1,17) NUM,IA,IB,IC - ENDDO - WRITE(IUO1,39) - ENDIF - CALL AMAS(NIV,ATOME,COORD,VALZ,IESURF,COUPUR,ROT,IRE,NATYP, - & NBZ,NAT2,NCOUCH,NMAX) - IF((IREL.GE.1).OR.(NRELA.GT.0)) THEN - CALL RELA(NBZ,NPLAN,NAT2,VALZ,VAL2,VAL,COORD,NATYP,REL, - & NCOUCH) - IF(IREL.EQ.1) THEN - DO JP=1,NPLAN - VAL(JP)=VAL2(JP) - ENDDO - ENDIF - ENDIF - ENDIF -C -C Storage of the extremal values of x and y for each plane. They define -C the exterior of the cluster when a new cluster has to be build to -C support a point-group -C - IF(I_GR.GE.1) THEN - IF((IREL.EQ.0).OR.(IBAS.EQ.0)) THEN - CALL ORDRE(NBZ,VALZ,NPLAN,VAL) - WRITE(IUO1,50) NPLAN - DO K=1,NPLAN - WRITE(IUO1,29) K,VAL(K) - X_MAX(K)=0. - X_MIN(K)=0. - Y_MAX(K)=0. - Y_MIN(K)=0. - ENDDO - ENDIF - DO JAT=1,NATCLU - X=COORD(1,JAT) - Y=COORD(2,JAT) - Z=COORD(3,JAT) - DO JPLAN=1,NPLAN - IF(ABS(Z-VAL(JPLAN)).LT.SMALL) THEN - X_MAX(JPLAN)=MAX(X,X_MAX(JPLAN)) - X_MIN(JPLAN)=MIN(X,X_MIN(JPLAN)) - Y_MAX(JPLAN)=MAX(Y,Y_MAX(JPLAN)) - Y_MIN(JPLAN)=MIN(Y,Y_MIN(JPLAN)) - ENDIF - ENDDO - ENDDO - ENDIF -C -C Instead of the symmetrization of the cluster (this version only) -C - N_PROT=NAT - NAT_ST=0 - DO JTYP=1,JATM - NB_AT=NATYP(JTYP) - IF(NB_AT.GT.NAT_EQ_M) GOTO 614 - DO JA=1,NB_AT - NAT_ST=NAT_ST+1 - NCORR(JA,JTYP)=NAT_ST - ENDDO - ENDDO - DO JC=1,3 - DO JA=1,NATCLU - SYM_AT(JC,JA)=COORD(JC,JA) - ENDDO - ENDDO -C -C Checking surface-like atoms for mean square displacements -C calculations -C - CALL CHECK_VIB(NAT2) -C -C.......... Set up of the variables used for an internal .......... -C.......... calculation of the mean free path and/or of .......... -C.......... the mean square displacements .......... -C - IF((IDCM.EQ.1).OR.(ILPM.EQ.1)) THEN - DO JTYP=1,NAT2 - XMT(JTYP)=XMAT(NZAT(JTYP)) - RHOT(JTYP)=RHOAT(NZAT(JTYP)) - ENDDO - XMTA=XMT(1) - RHOTA=RHOT(1) - NZA=NZAT(1) - ENDIF - IF(IDCM.GT.0) THEN - CALL CHNOT(3,VECBAS,VEC) - DO J=1,3 - VB1(J)=VEC(J,1) - VB2(J)=VEC(J,2) - VB3(J)=VEC(J,3) - ENDDO - CPR=1. - CALL PRVECT(VB2,VB3,VBS,CPR) - VM=PRSCAL(VB1,VBS) - QD=(6.*PI*PI*NAT/VM)**(1./3.) - ENDIF -C -C.......... Writing of the contents of the cluster, .......... -C.......... of the position of the different planes .......... -C.......... and of their respective absorbers in .......... -C.......... the control file IUO1 .......... -C - IF(I_AT.EQ.1) GOTO 153 - IF((IPRINT.EQ.2).AND.(IBAS.GT.0)) THEN - WRITE(IUO1,40) - NCA=0 - DO J=1,NAT - DO I=1,NMAX - NCA=NCA+1 - WRITE(IUO1,20) J,I - WRITE(IUO1,21) (ATOME(L,NCA),L=1,3) - K=IRE(NCA,1) - IF(K.EQ.0) THEN - WRITE(IUO1,22) - ELSE - WRITE(IUO1,23) (COORD(L,K),L=1,3),IRE(NCA,2) - ENDIF - ENDDO - ENDDO - WRITE(IUO1,41) - ENDIF - IF(IBAS.EQ.1) THEN - WRITE(IUO1,24) - NATCLU=0 - DO I=1,NAT - NN=NATYP(I) - NATCLU=NATCLU+NATYP(I) - WRITE(IUO1,26) NN,I - ENDDO - IF(IADS.EQ.1) NATCLU=NATCLU+NADS1+NADS2+NADS3 - WRITE(IUO1,782) NATCLU - IF(NATCLU.GT.NATCLU_M) GOTO 516 - IF(IPRINT.EQ.3) WRITE(IUO1,559) - IF(IPRINT.EQ.3) THEN - NBTA=0 - DO JT=1,NAT2 - NBJT=NATYP(JT) - DO JN=1,NBJT - NBTA=NBTA+1 - WRITE(IUO1,557) NBTA,COORD(1,NBTA),COORD(2,NBTA), - & COORD(3,NBTA),JT,JN,CHEM(JT) - ENDDO - ENDDO - ENDIF - ENDIF - 153 IF((ITEST.EQ.1).AND.(IBAS.GT.0)) THEN - CALL TEST(NIV,ROT,NATYP,NBZ,NAT2,IESURF,COUPUR,*56) - ENDIF - IF((IREL.EQ.0).OR.(IBAS.EQ.0)) THEN - CALL ORDRE(NBZ,VALZ,NPLAN,VAL) - IF(I_AT.EQ.0) WRITE(IUO1,50) NPLAN - DO K=1,NPLAN - IF(I_AT.EQ.0) WRITE(IUO1,29) K,VAL(K) - ENDDO - ENDIF -C - IF(SPECTRO.NE.'EIG') THEN - IF(I_AT.EQ.0) WRITE(IUO1,30) - IF((IPRINT.GT.0).AND.(I_AT.EQ.0)) THEN - WRITE(IUO1,31) (IEMET(J),J=1,NEMET) - ENDIF - ZEM=1.E+20 - DO L=1,NPLAN - Z=VAL(L) - DO JEMED=1,NEMET - CALL EMETT(JEMED,IEMET,Z,COORD,NATYP,EMET,NTEM,JNEM,* - & 93) - IF(I_AT.EQ.0) WRITE(IUO1,34) L,NTEM,EMET(1),EMET(2), - & EMET(3) - IF((IPHA.EQ.1).OR.(IPHA.EQ.2)) ZEM=EMET(3) - GO TO 33 - 93 IF(I_AT.EQ.0) WRITE(IUO1,94) L,NTEM - 33 CONTINUE - ENDDO - ENDDO - ENDIF -C -C.......... Loop on the electrons involved in the .......... -C.......... spectroscopy : N_EL = 1 for PHD, XAS .......... -C.......... or AED and N_EL = 2 for APC .......... -C - DO J_EL=1,N_EL -C -C.......... Writing the information on the spectroscopies .......... -C.......... in the control file IUO1 .......... -C - IF(SPECTRO.EQ.'EIG') WRITE(IUO1,252) - IF(SPECTRO.EQ.'XAS') GOTO 566 - IF((SPECTRO.EQ.'APC').AND.(J_EL.EQ.1)) THEN - IF(IPHI.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,236) - ELSE - WRITE(IUO1,248) - ENDIF - ENDIF - IF(ITHETA.EQ.1) WRITE(IUO1,245) - IF(I_TEST.EQ.1) WRITE(IUO1,234) - ENDIF -C -C---------- Photoelectron diffraction case (PHD) ---------- -C - IF((SPECTRO.EQ.'PHD').OR.(SPECTRO.EQ.'APC')) THEN - IF(SPECTRO.EQ.'PHD') THEN - IF(IPHI.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,35) - ELSE - WRITE(IUO1,246) - ENDIF - ENDIF - IF(ITHETA.EQ.1) WRITE(IUO1,44) - IF(IE.EQ.1) WRITE(IUO1,58) - IF(INITL.EQ.0) WRITE(IUO1,118) - IF(I_TEST.EQ.1) WRITE(IUO1,234) - ENDIF - IF((SPECTRO.EQ.'APC').AND.(J_EL.EQ.1)) THEN - WRITE(IUO1,418) - WRITE(IUO1,18) - ENDIF - IF(J_EL.EQ.2) GOTO 222 - IF(IPRINT.GT.0) THEN - WRITE(IUO1,92) - WRITE(IUO1,91) - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,335) - ELSE - WRITE(IUO1,336) - ENDIF - WRITE(IUO1,91) - IF(IPOTC.EQ.0) THEN - WRITE(IUO1,339) - ELSE - WRITE(IUO1,334) - ENDIF - WRITE(IUO1,91) - IF(INITL.NE.0) THEN - WRITE(IUO1,337) - WRITE(IUO1,91) - IF(IPOL.EQ.0) THEN - WRITE(IUO1,88) - ELSEIF(ABS(IPOL).EQ.1) THEN - WRITE(IUO1,87) - ELSEIF(IPOL.EQ.2) THEN - WRITE(IUO1,89) - ENDIF - WRITE(IUO1,91) - IF(IDICHR.GT.0) THEN - WRITE(IUO1,338) - ENDIF - WRITE(IUO1,91) - WRITE(IUO1,92) - WRITE(IUO1,90) - WRITE(IUO1,43) THLUM,PHILUM - IF((SPECTRO.EQ.'PHD').AND.(IMOD.EQ.1)) THEN - WRITE(IUO1,45) - ENDIF - ENDIF -C - IF(INITL.EQ.2) THEN - WRITE(IUO1,79) LI,LI-1,LI+1 - IF(I_SO.EQ.1) THEN - WRITE(IUO1,80) S_O - ENDIF - DO JE=1,NE - DO JEM=1,NEMET - JTE=IEMET(JEM) - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,111) JTE,RHOR(JE,JTE,NNL, - & 1,1),RHOR(JE,JTE,NNL,2,1) - IF(ITL.EQ.0) THEN - WRITE(IUO1,444) JTE,DLT(JE,JTE, - & NNL,1),DLT(JE,JTE,NNL,2) - ENDIF - ENDIF - ENDDO - ENDDO - ELSEIF(INITL.EQ.-1) THEN - WRITE(IUO1,82) LI,LI-1 - IF(I_SO.EQ.1) THEN - WRITE(IUO1,80) S_O - ENDIF - DO JE=1,NE - DO JEM=1,NEMET - JTE=IEMET(JEM) - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,113) JTE,RHOR(JE,JTE,NNL, - & 1,1) - IF(ITL.EQ.0) THEN - WRITE(IUO1,445) JTE,DLT(JE,JTE, - & NNL,1) - ENDIF - ENDIF - ENDDO - ENDDO - ELSEIF(INITL.EQ.1) THEN - WRITE(IUO1,82) LI,LI+1 - IF(I_SO.EQ.1) THEN - WRITE(IUO1,80) S_O - ENDIF - DO JE=1,NE - DO JEM=1,NEMET - JTE=IEMET(JEM) - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,113) JTE,RHOR(JE,JTE,NNL, - & 2,1) - IF(ITL.EQ.0) THEN - WRITE(IUO1,445) JTE,DLT(JE,JTE, - & NNL,2) - ENDIF - ENDIF - ENDDO - ENDDO - ENDIF -C - IF(I_AT.EQ.0) THEN - IF(INV.EQ.0) THEN - IF(NDIF.EQ.1) THEN - IF(ISPHER.EQ.1) THEN - WRITE(IUO1,83) - ELSEIF(ISPHER.EQ.0) THEN - WRITE(IUO1,84) - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,97) NDIF - ELSE - WRITE(IUO1,98) NDIF - ENDIF - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,122) - ELSE - WRITE(IUO1,120) - ENDIF - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,85) - ELSE - WRITE(IUO1,86) - ENDIF - ENDIF -C - ENDIF - 222 CONTINUE - ENDIF -C -C---------- Auger diffraction case (AED) ---------- -C - IF((SPECTRO.EQ.'AED').OR.(SPECTRO.EQ.'APC')) THEN - IF(SPECTRO.EQ.'AED') THEN - IF(IPHI_A.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,235) - ELSE - WRITE(IUO1,247) - ENDIF - ENDIF - IF(ITHETA_A.EQ.1) WRITE(IUO1,244) - IF(I_TEST_A.EQ.1) WRITE(IUO1,234) - ENDIF - IF((SPECTRO.EQ.'APC').AND.(J_EL.EQ.2)) THEN - WRITE(IUO1,419) - WRITE(IUO1,18) - ENDIF - IF((SPECTRO.EQ.'AED').OR.(J_EL.EQ.2)) THEN - IF(IPRINT.GT.0) THEN - WRITE(IUO1,92) - WRITE(IUO1,91) - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,335) - ELSE - WRITE(IUO1,336) - ENDIF - WRITE(IUO1,91) - IF(IPOTC_A.EQ.0) THEN - WRITE(IUO1,339) - ELSE - WRITE(IUO1,334) - ENDIF - WRITE(IUO1,91) - WRITE(IUO1,92) - WRITE(IUO1,95) AUGER - CALL AUGER_MULT - IF(I_MULT.EQ.0) THEN - WRITE(IUO1,154) - ELSE - WRITE(IUO1,155) MULTIPLET - ENDIF -C - DO JEM=1,NEMET - JTE=IEMET(JEM) - WRITE(IUO1,112) JTE - DO LE=LE_MIN,LE_MAX - WRITE(IUO1,119) LE - LA_MIN=L_BOUNDS(LE,1) - LA_MAX=L_BOUNDS(LE,2) - DO LA=LA_MIN,LA_MAX - IF(ISPIN.EQ.0) THEN - WRITE(IUO1,115) LA,RHOR_A(LE,JTE, - & LA,1,1),RHOR_A(LE,JTE,LA,2,1) - ENDIF - ENDDO - ENDDO - ENDDO -C - IF(I_AT.EQ.0) THEN - IF(INV.EQ.0) THEN - IF(NDIF.EQ.1) THEN - IF(ISPHER.EQ.1) THEN - WRITE(IUO1,83) - ELSEIF(ISPHER.EQ.0) THEN - WRITE(IUO1,84) - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,97) NDIF - ELSE - WRITE(IUO1,98) NDIF - ENDIF - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,122) - ELSE - WRITE(IUO1,120) - ENDIF - ENDIF - ELSE - IF(ISPHER.EQ.0) THEN - WRITE(IUO1,85) - ELSE - WRITE(IUO1,86) - ENDIF - ENDIF -C - ENDIF - ENDIF - ENDIF -C -C.......... Check of the dimensioning of the treatment routine .......... -C - CALL STOP_TREAT(NFICHLEC,NPLAN,NEMET,NE,NTHETA,NTHETA_A,NPHI, - & NPHI_A,ISOM,I_EXT,I_EXT_A,SPECTRO) -C -C.......... Call of the subroutine performing either .......... -C.......... the PhD, AED, EXAFS or APECS calculation .......... -C - 566 IF(ISPIN.EQ.0) THEN - IF(SPECTRO.EQ.'EIG') THEN - CALL EIGDIF_MI - ELSEIF(SPECTRO.EQ.'AED') THEN -c CALL AEDDIF_SE(NPLAN,VAL,ZEM,IPHA,NAT2,COORD,NATYP,RHOR_A, -c 1 NATCLU,NFICHLEC,JFICH,NP,LE_MIN,LE_MAX) - ELSEIF(SPECTRO.EQ.'XAS') THEN -c CALL XASDIF_SE(NPLAN,VAL,ZEM,IPHA,RHOR,NFICHLEC,JFICH,NP) - ELSEIF(SPECTRO.EQ.'APC') THEN -c IF(J_EL.EQ.1) THEN -c CALL PHDDIF_SE(NPLAN,VAL,ZEM,IPHA,NAT2,COORD,NATYP,RHOR, -c 1 NATCLU,NFICHLEC,JFICH,NP) -c ELSEIF(J_EL.EQ.2) THEN -c CALL AEDDIF_SE(NPLAN,VAL,ZEM,IPHA,NAT2,COORD,NATYP,RHOR_A, -c 1 NATCLU,NFICHLEC,JFICH,NP,LE_MIN,LE_MAX) -c ENDIF - ENDIF - ELSEIF(ISPIN.EQ.1) THEN -c IF(SPECTRO.EQ.'PHD') THEN -c CALL PHDDIF_SP(NPLAN,VAL,ZEM,IPHA,NAT2,COORD,NATYP,RHOR, -c 1 NATCLU,NFICHLEC,JFICH,NP) -c ELSEIF(SPECTRO.EQ.'AED') THEN -c CALL AEDDIF_SP -c ELSEIF(SPECTRO.EQ.'XAS') THEN -c CALL XASDIF_SP -c ENDIF - continue - ENDIF -C -C.......... End of the MS calculation : .......... -C.......... direct exit or treatment of the results .......... -C -C -C.......... End of the loop on the electrons .......... -C - ENDDO -C - IF(SPECTRO.EQ.'PHD') THEN - IF(IPHI.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,52) - ELSE - WRITE(IUO1,249) - ENDIF - ENDIF - IF(ITHETA.EQ.1) WRITE(IUO1,49) - IF(IE.EQ.1) WRITE(IUO1,59) - ELSEIF(SPECTRO.EQ.'XAS') THEN - WRITE(IUO1,51) - ELSEIF(SPECTRO.EQ.'AED') THEN - IF(IPHI_A.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,237) - ELSE - WRITE(IUO1,250) - ENDIF - ENDIF - IF(ITHETA_A.EQ.1) WRITE(IUO1,238) - ELSEIF(SPECTRO.EQ.'APC') THEN - IF(IPHI.EQ.1) THEN - IF(STEREO.EQ.' NO') THEN - WRITE(IUO1,239) - ELSE - WRITE(IUO1,251) - ENDIF - ENDIF - IF(ITHETA.EQ.1) WRITE(IUO1,240) - ELSEIF(SPECTRO.EQ.'EIG') THEN - WRITE(IUO1,253) - ENDIF -C - CLOSE(ICOM) - IF((NFICHLEC.GT.1).AND.(ISOM.NE.0)) THEN - WRITE(IUO1,562) - ENDIF - IF(ISOM.EQ.0) CLOSE(IUO2) - IF((ISOM.EQ.0).AND.(NFICHLEC.NE.1)) CLOSE(IUO1) -C -C.......... End of the loop on the data files .......... -C - ENDDO -C - IF(ISOM.NE.0) THEN - JFF=1 - IF(ISPIN.EQ.0) THEN - IF((SPECTRO.EQ.'PHD').OR.(SPECTRO.EQ.'AED')) THEN -c CALL TREAT_PHD(ISOM,NFICHLEC,JFF,NP) - ELSEIF(SPECTRO.EQ.'XAS') THEN -c CALL TREAT_XAS(ISOM,NFICHLEC,NP) - ENDIF - ELSEIF(ISPIN.EQ.1) THEN -c IF((SPECTRO.EQ.'PHD').OR.(SPECTRO.EQ.'AED')) THEN -c CALL TREAT_PHD_SP(ISOM,NFICHLEC,JFF,NP) -c ELSEIF(SPECTRO.EQ.'XAS') THEN -c CALL TREAT_XAS_SP(ISOM,NFICHLEC,NP) -c ENDIF - continue - ENDIF - ENDIF -C - IF((ISOM.NE.0).OR.(NFICHLEC.EQ.1)) CLOSE(IUO1) - IF(ISOM.NE.0) CLOSE(IUO2) - STOP -C - 1 WRITE(IUO1,60) - STOP - 2 WRITE(IUO1,61) - STOP - 55 WRITE(IUO1,65) - STOP - 56 WRITE(IUO1,64) - STOP - 74 WRITE(IUO1,75) - STOP - 99 WRITE(IUO1,100) - STOP - 180 WRITE(IUO1,181) - STOP - 182 WRITE(IUO1,183) - STOP - 184 WRITE(IUO1,185) - STOP - 504 WRITE(IUO1,505) - STOP - 510 WRITE(IUO1,511) IUI4 - STOP - 514 WRITE(IUO1,515) - STOP - 516 WRITE(IUO1,517) - STOP - 518 WRITE(IUO1,519) - WRITE(IUO1,889) - STOP - 520 WRITE(IUO1,521) - STOP - 540 WRITE(IUO1,541) - STOP - 550 WRITE(IUO1,551) - STOP - 570 WRITE(IUO1,571) - STOP - 580 WRITE(IUO1,581) - STOP - 590 WRITE(IUO1,591) - STOP - 600 WRITE(IUO1,601) - STOP - 602 WRITE(IUO1,603) - STOP - 604 WRITE(IUO1,605) - STOP - 606 WRITE(IUO1,607) - STOP - 608 WRITE(IUO1,609) - STOP - 610 WRITE(IUO1,611) - STOP - 614 WRITE(IUO1,615) NB_AT - STOP - 620 WRITE(IUO1,621) LE_MAX - STOP - 630 WRITE(IUO1,631) - STOP - 890 WRITE(IUO1,891) - STOP - 895 WRITE(IUO1,896) JA1,JA2 -C - 3 FORMAT(5(5X,I4)) - 7 FORMAT(3X,F9.4,1X,F9.4,5X,F12.9,5X,F12.9) - 9 FORMAT(3X,F9.4,1X,F9.4,5X,E12.6,5X,E12.6) - 17 FORMAT(12X,'ATOM NUMBER ',I4,10X,'CORRESPONDING TRANSLATIONS ',': - & (',I3,',',I3,',',I3,')') - 18 FORMAT(' ',/) - 20 FORMAT(/,7X,'ATOM OF TYPE ',I2,' AND OF NUMBER ',I5) - 21 FORMAT(17X,'COORDINATES IN THE TOTAL CLUSTER : (',F7.3,',',F7.3, - &',',F7.3,')') - 22 FORMAT(22X,'THIS ATOM HAS BEEN SUPRESSED IN THE REDUCED CLUSTER') - 23 FORMAT(17X,'COORDINATES IN THE REDUCED CLUSTER :(',F7.3,',',F7.3, - &',',F7.3,')',5X,'NEW NUMBER : ',I4) - 24 FORMAT(///,29X,'CONTENTS OF THE REDUCED CLUSTER :',/) - 26 FORMAT(28X,I4,' ATOMS OF TYPE ',I2) - 29 FORMAT(/,20X,'THE Z POSITION OF PLANE ',I3,' IS : ',F6.3) - 30 FORMAT(///,23X,'THE ABSORBING ATOMS ARE OF TYPE :',/) - 31 FORMAT(38X,10(I2,3X),//) - 34 FORMAT(//,2X,'PLANE No ',I3,3X,'THE ABSORBER OF TYPE ', I2,' IS - &POSITIONED AT (',F7.3,',',F7.3,',',F7.3,')') - 35 FORMAT(/////,'########## BEGINNING ', 'OF THE AZIMUTHAL - &PHOTOELECTRON DIFFRACTION CALCULATION #####', '#####',/////) - 36 FORMAT(/////,'########## BEGINNING ', 'OF THE - &EXAFS CALCULATION ##########',/////) - 37 FORMAT(/////,'++++++++++++++++++++', ' NUMBERING OF THE - &ATOMS GENERATED +++++++++++++++++++') - 38 FORMAT(///,30X,'TRANSLATION LEVEL : ',I2,///) - 39 FORMAT(///,'++++++++++++++++++++++++++++++++++++++++++++++++', - & '++++++++++++++++++++++++++++++++',/////) - 40 FORMAT(/////,'======================', ' CONTENTS OF THE - &REDUCED CLUSTER ======================',///) - 41 FORMAT(///,'==================================================== - &','============================',/////) - 43 FORMAT(14X,'TH_LIGHT = ',F6.2,' DEGREES',5X,'PHI_LIGHT = ',F6.2, - &' DEGREES') - 44 FORMAT(/////,'########## BEGINNING ', 'OF THE POLAR - &PHOTOELECTRON DIFFRACTION CALCULATION #####', '#####',/////) - 45 FORMAT(14X,' (WHEN THE DETECTOR IS ALONG ','THE NORMAL TO THE - &SURFACE)') - 49 FORMAT(/////,'########## END OF THE ', 'POLAR PHOTOELECTRON - &DIFFRACTION CALCULATION ##########') - 50 FORMAT(///,22X,'THE CLUSTER IS COMPOSED OF ',I2,' PLANES :') - 51 FORMAT(/////,'########## END OF THE ', 'EXAFS - &CALCULATION ##########') - 52 FORMAT(/////,'########## END OF THE ', 'AZIMUTHAL PHOTOELECTRON - &DIFFRACTION CALCULATION #####','#####') - 57 FORMAT(///,27X,'CALCULATION OF THE SCATTERING FACTOR DONE') - 58 FORMAT(/////,'########## BEGINNING ', 'OF THE FINE - &STRUCTURE OSCILLATIONS CALCULATION #####', '#####',/////) - 59 FORMAT(/////,'########## END OF THE ', 'FINE STRUCTURE - &OSCILLATIONS CALCULATION #####','#####') - 60 FORMAT(///,'<<<<<<<<<< (NAT,NE,NEMET) > (NATP_M,NE_M,','NEMET_M) - & - CHECK THE DIMENSIONING >>>>>>>>>>') - 61 FORMAT(///,22X,' <<<<<<<<<< THIS STRUCTURE DOES NOT EXIST ', - &' >>>>>>>>>>') - 64 FORMAT(///,4X,' <<<<<<<<<< NIV IS TOO SMALL, THE REDUCED ', - &'CLUSTER HAS NOT CONVERGED YET >>>>>>>>>>') - 65 FORMAT(///,4X,' <<<<<<<<<< ONLY ONE OF THE VALUES IPHI,ITHETA ', - & 'ET IE CAN BE EQUAL TO 1 >>>>>>>>>>') - 75 FORMAT(///,8X,' <<<<<<<<<< CHANGE THE DIMENSIONING OF PCREL ', - & 'IN MAIN ET READ_DATA >>>>>>>>>>') - 79 FORMAT(//,18X,'INITIAL STATE L = ',I1,5X,'FINAL STATES L = ', - & I1,',',I1,/) - 80 FORMAT(15X,'(SPIN-ORBIT COMPONENT OF THE INITIAL CORE STATE : ', - &A3,')',//) - 81 FORMAT(18X,'(BOTH SPIN-ORBIT COMPONENTS TAKEN INTO ACCOUNT)') - 82 FORMAT(//,21X,'INITIAL STATE L = ',I1,5X,'FINAL STATE L = ',I1) - 83 FORMAT(//,32X,'(SPHERICAL WAVES)') - 84 FORMAT(//,34X,'(PLANE WAVES)') - 85 FORMAT(//,26X,'(PLANE WAVES - ATOMIC CASE)') - 86 FORMAT(//,24X,'(SPHERICAL WAVES - ATOMIC CASE)') - 87 FORMAT(24X,'+ LINEARLY POLARIZED LIGHT +') - 88 FORMAT(24X,'+ NON POLARIZED LIGHT +') - 89 FORMAT(24X,'+ CIRCULARLY POLARIZED LIGHT +') - 90 FORMAT(////,31X,'POSITION OF THE LIGHT :',/) - 91 FORMAT(24X,'+',35X,'+') - 92 FORMAT(24X,'+++++++++++++++++++++++++++++++++++++') - 94 FORMAT(//,2X,'PLANE No ',I3,3X,'NO ABSORBER OF TYPE ',I2, ' IS - &PRESENT IN THIS PLANE') - 95 FORMAT(////,31X,'AUGER LINE :',A6,//) - 97 FORMAT(///,19X,'(PLANE WAVES MULTIPLE SCATTERING - ORDER ',I1,') - &') - 98 FORMAT(///,17X,'(SPHERICAL WAVES MULTIPLE SCATTERING - ORDER ', - &I1,')') - 100 FORMAT(///,8X,'<<<<<<<<<< WRONG NAME FOR THE INITIAL STATE',' - &>>>>>>>>>>') - 101 FORMAT(24X,I3,24X,I3) - 102 FORMAT(A1) - 103 FORMAT(31X,F7.2) - 104 FORMAT(29X,F8.5,4X,F8.5,7X,F8.5,4X,F8.5) - 105 FORMAT(1X,E12.5,1X,E12.5,2X,E12.5,1X,E12.5,4X,E12.5,1X,E12.5,2X, - &E12.5,1X,E12.5,2X,E12.5,1X,E12.5,4X,A9) - 106 FORMAT(12X,I3,12X,I3,12X,I3) - 107 FORMAT(5X,I2,5X,I2,5X,I2) - 108 FORMAT(19X,I2,8X,F8.5,1X,F8.5,4X,F8.5,1X,F8.5) - 109 FORMAT(5X,I2,12X,I2,11X,I2) - 110 FORMAT(16X,'RADIAL MATRIX ELEMENTS FOR THE ABSORBER OF TYPE ',I2, - &' :',/,22X,'(THE SPIN DOUBLET IS GIVEN AS : OUT/IN)',//) - 111 FORMAT(6X,'RADIAL MATRIX ELEMENTS FOR THE ABSORBER OF TYPE ',I2, - &' : (',F8.5,',',F8.5,')',/,59X,'(',F8.5,',',F8.5,')') - 112 FORMAT(6X,'RADIAL MATRIX ELEMENTS FOR THE ABSORBER OF TYPE ',I2, - &' : ',/,8X,'(LE : ALLOWED VALUES FOR ESCAPING AUGER',' ELECTRON) - &',/,8X,'(L : INTERNAL VALUE THAT WILL BE SUMMED ON)',//) - 113 FORMAT(6X,'RADIAL MATRIX ELEMENT FOR THE ABSORBER OF ', - * 'TYPE ',I2,' : (',F8.5,',',F8.5,')') - 114 FORMAT(/) - 115 FORMAT(15X,'L = ',I2,5X,'(',F8.5,',',F8.5,')',5X,'(',F8.5,',',F8. - &5,')') - 117 FORMAT(12X,I2,5X,I2) - 118 FORMAT(/,37X,'AUGER ELECTRON DIFFRACTION',/) - 119 FORMAT(10X,'LE = ',I2,11X,'DIRECT INTEGRAL',8X,'EXCHANGE - &INTEGRAL') - 120 FORMAT(///,15X,'(SPHERICAL WAVES MULTIPLE SCATTERING - MATRIX ', - &'INVERSION)') - 122 FORMAT(///,17X,'(PLANE WAVES MULTIPLE SCATTERING - MATRIX ', - &'INVERSION)') - 125 FORMAT(11X,A2,5X,I2,3F10.4,12X,I4) - 154 FORMAT(///,20X,'CALCULATION MADE FOR THE FULL AUGER LINE',' ',/, - &' ',/,' ') - 155 FORMAT(///,20X,'CALCULATION MADE FOR THE ',A3,' MULTIPLET ', - &'LINE',' ',/,' ',/,' ') - 181 FORMAT(///,'<<<<<<<<<< NAT OR NE DIFFERENT BETWEEN THE INPUT ', - &'AND PHASE SHIFTS FILES >>>>>>>>>>') - 183 FORMAT(///,'<<<<<<<<<< NAT OR NE DIFFERENT BETWEEN THE INPUT ', - &'AND RADIAL MATRIX ELEMENTS FILES >>>>>>>>>>') - 185 FORMAT(///,'<<<<<<<<<< LMAX > NL_M-1 IN THE PHASE SHIFTS ', - &'FILE >>>>>>>>>>') - 234 FORMAT(' -----> TEST CALCULATION : NO EXCITATION ','MATRIX - &ELEMENTS TAKEN INTO ACCOUNT <-----',///) - 235 FORMAT(/////,'########## BEGINNING ', 'OF THE AZIMUTHAL - &AUGER DIFFRACTION CALCULATION #####', '#####',/////) - 236 FORMAT(/////,'########## BEGINNING ', 'OF THE AZIMUTHAL - &APECS DIFFRACTION CALCULATION #####', '#####',/////) - 237 FORMAT(/////,'########## END ', 'OF THE AZIMUTHAL AUGER - &DIFFRACTION CALCULATION #####', '#####',/////) - 238 FORMAT(/////,6X,'########## END ', 'OF THE POLAR AUGER - &DIFFRACTION CALCULATION #####', '#####',/////) - 239 FORMAT(/////,'########## END ', 'OF THE AZIMUTHAL APECS - &DIFFRACTION CALCULATION #####', '#####',/////) - 240 FORMAT(/////,6X,'########## END ', 'OF THE POLAR APECS - &DIFFRACTION CALCULATION #####', '#####',/////) - 244 FORMAT(/////,6X,'########## BEGINNING ', 'OF THE POLAR AUGER - &DIFFRACTION CALCULATION #####', '#####',/////) - 245 FORMAT(/////,6X,'########## BEGINNING ', 'OF THE POLAR APECS - &DIFFRACTION CALCULATION #####', '#####',/////) - 246 FORMAT(/////,'########## BEGINNING ', 'OF THE FULL ANGLE - &PHOTOELECTRON DIFFRACTION CALCULATION ','##########',/////) - 247 FORMAT(/////,'########## BEGINNING ', 'OF THE FULL ANGLE - &AUGER DIFFRACTION CALCULATION ', '##########',/////) - 248 FORMAT(/////,'########## BEGINNING ', 'OF THE FULL ANGLE - &APECS DIFFRACTION CALCULATION ', '##########',/////) - 249 FORMAT(/////,'########## END OF THE ', 'FULL ANGLE PHOTOELECTRON - &DIFFRACTION CALCULATION #####','#####') - 250 FORMAT(/////,'########## END ', 'OF THE FULL ANGLE AUGER - &DIFFRACTION CALCULATION #####', '#####',/////) - 251 FORMAT(/////,'########## END ', 'OF THE FULL ANGLE APECS - &DIFFRACTION CALCULATION #####', '#####',/////) - 252 FORMAT(/////,'########## BEGINNING ', 'OF THE MULTIPLE - &SCATTERING EIGENVALUE CALCULATION #####', '#####',/////) - 253 FORMAT(/////,'########## END ', 'OF THE MULTIPLE SCATTERING - &EIGENVALUE CALCULATION #####', '#####',/////) - 334 FORMAT(24X,'+ COMPLEX POTENTIAL CALCULATION +') - 335 FORMAT(24X,'+ STANDARD +') - 336 FORMAT(24X,'+ SPIN-POLARIZED +') - 337 FORMAT(24X,'+ WITH +') - 338 FORMAT(24X,'+ IN DICHROIC MODE +') - 339 FORMAT(24X,'+ REAL POTENTIAL CALCULATION +') - 418 FORMAT(///,9X,'------------------------ FIRST ELECTRON : ','---- - &--------------------') - 419 FORMAT(///,9X,'------------------------ SECOND ELECTRON : ','---- - &--------------------') - 420 FORMAT(///,9X,'----------------------------------------------','- - &---------------------') - 444 FORMAT(12X,'PHASE SHIFTS FOR THE ABSORBER OF TYPE ',I2,' : ','( - &',F8.5,',',F8.5,')',/,56X,'(',F8.5,',',F8.5,')') - 445 FORMAT(12X,'PHASE SHIFT FOR THE ABSORBER OF TYPE ',I2,' : (',F8. - &5,',',F8.5,')') - 505 FORMAT(///,'<<<<<<<<<< LI IS LARGER THAN LI_M - ','CHECK THE - &DIMENSIONING >>>>>>>>>>') - 511 FORMAT(///,'<<<<<<<<<< NATCLU_M IN THE .inc FILE IS NOT ', - &'CONSISTENT WITH THE NUMBER OF ATOMS READ FROM UNIT ',I2,' - &>>>>>>>>>>') - 515 FORMAT(///,'<<<<<<<<<< INCOMPATIBILITY BETWEEN THE VALUES OF ', - &'NAT IN THE DATA AND CLUSTER FILES >>>>>>>>>>') - 517 FORMAT(///,'<<<<<<<<<< THERE ARE MISSING VALUES FOR THFWD AND ', - &'IBWD >>>>>>>>>>') - 519 FORMAT(///,'<<<<<<<<<< NATCLU_M IN THE .inc FILE IS NOT',' - &CONSISTENT WITH THE NUMBER OF ATOMS GENERATED BY THE ','CODE - &>>>>>>>>>>') - 521 FORMAT(///,'<<<<<<<<<< SPIN-ORBIT COMPONENT NOT CONSISTENT - &WITH',' THE VALUE OF LI >>>>>>>>>>') - 530 FORMAT(3X,F9.4,3X,F9.4,3X,F9.4) - 535 FORMAT(29X,F8.5,1X,F8.5) - 541 FORMAT(///,'<<<<<<<<<< THE NUMBER OF LINES THFWD DOES NOT ', - &'CORRESPOND TO NAT >>>>>>>>>>') - 543 FORMAT(5X,F12.9,5X,F12.9) - 549 FORMAT(//,14X,' No ',10X,'COORDINATES',9X,'TYPE',2X,'SNo',2X, - &'SYM',/) - 551 FORMAT(///,'<<<<<<<<<< THE NUMBER OF LINES UJ2 DOES NOT ', - &'CORRESPOND TO NAT >>>>>>>>>>') - 555 FORMAT(4(7X,I2)) - 556 FORMAT(28X,4(I2,5X)) - 557 FORMAT(13X,I4,3X,'(',F7.3,',',F7.3,',',F7.3,')',2X,I4,2X,I4,3X, - &A2) - 558 FORMAT(/////,18X,'CONTENTS OF THE CLUSTER READ FROM UNIT ',I2,' : - & ',/,20X,'READ IN ',A30,//,15X,'No',13X,'(X,Y,Z)',10X,'CLASS',1X, - &'ATOM',/) - 559 FORMAT(/////,25X,'CONTENTS OF THE CLUSTER GENERATED : ',//,14X,' - &No ',10X,'COORDINATES',9X,'TYPE',2X,'SNo',2X,'SYM',/) - 560 FORMAT(////,12X,'MAXIMAL VALUES OF L FOR THE ',I3,' PROTOTYPICAL - &ATOMS : ',//) - 561 FORMAT(////,18X,'MAXIMAL VALUE OF L FOR THE ','PROTOTYPICAL ATOM - &: ',//) - 562 FORMAT(///,'oooooooooooooooo',12X,'END OF THE INPUT DATA FILE', - &13X,'oooooooooooooooo',///) - 563 FORMAT(//,20X,'ENERGY POINT No ',I3,' :',/) - 571 FORMAT(///,'<<<<<<<<<< THE NUMBER OF LINES ATBAS DOES NOT ', - &'CORRESPOND TO NAT >>>>>>>>>>') - 581 FORMAT(///,'<<<<<<<<<< LI OR IMOD NOT CONSISTENT BETWEEN ','PHD - &AND AED FOR COINCIDENCE CALCULATION >>>>>>>>>>') - 591 FORMAT(///,'<<<<<<<<<< THE EXTERNAL DIRECTIONS FILE IS ','NOT - &CONSISTENT WITH THE INPUT DATA FILE >>>>>>>>>>') - 601 FORMAT(///,'<<<<<<<<<< NO_ST_M IS TOO SMALL IN THE .inc FILE ', - &'>>>>>>>>>>',//) - 603 FORMAT(///,'<<<<<<<<<< NSPIN_M OR NSPIN2_M IS TOO SMALL IN THE - &','.inc FILE >>>>>>>>>>',//) - 605 FORMAT(///,'<<<<<<<<<< NT_M IS TOO SMALL IN THE .inc FILE ', - &'>>>>>>>>>>',//) - 607 FORMAT(///,'<<<<<<<<<< THE INITIAL STATE LI IN THE INPUT DATA - &','FILE IS DIFFERENT FROM THAT IN THE RADIAL MATRIX ','ELEMENTS - &FILE >>>>>>>>>>',//) - 609 FORMAT(///,'<<<<<<<<<< THE TWO TL FILE ARE NOT COMPATIBLE ', - &'>>>>>>>>>>',//) - 611 FORMAT(///,3X,'<<<<<<<<<< THE RADIAL FILE FOR THE AUGER ', - &'ELECTRON IS NOT COMPATIBLE >>>>>>>>>>',/,3X,'<<<<<<<<<< ', - &17X,'WITH THE INPUT DATA FILE ',16X,'>>>>>>>>>>',//) - 613 FORMAT(///,'<<<<<<<<<< NATP_M SHOULD BE AT LEAST ',I3,' IN ', - &'THE DIMENSIONNING FILE >>>>>>>>>>',//) - 615 FORMAT(///,'<<<<<<<<<< NAT_EQ_M SHOULD BE AT LEAST ',I3,' IN ', - &'THE DIMENSIONNING FILE >>>>>>>>>>',//) - 621 FORMAT(///,'<<<<<<<<<< LI_M SHOULD BE AT LEAST ',I3,' IN ', - &'THE DIMENSIONNING FILE >>>>>>>>>>',//) - 631 FORMAT(///,'<<<<<<<<<< EXCURSIONS OF ANGLES SHOULD ',' BE - &IDENTICAL >>>>>>>>>>',/,'<<<<<<<<<< ','FOR BOTH - &ELECTRONS IN CLUSTER ROTATION MODE',' >>>>>>>>>>',//) - 776 FORMAT(I2) - 777 FORMAT(A24) - 778 FORMAT(30X,I1) - 779 FORMAT(11X,A2,5X,I2,3F10.4,I5) - 782 FORMAT(/////,22X,'THE CLUSTER GENERATED CONSISTS OF : ',I4,' - &ATOMS') - 889 FORMAT(/////,'<<<<<<<<<< DECREASE NIV OR INCREASE',' NATCLU_M - &>>>>>>>>>>') - 891 FORMAT(/////,'<<<<<<<<<< WRONG NAME FOR THE COORDINATES ''', - &'UNITS >>>>>>>>>>') - 896 FORMAT(///,10X,'<<<<<<<<<< ERROR IN THE COORDINATES OF THE',' - &ATOMS >>>>>>>>>>',/,10X,'<<<<<<<<<< ATOMS ',I4,' AND ',I4,' - &ARE IDENTICAL >>>>>>>>>>') -C - END -C -C======================================================================= -C - SUBROUTINE DWSPH(JTYP,JE,X,TLT,ISPEED) -C -C This routine recomputes the T-matrix elements taking into account the -C mean square displacements. -C -C When the argument X is tiny, no vibrations are taken into account -C -C Last modified : 25 Apr 2013 -C -C INCLUDE 'spec.inc' - USE DIM_MOD - USE TRANS_MOD -C - DIMENSION GNT(0:N_GAUNT) -C - COMPLEX TLT(0:NT_M,4,NATM,NE_M),SL1,ZEROC -C - COMPLEX*16 FFL(0:2*NL_M) -C - DATA PI4,EPS /12.566371,1.0E-10/ -C - ZEROC=(0.,0.) -C - IF(X.GT.EPS) THEN -C -C Standard case: vibrations -C - IF(ISPEED.LT.0) THEN - NSUM_LB=ABS(ISPEED) - ENDIF -C - COEF=PI4*EXP(-X) - NL2=2*LMAX(JTYP,JE)+2 - IBESP=5 - MG1=0 - MG2=0 -C - CALL BESPHE(NL2,IBESP,X,FFL) -C - DO L=0,LMAX(JTYP,JE) - XL=FLOAT(L+L+1) - SL1=ZEROC -C - DO L1=0,LMAX(JTYP,JE) - XL1=FLOAT(L1+L1+1) - CALL GAUNT(L,MG1,L1,MG2,GNT) - L2MIN=ABS(L1-L) - IF(ISPEED.GE.0) THEN - L2MAX=L1+L - ELSEIF(ISPEED.LT.0) THEN - L2MAX=L2MIN+2*(NSUM_LB-1) - ENDIF - SL2=0. -C - DO L2=L2MIN,L2MAX,2 - XL2=FLOAT(L2+L2+1) - C=SQRT(XL1*XL2/(PI4*XL)) - SL2=SL2+C*GNT(L2)*REAL(DREAL(FFL(L2))) - ENDDO -C - SL1=SL1+SL2*TL(L1,1,JTYP,JE) - ENDDO -C - TLT(L,1,JTYP,JE)=COEF*SL1 -C - ENDDO -C - ELSE -C -C Argument X tiny: no vibrations -C - DO L=0,LMAX(JTYP,JE) -C - TLT(L,1,JTYP,JE)=TL(L,1,JTYP,JE) -C - ENDDO -C - ENDIF -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE EIG_MAT_MS(JE,E_KIN) -C -C This subroutine stores the G_o T kernel matrix and computes -C its eigenvalues to check the larger ones. If one or more -C of these eigenvalues are larger than or equal to 1, no -C series expansion can converge -C -C Last modified : 20 Jul 2009 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - USE COOR_MOD - USE OUTFILES_MOD - USE OUTUNITS_MOD - USE TRANS_MOD -C - PARAMETER(NLTWO=2*NL_M) -C - CHARACTER*24 OUTFILE,PATH -C - COMPLEX*16 HL1(0:NLTWO),SM(LINMAX*NATCLU_M,LINMAX*NATCLU_M) - COMPLEX*16 SUM_L,IC,ZEROC,WORK(32*LINMAX*NATCLU_M) - COMPLEX*16 YLM(0:NLTWO,-NLTWO:NLTWO),TLK,EXPKJ - COMPLEX*16 W(LINMAX*NATCLU_M) - COMPLEX*16 VL(1,1),VR(1,1) -C - DOUBLE PRECISION RWORK(2*LINMAX*NATCLU_M) -C - REAL*8 PI,ATTKJ,GNT(0:N_GAUNT),XKJ,YKJ,ZKJ,RKJ,ZDKJ,KRKJ -C - REAL W1(LINMAX*NATCLU_M),W2(LINMAX*NATCLU_M) -C - DATA PI,SMALL /3.1415926535898D0,0.0001/ -C - IC=(0.D0,1.D0) - ZEROC=(0.D0,0.D0) - IBESS=3 - NPRINT=10 -C - WRITE(IUO1,5) - IF(JE.EQ.1) THEN -C -C Name of the second output file where all the eigenvalues -C will be written -C - IOUT2=55 - IOUT3=60 - N_DOT=1 - DO J_CHAR=1,24 - IF(OUTFILE2(J_CHAR:J_CHAR).EQ.'.') GOTO 888 - N_DOT=N_DOT+1 - ENDDO - 888 CONTINUE - OUTFILE=OUTFILE2(1:N_DOT)//'egv' - PATH=OUTFILE2(1:N_DOT)//'pth' - OPEN(UNIT=IOUT2, FILE=OUTFILE, STATUS='UNKNOWN') - OPEN(UNIT=IOUT3, FILE=PATH, STATUS='UNKNOWN') - ENDIF -C -C Construction of the multiple scattering kernel matrix G_o T. -C Elements are stored using a linear index LINJ -C representing (J,LJ) -C - JLIN=0 - DO JTYP=1,N_PROT - NBTYPJ=NATYP(JTYP) - LMJ=LMAX(JTYP,JE) - DO JNUM=1,NBTYPJ - JATL=NCORR(JNUM,JTYP) - XJ=SYM_AT(1,JATL) - YJ=SYM_AT(2,JATL) - ZJ=SYM_AT(3,JATL) -C - DO LJ=0,LMJ - DO MJ=-LJ,LJ - JLIN=JLIN+1 -C - KLIN=0 - DO KTYP=1,N_PROT - NBTYPK=NATYP(KTYP) - LMK=LMAX(KTYP,JE) - DO KNUM=1,NBTYPK - KATL=NCORR(KNUM,KTYP) - IF(KATL.NE.JATL) THEN - XKJ=DBLE(SYM_AT(1,KATL)-XJ) - YKJ=DBLE(SYM_AT(2,KATL)-YJ) - ZKJ=DBLE(SYM_AT(3,KATL)-ZJ) - RKJ=DSQRT(XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ) - KRKJ=DBLE(VK(JE))*RKJ - ATTKJ=DEXP(-DIMAG(DCMPLX(VK(JE)))* - & RKJ) - EXPKJ=(XKJ+IC*YKJ)/RKJ - ZDKJ=ZKJ/RKJ - CALL SPH_HAR2(2*NL_M,ZDKJ,EXPKJ,YLM, - & LMJ+LMK) - CALL BESPHE2(LMJ+LMK+1,IBESS,KRKJ, - & HL1) - ENDIF -C - DO LK=0,LMK - L_MIN=ABS(LK-LJ) - L_MAX=LK+LJ - TLK=DCMPLX(TL(LK,1,KTYP,JE)) - DO MK=-LK,LK - KLIN=KLIN+1 - SM(KLIN,JLIN)=ZEROC - SUM_L=ZEROC - IF(KATL.NE.JATL) THEN - CALL GAUNT2(LK,MK,LJ,MJ,GNT) -C - DO L=L_MIN,L_MAX,2 - M=MJ-MK - IF(ABS(M).LE.L) THEN - SUM_L=SUM_L+(IC**L)* - & HL1(L)*YLM(L,M)*GNT(L) - ENDIF - ENDDO - SUM_L=SUM_L*ATTKJ*4.D0*PI*IC - ELSE - SUM_L=ZEROC - ENDIF -C - SM(KLIN,JLIN)=TLK*SUM_L -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - N_DIM=LINMAX*NATCLU_M -C -C Eigenvalues of the kernel multiple scattering matrix SM -C - CALL ZGEEV('N','N',JLIN,SM,N_DIM,W,VL,1,VR,1,WORK,34*N_DIM,RWORK, - &INFO) - IF(INFO.NE.0) THEN - WRITE(IUO1,*) ' ' - WRITE(IUO1,*) ' ---> WORK(1),INFO =',WORK(1),INFO - WRITE(IUO1,*) ' ' - ENDIF -C - N_EIG=0 -C - WRITE(IOUT2,75) - WRITE(IOUT2,110) - WRITE(IOUT2,80) E_KIN - WRITE(IOUT2,110) - WRITE(IOUT2,75) - WRITE(IOUT2,105) - WRITE(IOUT2,75) - XMAX_L=0. - N_XMAX=0 -C - DO LIN=1,JLIN - EIG=REAL(CDABS(W(LIN))) - WRITE(IOUT2,100) REAL(DBLE(W(LIN))),REAL(DIMAG(W(LIN))),EIG - IF((EIG-XMAX_L).GT.0.0001) N_XMAX=LIN - XMAX_L=MAX(XMAX_L,EIG) - W1(LIN)=EIG - IF(EIG.GT.1.000) THEN - N_EIG=N_EIG+1 - ENDIF - ENDDO -C - WRITE(IOUT2,75) - WRITE(IOUT2,85) XMAX_L - WRITE(IOUT2,90) N_XMAX - WRITE(IOUT2,95) W(N_XMAX) - WRITE(IOUT2,75) -C - CALL ORDRE(JLIN,W1,NFIN,W2) -C -C - WRITE(IUO1,10) - WRITE(IUO1,10) - WRITE(IUO1,15) W2(1) - WRITE(IUO1,20) W2(NFIN) - WRITE(IUO1,10) - WRITE(IUO1,10) - IF(N_EIG.GE.1) THEN - IF(N_EIG.EQ.1) THEN - WRITE(IUO1,25) JLIN - ELSE - WRITE(IUO1,30) N_EIG,JLIN - ENDIF - ENDIF -C - WRITE(IUO1,65) N_XMAX - WRITE(IUO1,70) W(N_XMAX) - WRITE(IUO1,10) - WRITE(IOUT3,100) REAL(DBLE(W(N_XMAX))),REAL(DIMAG(W(N_XMAX))) -C - NPR=NPRINT/5 - WRITE(IUO1,10) - WRITE(IUO1,10) - WRITE(IUO1,35) 5*NPR - WRITE(IUO1,10) -C - DO JP=0,NPR-1 - J=5*JP - WRITE(IUO1,40) W2(J+1),W2(J+2),W2(J+3),W2(J+4),W2(J+5) - ENDDO - WRITE(IUO1,10) - WRITE(IUO1,10) - WRITE(IUO1,45) W2(1) - WRITE(IUO2,*) E_KIN,W2(1) - IF(N_EIG.EQ.0) THEN - WRITE(IUO1,50) - ELSE - WRITE(IUO1,55) - ENDIF - WRITE(IUO1,10) - WRITE(IUO1,10) - WRITE(IUO1,60) -C - RETURN -C - 5 FORMAT(/,11X,'----------------- EIGENVALUE ANALYSIS ','--------- - &--------') - 10 FORMAT(11X,'-',54X,'-') - 15 FORMAT(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-') - 20 FORMAT(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-') - 25 FORMAT(11X,'-',6X,'1 EIGENVALUE IS > 1 ON A TOTAL OF ',I8,6X,'-') - 30 FORMAT(11X,'-',4X,I5,' EIGENVALUES ARE > 1 ON A TOTAL OF ',I8,2X, - &'-') - 35 FORMAT(11X,'-',11X,'THE ',I3,' LARGER EIGENVALUES ARE :',11X,'-') - 40 FORMAT(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-') - 45 FORMAT(11X,'-',5X,'SPECTRAL RADIUS OF THE KERNEL MATRIX :',F6.3, - &5X,'-') - 50 FORMAT(11X,'-',5X,'---> THE MULTIPLE SCATTERING SERIES ', - &'CONVERGES',4X,'-') - 55 FORMAT(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/ - &,11X,'-',18X,'SCATTERING SERIES',19X,'-') - 60 FORMAT(11X,'----------------------------------------','---------- - &------',/) - 65 FORMAT(11X,'-',5X,' LABEL OF LARGEST EIGENVALUE : ',I5,8X,'- - &') - 70 FORMAT(11X,'-',5X,' LARGEST EIGENVALUE : ','(',F6.3,',',F6.3, - &')',8X,'-') - 75 FORMAT(' ') - 80 FORMAT(' KINETIC ENERGY : ',F7.2,' eV') - 85 FORMAT(' LARGEST MODULUS OF EIGENVALUE : ',F6.3) - 90 FORMAT(' LABEL OF LARGEST EIGENVALUE : ',I5) - 95 FORMAT(' LARGEST EIGENVALUE : (',F6.3,',',F6.3,') - &') - 100 FORMAT(5X,F7.3,2X,F7.3,2X,F6.3) - 105 FORMAT(7X,'EIGENVALUES :',3X,'MODULUS :') - 110 FORMAT(2X,'-------------------------------') -C - END -C -C======================================================================= -C - SUBROUTINE FACDIF1(VKE,RJ,RJK,THRJ,PHIRJ,BETA,GAMMA,L,M,FSPH,JAT, - &JE,*) -C -C This routine computes a spherical wave scattering factor -C -C Last modified : 03/04/2006 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - USE APPROX_MOD - USE EXPFAC_MOD - USE TRANS_MOD - USE TYPCAL_MOD I2 => IPHI, I3 => IE, I4 => ITHETA, I5 => IMOD, - &I6 => IPOL, I7 => I_CP, I8 => I_EXT, I9 => I_TEST - - DIMENSION PLMM(0:100,0:100) - DIMENSION D(1-NL_M:NL_M-1,1-NL_M:NL_M-1,0:NL_M-1) -C - COMPLEX HLM(0:NO_ST_M,0:NL_M-1),HLN(0:NO_ST_M,0:NL_M-1),FSPH,RHOJ - COMPLEX HLM1,HLM2,HLM3,HLM4,ALMU,BLMU,SLP,SNU,SMU,VKE - COMPLEX RHOJK -C -C - DATA PI/3.141593/ -C - A=1. - INTER=0 - IF(ITL.EQ.1) VKE=VK(JE) - RHOJ=VKE*RJ - RHOJK=VKE*RJK - HLM1=(1.,0.) - HLM2=(1.,0.) - HLM3=(1.,0.) - HLM4=(1.,0.) - IEM=1 - CSTH=COS(BETA) - IF((IFTHET.EQ.0).OR.(THRJ.LT.0.0001)) THEN - INTER=1 - BLMU=SQRT(4.*PI/FLOAT(2*L+1))*CEXP((0.,-1.)*M*(PHIRJ-PI)) - ENDIF - CALL PLM(CSTH,PLMM,LMAX(JAT,JE)) - IF(ISPHER.EQ.0) NO1=0 - IF(ISPHER.EQ.1) THEN - IF(NO.EQ.8) THEN - NO1=LMAX(JAT,JE)+1 - ELSE - NO1=NO - ENDIF - CALL POLHAN(ISPHER,NO1,LMAX(JAT,JE),RHOJ,HLM) - IF(IEM.EQ.0) THEN - HLM4=HLM(0,L) - ENDIF - IF(RJK.GT.0.0001) THEN - NDUM=0 - CALL POLHAN(ISPHER,NDUM,LMAX(JAT,JE),RHOJK,HLN) - ENDIF - CALL DJMN(THRJ,D,L) - A1=ABS(D(0,M,L)) - IF(((A1.LT.0.0001).AND.(IFTHET.EQ.1)).AND.(INTER.EQ.0)) - & RETURN 1 - ENDIF - MUMAX=MIN0(L,NO1) - SMU=(0.,0.) - DO 10 MU=0,MUMAX - IF(MOD(MU,2).EQ.0) THEN - B=1. - ELSE - B=-1. - IF(SIN(BETA).LT.0.) THEN - A=-1. - ENDIF - ENDIF - IF(ISPHER.LE.1) THEN - ALMU=(1.,0.) - C=1. - ENDIF - IF(ISPHER.EQ.0) GOTO 40 - IF(INTER.EQ.0) BLMU=CMPLX(D(M,0,L)) - IF(MU.GT.0) THEN - C=B*FLOAT(L+L+1)/EXPF(MU,L) - ALMU=(D(M,MU,L)*CEXP((0.,-1.)*MU*GAMMA)+B* - * CEXP((0.,1.)*MU*GAMMA)*D(M,-MU,L))/BLMU - ELSE - C=1. - ALMU=CMPLX(D(M,0,L))/BLMU - ENDIF - 40 SNU=(0.,0.) - NU1=INT(0.5*(NO1-MU)+0.0001) - NUMAX=MIN0(NU1,L-MU) - DO 20 NU=0,NUMAX - SLP=(0.,0.) - LPMIN=MAX0(MU,NU) - DO 30 LP=LPMIN,LMAX(JAT,JE) - IF(ISPHER.EQ.1) THEN - HLM1=HLM(NU,LP) - IF(RJK.GT.0.0001) HLM3=HLN(0,LP) - ENDIF - SLP=SLP+FLOAT(2*LP+1)*TL(LP,1,JAT,JE)*HLM1*PLMM(LP, - & MU)*HLM3 - 30 CONTINUE - IF(ISPHER.EQ.1) THEN - HLM2=HLM(MU+NU,L) - ENDIF - SNU=SNU+SLP*HLM2 - 20 CONTINUE - SMU=SMU+SNU*C*ALMU*A*B - 10 CONTINUE - FSPH=SMU/(VKE*HLM4) -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE PLOTFD(A,LMX,ITL,NL,NAT,NE) -C -C This routine prepares the output for a plot of the scattering factor -C -C INCLUDE 'spec.inc' -C -C - USE APPROX_MOD - USE FDIF_MOD - USE INIT_L_MOD L => LI, I2 => INITL, I3 => NNL, I4 => LF1, I5 => - &LF2, I10 => ISTEP_LF - USE INIT_J_MOD - USE OUTFILES_MOD - USE OUTUNITS_MOD - USE PARCAL_MOD N3 => NPHI, N4 => NE, N5 => NTHETA, N6 => NEPS - USE TYPCAL_MOD I7 => IFTHET, I8 => IMOD, I9 => IPOL, I12 => I_CP, - & I13 => I_EXT, I14 => I_TEST - USE VALIN_MOD U1 => THLUM, U2 => PHILUM, U3 => ELUM, N7 => NONVOL - USE VALFIN_MOD -C - DIMENSION LMX(NATM,NE_M) -C - COMPLEX FSPH,VKE -C - DATA PI,CONV/3.141593,0.512314/ -C - OPEN(UNIT=IUO3, FILE=OUTFILE3, STATUS='UNKNOWN') - IF(ISPHER.EQ.0) THEN - L=0 - LMAX=0 - ELSE - LMAX=L - ENDIF - PHITOT=360. - THTOT=360.*ITHETA*(1-IPHI)+180.*ITHETA*IPHI - NPHI=(NFTHET+1)*IPHI+(1-IPHI) - NTHT=(NFTHET+1)*ITHETA*(1-IPHI)+(NFTHET/2+1)*ITHETA*IPHI+ - * (1-ITHETA) - NE=NFTHET*IE + (1-IE) - WRITE(IUO3,1) ISPHER,NL,NAT,L,NTHT,NPHI,NE,E0,EFIN - DO 10 JT=1,NTHT - DTHETA=THETA1+FLOAT(JT-1)*THTOT/FLOAT(MAX0(NTHT-1,1)) - RTHETA=DTHETA*PI/180. - TEST=SIN(RTHETA) - IF(TEST.GE.0.) THEN - POZ=PI - EPS=1. - ELSE - POZ=0. - EPS=-1. - ENDIF - BETA=RTHETA*EPS - IF(ABS(TEST).LT.0.0001) THEN - NPHIM=1 - ELSE - NPHIM=NPHI - ENDIF - DO 20 JP=1,NPHIM - DPHI=PHI1+FLOAT(JP-1)*PHITOT/FLOAT(MAX0(NPHI-1,1)) - RPHI=DPHI*PI/180. - GAMMA=POZ-RPHI - DO 30 JE=1,NE - IF(NE.EQ.1) THEN - ECIN=E0 - ELSE - ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1) - ENDIF - IF(ITL.EQ.0) VKE=SQRT(ECIN-ABS(VINT))*CONV*A*(1.,0.) - DO 40 JAT=1,NAT - IF(L.GT.LMX(JAT,JE)) GOTO 90 - DO 50 M=-LMAX,LMAX - CALL FACDIF1(VKE,R1,R2,THETA0,PHI0,BETA, - & GAMMA,L,M,FSPH,JAT,JE,*60) - GOTO 70 - 60 WRITE(IUO1,80) - STOP - 70 REFTH=REAL(FSPH) - XIMFTH=AIMAG(FSPH) - WRITE(IUO3,5) JE,JAT,L,M,REFTH,XIMFTH,DTHETA, - & DPHI,ECIN - 50 CONTINUE - GOTO 40 - 90 WRITE(IUO1,100) JAT - STOP - 40 CONTINUE - 30 CONTINUE - 20 CONTINUE - 10 CONTINUE - CLOSE(IUO3) - 1 FORMAT(5X,I1,2X,I2,2X,I4,2X,I2,2X,I3,2X,I3,2X,I3,2X,F8.2,2X,F8.2) - 5 FORMAT(1X,I3,1X,I4,1X,I2,1X,I3,1X,F6.3,1X,F6.3,1X,F6.2,1X,F6.2, - &1X,F8.2) - 80 FORMAT(15X,'<<<<< WRONG VALUE OF THETA0 : THE DENOMINATOR ','IS - &ZERO >>>>>') - 100 FORMAT(15X,'<<<<< THE VALUE OF L EST IS TOO LARGE FOR ATOM',' : - &',I2,' >>>>>') -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE EIGDIF_MI -C -C This subroutine computes the eigenvalues of the -C multiple scattering matrix to obtain its spectral radius -C -C The eigenvalue calculation is performed using the LAPACK -C eigenvalue routines for a general complex matrix (I_PWM = 0) -C or using the power method (I_PWM > 0) -C -C Last modified : 26 Apr 2013 -C -C INCLUDE 'spec.inc' - USE DIM_MOD - USE CONVTYP_MOD - USE COOR_MOD NTCLU => NATCLU, NTP => NATYP - USE DEBWAL_MOD - USE EIGEN_MOD NE => NE_EIG, E0 => E0_EIG, EFIN => EFIN_EIG - USE OUTFILES_MOD - USE OUTUNITS_MOD - USE RESEAU_MOD - USE TESTS_MOD - USE TRANS_MOD - USE VALIN_MOD E1 => E0, PHLUM => PHILUM -C - COMPLEX IC,ONEC - COMPLEX TLT(0:NT_M,4,NATM,NE_M) -C -C - DATA CONV /0.512314/ -C - IC=(0.,1.) - ONEC=(1.,0.) -C - OPEN(UNIT=IUO2, FILE=OUTFILE2, STATUS='UNKNOWN') -C -C Loop over the energies -C - DO JE=1,NE - IF(NE.GT.1) THEN - ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1) - ELSEIF(NE.EQ.1) THEN - ECIN=E0 - ENDIF - CALL LPM(ECIN,XLPM,*6) - XLPM1=XLPM/A - IF(IPRINT.GT.0) WRITE(IUO1,56) A,XLPM1 - IF(ITL.EQ.0) THEN - VK(JE)=SQRT(ECIN+VINT)*CONV*A*ONEC - VK2(JE)=CABS(VK(JE)*VK(JE)) - ENDIF - GAMMA=1./(2.*XLPM1) - IF(IPOTC.EQ.0) THEN - VK(JE)=VK(JE)+IC*GAMMA - ENDIF - IF(I_MFP.EQ.0) THEN - VK(JE)=CMPLX(REAL(VK(JE))) - VK2(JE)=CABS(VK(JE)*VK(JE)) - ENDIF - IF(I_VIB.EQ.1) THEN - IF(IDCM.GE.1) WRITE(IUO1,22) - DO JAT=1,N_PROT - IF(IDCM.EQ.0) THEN - XK2UJ2=VK2(JE)*UJ2(JAT) - ELSE - XK2UJ2=VK2(JE)*UJ_SQ(JAT) - WRITE(IUO1,23) JAT,UJ_SQ(JAT)*A*A - ENDIF - CALL DWSPH(JAT,JE,XK2UJ2,TLT,I_VIB) - DO LAT=0,LMAX(JAT,JE) - TL(LAT,1,JAT,JE)=TLT(LAT,1,JAT,JE) - ENDDO - ENDDO - ENDIF -C -C Eigenvalue calculation -C - IF(I_PWM.EQ.0) THEN - CALL EIG_MAT_MS(JE,ECIN) - ELSE - CALL SPEC_RAD_POWER(JE,ECIN) - ENDIF -C -C -C End of the loop on the energy -C - ENDDO - GOTO 7 -C - 6 WRITE(IUO1,55) -C - 22 FORMAT(16X,'INTERNAL CALCULATION OF MEAN SQUARE DISPLACEMENTS',/, - &25X,' BY DEBYE UNCORRELATED MODEL:',/) - 23 FORMAT(21X,'ATOM TYPE ',I5,' MSD = ',F8.6,' ANG**2') - 55 FORMAT(///,12X,' <<<<<<<<<< THIS VALUE OF ILPM IS NOT', - &'AVAILABLE >>>>>>>>>>') - 56 FORMAT(4X,'LATTICE PARAMETER A = ',F6.3,' ANGSTROEMS',4X, - *'MEAN FREE PATH = ',F6.3,' * A',//) -C - 7 RETURN -C - END -C -C======================================================================= -C - SUBROUTINE SPEC_RAD_POWER(JE,E_KIN) -C -C This subroutine stores the G_o T kernel matrix and computes -C its spectral radius using the power method. In addition, -C a scalar convergence acceleration method can be used to -C refine the results. -C -C Author : D. Sebilleau -C -C Last modified : 4 Feb 2013 -C -C INCLUDE 'spec.inc' - USE DIM_MOD - USE CONVTYP_MOD - USE COOR_MOD - USE OUTFILES_MOD - USE OUTUNITS_MOD - USE TRANS_MOD -C -C - PARAMETER(NLTWO=2*NL_M) -C - INTEGER NN(0:N_ORD_M) -C - REAL R_SN(N_ORD_M),R_SN1(N_ORD_M),R_SN2(N_ORD_M) -C - COMPLEX*16 TLK,SUM_L,EXPKJ,ONEC,MULT - COMPLEX*16 YLM(0:NLTWO,-NLTWO:NLTWO) - COMPLEX*16 X(LINMAX*NATCLU_M),Y(LINMAX*NATCLU_M) - COMPLEX*16 AX(LINMAX*NATCLU_M),AY(LINMAX*NATCLU_M),SUM_AX,SUM_AY - COMPLEX*16 HL1(0:NLTWO),A(LINMAX*NATCLU_M,LINMAX*NATCLU_M) - COMPLEX*16 ZEROC,IC,JC,JMULT - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M) -C - REAL*8 PI,ATTKJ,GNT(0:N_GAUNT),XKJ,YKJ,ZKJ,RKJ,ZDKJ,KRKJ - REAL*8 NORMX2,NORMY2,NORMAX2,NORMAY2,RLIN -C -C - CHARACTER*10 METH - CHARACTER*24 OUTFILE,PATH -C - DATA PI /3.1415926535898D0/ -C - ONEC=(1.D0,0.D0) - IC=(0.0D0,1.0D0) - ZEROC=(0.D0,0.D0) - IBESS=3 - N_CALL=0 - N_ACC=0 - K_ACC=0 - I_CONV=0 -C - IF(N_MAX.GT.N_ORD_M) THEN - WRITE(IUO1,99) N_MAX - STOP - ENDIF -C -C Starting vector for power method: -C -C (JC^0,JC^1,JC^2,JC^3, ... ,JC^N)^EXPO -C - JC=DCOS(PI/6.D0)+IC*DSIN(PI/6.D0) - JMULT=JC**EXPO -C - WRITE(IUO1,6) -C -C Construction of the the multiple scattering kernel matrix G_o T -C -C Elements are stored in A(KLIN,JLIN) using linear indices -C JLIN and KLIN representing (J,LJ) and (K,LK) -C - JLIN=0 - DO JTYP=1,N_PROT - NBTYPJ=NATYP(JTYP) - LMJ=LMAX(JTYP,JE) - DO JNUM=1,NBTYPJ - JATL=NCORR(JNUM,JTYP) - XJ=SYM_AT(1,JATL) - YJ=SYM_AT(2,JATL) - ZJ=SYM_AT(3,JATL) -C - DO LJ=0,LMJ - DO MJ=-LJ,LJ - JLIN=JLIN+1 -C - KLIN=0 - DO KTYP=1,N_PROT - NBTYPK=NATYP(KTYP) - LMK=LMAX(KTYP,JE) - DO KNUM=1,NBTYPK - KATL=NCORR(KNUM,KTYP) - IF(KATL.NE.JATL) THEN - XKJ=DBLE(SYM_AT(1,KATL)-XJ) - YKJ=DBLE(SYM_AT(2,KATL)-YJ) - ZKJ=DBLE(SYM_AT(3,KATL)-ZJ) - RKJ=DSQRT(XKJ*XKJ+YKJ*YKJ+ZKJ*ZKJ) - KRKJ=DBLE(VK(JE))*RKJ - ATTKJ=DEXP(-DIMAG(DCMPLX(VK(JE)))* - & RKJ) - EXPKJ=(XKJ+IC*YKJ)/RKJ - ZDKJ=ZKJ/RKJ - CALL SPH_HAR2(2*NL_M,ZDKJ,EXPKJ,YLM, - & LMJ+LMK) - CALL BESPHE2(LMJ+LMK+1,IBESS,KRKJ, - & HL1) - ENDIF -C - DO LK=0,LMK - L_MIN=ABS(LK-LJ) - L_MAX=LK+LJ - TLK=DCMPLX(TL(LK,1,KTYP,JE)) - DO MK=-LK,LK - KLIN=KLIN+1 - A(KLIN,JLIN)=ZEROC - SUM_L=ZEROC - IF(KATL.NE.JATL) THEN - CALL GAUNT2(LK,MK,LJ,MJ,GNT) -C - DO L=L_MIN,L_MAX,2 - M=MJ-MK - IF(ABS(M).LE.L) THEN - SUM_L=SUM_L+(IC**L)* - & HL1(L)*YLM(L,M)*GNT(L) - ENDIF - ENDDO - SUM_L=SUM_L*ATTKJ*4.D0*PI*IC - ELSE - SUM_L=ZEROC - ENDIF -C - IF(KATL.NE.JATL) THEN - A(KLIN,JLIN)=TLK*SUM_L - ENDIF -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - ENDDO - ENDDO -C - NLIN=JLIN -C -C Power method approximation of the spectral radius : -C -C SR(A) = lim ||A x||^n / ||x||^n -C n --> +infty -C for any starting vector x -C - RLIN=DFLOAT(NLIN) -C -C Initialization the vectors and their squared norm -C - MULT=JMULT -C - AX(1)=ONEC - AY(1)=ONEC -C - NORMAX2=1.D0 - NORMAY2=1.D0 -C - DO ILIN=2,NLIN -C - AX(ILIN)=AX(ILIN-1)*MULT - AY(ILIN)=AY(ILIN-1)*MULT - NORMAX2=NORMAX2+DREAL(AX(ILIN)*DCONJG(AX(ILIN))) - NORMAY2=NORMAY2+DREAL(AY(ILIN)*DCONJG(AY(ILIN))) -C - ENDDO -C -C Computation of the vectors and their squared norm -C -C X ---> A*X -C - 120 IF(N_CALL.EQ.0) THEN - J_INI=1 - ELSE - J_INI=N_CALL*N_ITER+1 - ENDIF - J_FIN=MIN((N_CALL+1)*N_ITER,N_MAX) -C - IF(J_INI.GT.J_FIN) GOTO 112 -C - DO JORD=J_INI,J_FIN -C - NORMX2=NORMAX2 - NORMY2=NORMAY2 -C - DO JLIN=1,NLIN -C - X(JLIN)=AX(JLIN) - Y(JLIN)=AY(JLIN) -C - ENDDO -C - NORMAX2=0.D0 - NORMAY2=0.D0 -C - DO ILIN=1,NLIN -C - SUM_AX=ZEROC - SUM_AY=ZEROC -C - DO JLIN=1,NLIN -C - SUM_AX=SUM_AX+A(ILIN,JLIN)*X(JLIN) - SUM_AY=SUM_AY+A(JLIN,ILIN)*Y(JLIN) -C - ENDDO -C - AX(ILIN)=SUM_AX - AY(ILIN)=SUM_AY - NORMAX2=NORMAX2+DREAL(SUM_AX*DCONJG(SUM_AX)) - NORMAY2=NORMAY2+DREAL(SUM_AY*DCONJG(SUM_AY)) -C - ENDDO -C - R_SN1(JORD)=SQRT(REAL(NORMAX2/NORMX2)) - R_SN2(JORD)=SQRT(REAL(NORMAY2/NORMY2)) - R_SN(JORD)=SQRT(R_SN1(JORD)*R_SN2(JORD)) -C - ENDDO -C - IF(J_INI.EQ.1) THEN - WRITE(IUO1,10) - WRITE(IUO1,46) NLIN - WRITE(IUO1,10) - WRITE(IUO1,61) - ENDIF -C - WRITE(IUO1,10) - WRITE(IUO1,44) R_SN(JORD-1) - WRITE(IUO1,47) JORD-1 -C - 112 IF(I_ACC.GE.1) THEN -C -C Convergence acceleration on the results whenever necessary -C - CALL ACC_CONV(N_CALL,J_INI,J_FIN,N_ACC,K_ACC,K_MAX,I_CONV, - & R_SN,SN,METH) - IF((I_CONV.EQ.1).AND.(I_ACC.EQ.1)) GOTO 111 - IF((I_CONV.EQ.0).OR.(I_ACC.EQ.2)) GOTO 120 -C - 111 WRITE(IUO1,10) - WRITE(IUO1,61) - WRITE(IUO1,10) - WRITE(IUO1,52) METH - WRITE(IUO1,48) CDABS(SN(N_ACC,K_ACC)) - WRITE(IUO1,49) N_ACC - WRITE(IUO1,51) K_ACC -C - IF(I_ACC.EQ.2) THEN - OPEN(UNIT=35, FILE='div/n_k_table.lis', STATUS='unknown') - DO L=0,K_MAX - NN(L)=L - ENDDO - WRITE(35,*) ' ' - WRITE(35,*) ' (N,K) TABLE OF ',METH,' - & METHOD' - WRITE(35,*) ' ' - WRITE(35,198) (NN(K),K=0,K_MAX) - WRITE(35,199) - WRITE(35,*) ' ' - DO N=0,K_MAX - WRITE(35,200) N,(CDABS(SN(N,K)), K=0,K_MAX-N) - WRITE(35,*) ' ' - ENDDO - ENDIF -C - WRITE(IUO1,10) - WRITE(IUO1,60) - IF(I_ACC.EQ.2) THEN - WRITE(IUO1,210) - ENDIF -C - ENDIF -C - WRITE(IUO2,*) E_KIN,R_SN(JORD-1) -C - RETURN -C - 5 FORMAT(/,11X,'----------------- EIGENVALUE ANALYSIS ','--------- - &--------') - 6 FORMAT(/,11X,'---------------- POWER METHOD ANALYSIS ','-------- - &--------') - 10 FORMAT(11X,'-',54X,'-') - 15 FORMAT(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-') - 20 FORMAT(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-') - 25 FORMAT(11X,'-',6X,'1 EIGENVALUE IS > 1 ON A TOTAL OF ',I8,6X,'-') - 30 FORMAT(11X,'-',4X,I5,' EIGENVALUES ARE > 1 ON A TOTAL OF ',I8,2X, - &'-') - 35 FORMAT(11X,'-',11X,'THE ',I3,' LARGER EIGENVALUES ARE :',11X,'-') - 40 FORMAT(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-') - 44 FORMAT(11X,'-',4X,'SPECTRAL RADIUS BY THE POWER METHOD : ',F8.5, - &3X,'-') - 45 FORMAT(11X,'-',4X,'SPECTRAL RADIUS OF THE KERNEL MATRIX : ',F8.5, - &3X,'-') - 46 FORMAT(11X,'-',4X,'MATRIX SIZE : ',I8, - &3X,'-') - 47 FORMAT(11X,'-',4X,'POWER METHOD TRUNCATION ORDER : ',I8, - &3X,'-') - 48 FORMAT(11X,'-',4X,'SPECTRAL RADIUS BY ACCELERATION : ',F8.5, - &3X,'-') - 49 FORMAT(11X,'-',4X,'ACCELERATION TRUNCATION ORDER N : ',I8, - &3X,'-') - 50 FORMAT(11X,'-',4X,'---> THE MULTIPLE SCATTERING SERIES ', - &'CONVERGES ',4X,'-') - 51 FORMAT(11X,'-',4X,'ACCELERATION TRUNCATION ORDER K : ',I8, - &3X,'-') - 52 FORMAT(11X,'-',4X,'ACCELERATION METHOD : ',A10, - &1X,'-') - 55 FORMAT(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/ - &,11X,'-',18X,'SCATTERING SERIES',19X,'-') - 60 FORMAT(11X,'----------------------------------------','---------- - &------',/) - 61 FORMAT(11X,'----------------------------------------','---------- - &------') - 65 FORMAT(11X,'-',4X,' LABEL OF LARGEST EIGENVALUE : ',I5,8X, - &'-') - 70 FORMAT(11X,'-',4X,' LARGEST EIGENVALUE : ','(',F6.3,',',F6. - &3,')',8X,'-') - 75 FORMAT(' ') - 80 FORMAT(' KINETIC ENERGY : ',F7.2,' eV') - 85 FORMAT(' LARGEST MODULUS OF EIGENVALUE : ',F6.3) - 90 FORMAT(' LABEL OF LARGEST EIGENVALUE : ',I5) - 95 FORMAT(' LARGEST EIGENVALUE : (',F6.3,',',F6.3,') - &') - 99 FORMAT(///,12X,' <<<<<<<<<< DIMENSION OF N_ORD_M IN INCLUDE',' - &FILE >>>>>>>>>>',/,12X,' <<<<<<<<<< SHOULD BE AT LEAST ', - &I5, 6X,' >>>>>>>>>>',///) - 100 FORMAT(5X,F7.3,2X,F7.3,2X,F6.3) - 105 FORMAT(7X,'EIGENVALUES :',3X,'MODULUS :') - 110 FORMAT(2X,'-------------------------------') - 198 FORMAT(' K',50(I3,4X)) - 199 FORMAT(' N') - 200 FORMAT(I3,50(2X,F5.3)) - 210 FORMAT(//,' ---> THE (N,K) TABLE HAS BEEN WRITTEN ', - &'IN /div/n_k_table.lis',//) -C - END -C -C======================================================================= -C - SUBROUTINE ACC_CONV(N_CALL,J_INI,J_FIN,N_ACC,K_ACC,K_MAX,I_CONV, - &R_SN,SN,METH) -C -C Use of the various acceleration scheme on a scalar series -C -C AITK : Aitken, modified Aitken -C RICH : Richardson -C SALZ : Salzer -C EPSI : epsilon -C EPSG : generalized epsilon -C RHOA : rho -C THET : theta -C LEGE : Legendre-Toeplitz -C CHEB : Chebyshev-Toeplitz -C OVER : Overholt -C DURB : Durbin -C DLEV : Levin d-transform !! differ by the -C TLEV : Levin t-transform !! choice of w_n -C ULEV : Levin u-transform !! the remainder -C VLEV : Levin v-transform !! estimate given -C ELEV : Levin e-transform !! by I_WN -C EULE : generalized Euler transform -C GBWT : Germain-Bronne-Wimp transform -C VARI : various algorithms : Weniger, BDG, iterated rho -C ITHE : iterated theta : Lubkin, iterated theta, modified Aitken 2 -C EALG : E-algorithm -C -C SN : series to be accelerated -C XN : auxilliary series (chosen with I_XN) : interpolation points -C GN : auxilliary series (chosen with I_GN) ---> E algorithm -C -C N_DIV : number of SN(N+M,L), M > 0 and L < K+1, necessary to compute S(N,K+1) -C example: iterated Aitken method, we need SN(N,K), SN(N+1,K) and SN(N+2,K) -C so N_DIV=2 -C Because of this only (N_TABLE/N_DIV,N_TABLE/N_DIV) (n,k) tables are -C meaningful. Hence the K_MAX size -C -C COL_M : type of columns meaningful in the (n,k) table (example: even columns -C for the epsilon algorithm) -C -C -C Author : D. Sebilleau -C -C Last modified : 1 Mar 2013 -C -C -C INCLUDE 'spec.inc' - USE DIM_MOD - USE CONVACC_MOD L => LEVIN - USE CONVTYP_MOD -C - PARAMETER (N_METH=24) -C - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M) -C - COMPLEX*16 ZEROC,ONEC -C - REAL*8 R2_SN(N_ORD_M),RHO -C - REAL R_SN(N_ORD_M) -C - INTEGER N_DIV(N_METH) -C - CHARACTER*3 COL_M(N_METH) - CHARACTER*4 SCHE(N_METH) - CHARACTER*10 NAME(N_METH),METH -C -C - DATA SCHE /'AITK','RICH','EPSI','RHOA','THET','LEGE','CHEB', - &'OVER','DURB','DLEV','TLEV','ULEV','VLEV','ELEV','EULE','GBWT', - &'EALG','SALZ','VARI','ITHE','EPSG','NONE','NONE','NONE'/ - DATA NAME /' AITKEN ','RICHARDSON',' EPSILON ',' RHO ',' - & THETA ',' LEGENDRE ','CHEBYSHEV ',' OVERHOLT ',' DURBIN ',' - &D LEVIN ',' T LEVIN ',' U LEVIN ',' V LEVIN ',' E LEVIN ',' - & EULER ',' GBW ',' E ',' SALZER ',' VARIA ', - &'ITER THETA',' EPSILON G',' ',' ',' '/ - DATA COL_M /'ALL','ALL','EVE','EVE','EVE','ALL','ALL','ALL', - &'ALL','ALL','ALL','ALL','ALL','ALL','ALL','ALL','ALL','ALL', - &'ALL','ALL','EVE','ALL','ALL','ALL'/ - DATA N_DIV /2,1,1,2,2,1,1,1,4,1,1,1,1,1,1,1,1,1,2,4,1,1,1,1/ -C - J_NAME=0 - I_COL=0 -C - ZEROC=(0.D0,0.D0) - ONEC=(1.D0,0.D0) -C - DO J=J_INI,J_FIN - R2_SN(J)=DBLE(R_SN(J)) - ENDDO -C -C Finding the name of the method -C - DO JM=1,N_METH -C - IF(METHOD.EQ.SCHE(JM)) THEN - J_NAME=JM - K_MAX=N_MAX/N_DIV(JM) - N_COEF=N_DIV(JM) - IF(COL_M(JM).EQ.'EVE') I_COL=2 - ENDIF -C - ENDDO -C -C Initialization of array SN -C - DO N=J_INI,J_FIN -C - DO K=-1,J_FIN -C - SN(N,K)=ZEROC -C - ENDDO -C - ENDDO -C -C Initialisation of the schemes : -C -C -- SN(N,0) is the series to be accelerated -- -C - SN(0,-1)=ZEROC - SN(0,0)=ONEC -C - DO N=J_INI,J_FIN -C - SN(N,-1)=ZEROC - SN(N,0)=R2_SN(N) -C - ENDDO -C - CALL ACC_SCAL(J_INI,J_FIN,METHOD,SN) -C -C Check for convergence : all results equal within ACC -C in a N_TABLE x N_TABLE square -C - IF(I_CONV.EQ.0) THEN - CALL CHECK_CONV(J_INI,J_FIN,N_TABLE,I_CONV,N_ACC,K_ACC,K_MAX, - & N_COEF,I_COL,ACC,RHO,SN) - ENDIF -C - IF(METHOD(2:4).EQ.'LEV') THEN - METH=NAME(J_NAME)(1:9)//CHAR(48+I_WN) - ELSE - METH=NAME(J_NAME) - ENDIF -C -C Incrementation of the number of calls to this subroutine -C if convergence has not been achieved -C - IF((I_CONV.EQ.0).OR.(ABS(I_ACC).EQ.2)) THEN - N_CALL=N_CALL+1 - ELSE - GOTO 10 - ENDIF -C -C Printing the results in the check file -C - 15 FORMAT(2X,'*',3X,I3,5X,'(',E12.6,',',E12.6,')',3X,'(',E12.6,',', - &E12.6,')',5X,'*') - 16 FORMAT(2X,'*',3X,I3,5X,'(',E12.6,',',E12.6,')',35X,'*') - 17 FORMAT(2X,'*',3X,I3,' --->','(',E12.6,',',E12.6,')',35X,'*') - 18 FORMAT(2X,'*',3X,I3,' --->','(',E12.6,',',E12.6,')',3X,'(',E12.6, - &',',E12.6,')',5X,'*') - 19 FORMAT(2X,'*',3X,I3,25X,'(',E12.6,',',E12.6,')',20X,'*') - 21 FORMAT(2X,'*',3X,I3,25X,'(',E12.6,',',E12.6,')',' <--- - &convergence',2X,'*') - 25 FORMAT(2X,'*',3X,I3,5X,'(',E12.6,',',E12.6,')',3X,'(',E12.6,',', - &E12.6,')',3X,' * <--- convergence') - 35 FORMAT(2X,'***************************************','************ - &************************') - 45 FORMAT(2X,'* ',' - & *') - 65 FORMAT(2X,'* Exact result S = (',E12.6,',',E12.6,') - & *') - 75 FORMAT(2X,'* Order Taylor ',20X,A10,' - & *') - 105 FORMAT(2X,'* Convergence ',A4,26X,A4,14X,'*',/2X,'* - & order ',60X,'*') - 133 FORMAT(//,5X,'<<<<<<<<<< THIS METHOD IS NOT IMPLEMENTED ', - &'>>>>>>>>>>',//) -C - 10 RETURN -C - END -C -C======================================================================= -C - SUBROUTINE ACC_SCAL(J_INI,J_FIN,METHOD,SN) -C -C This subroutine computes the scalar convergence acceleration -C for various methods -C -C Note: let us consider that the (n,k) space pattern indicates -C that we need S(n,k),...,S(n+p,k) in order to compute -C S(n,k+1). We call N the maximum value of n for which -C the S(n,0) are known. Then, to compute S(n,k+1), we -C we need to know up to S(n+pk+p,0). This means that -C the value of n is limited to N-pk-p. -C -C Author : D. Sebilleau -C -C Last modified : 14 Mar 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - USE CONVACC_MOD L => LEVIN -C - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M),XN(-1:N_ORD_M) - COMPLEX*16 NK(-1:3*N_ORD_M,-1:N_ORD_M) - COMPLEX*16 SI(-1:N_ORD_M),TN(-1:N_ORD_M,-1:N_ORD_M) - COMPLEX*16 GN(-1:N_ORD_M,-1:N_ORD_M,-1:N_ORD_M) - COMPLEX*16 DELTA1,DELTA2,DELTA3,DELTA4,DELTA5 - COMPLEX*16 NUM,DEN,FIR - COMPLEX*16 VAR,A,B,C,ZEROC,ONEC,XK -C - REAL*8 EPS,EPSS,MUL -C - CHARACTER*4 METHOD -C - DATA EPS,EPSS /1.D-12,1.D-150/ -C - ZEROC=(0.D0,0.D0) - ONEC=(1.D0,0.D0) -C - N_APP=J_FIN - N_INI=0 -C - IF(METHOD.EQ.'AITK') THEN -C -C The iterated Aitken and modified Aitken schemes -C -C H. H. H. Homeier, Num. Alg. 8, 47 (1994) -C -C I_VA = 1 : iterated Aitken (p 49 eq 5) -C I_VA = 2 : modified iterated Aitken (p 49 eq 9) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - IF(I_VA.EQ.1) THEN - XK=ONEC - ELSEIF(I_VA.EQ.2) THEN - XK=ONEC*DFLOAT(K+K+1)/DFLOAT(K+1) - ENDIF -C - DO N=N_INI,N_APP-2*K-2 -C - DELTA1=(SN(N+1,K)-SN(N,K)) - DELTA2=(SN(N+2,K)-SN(N+1,K)) - IF(CDABS(DELTA2).LT.EPSS) THEN - DELTA3=ZEROC - ELSE - DELTA3=DELTA1*DELTA1/(DELTA2-DELTA1) - ENDIF - SN(N,K+1)=SN(N,K)-XK*DELTA3 - N_COUNT=N_COUNT+1 -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'RICH') THEN -C -C The Richardson scheme -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-2 -C - DELTA1=SN(N+1,0)-SN(N,0) - DELTA2=SN(N+K+2,0)-SN(N+K+1,0) - IF(CDABS(DELTA2-DELTA1).LT.EPSS) GOTO 10 - SN(N,K+1)=(SN(N+1,K)*DELTA1-SN(N,K)*DELTA2)/(DELTA1- - & DELTA2) - N_COUNT=N_COUNT+1 -C - 10 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'SALZ') THEN -C -C The Salzer scheme -C -C A. Hautot, http://physinfo.org/Acc_Conv/Acc_Conv_Part4.pdf (p 12) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - DELTA1=DFLOAT(N+K+2) - DELTA2=DFLOAT(N+1) - IF(CDABS(DELTA2-DELTA1).LT.EPSS) GOTO 11 - SN(N,K+1)=(SN(N+1,K)*DELTA1-SN(N,K)*DELTA2)/(DELTA1- - & DELTA2) - N_COUNT=N_COUNT+1 -C - 11 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'GBWT') THEN -C -C The GBW scheme -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C (eq 6.1-5 p33) -C - N_COUNT=0 -C - CALL INTERP_POINTS(I_XN,2*N_APP,ALPHA,BETA,SN,XN) -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - DELTA1=XN(N) - DELTA2=XN(N+K+1) - IF(CDABS(DELTA1-DELTA2).LT.EPSS) GOTO 80 -C - SN(N,K+1)=(DELTA1*SN(N+1,K)-DELTA2*SN(N,K))/(DELTA1- - & DELTA2) - N_COUNT=N_COUNT+1 -C - 80 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'EPSI') THEN -C -C The epsilon scheme -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) (p 21 eq 4-2-1) -C A. Salam, J. Comput. Appl. Math 46, 455 (1993) (p 456) -C - N_COUNT=0 -C - CALL INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - DELTA1=(SN(N+1,K)-SN(N,K)) - IF(CDABS(DELTA1).LT.EPSS) GOTO 20 - IF(I_VA.EQ.1) THEN - VAR=ONEC - ELSEIF(I_VA.EQ.2) THEN - VAR=XN(N+1)-XN(N) - ELSEIF(I_VA.EQ.3) THEN - VAR=XN(N+K+1)-XN(N+K) - ENDIF - SN(N,K+1)=SN(N+1,K-1)+VAR/DELTA1 - N_COUNT=N_COUNT+1 -C - 20 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'EPSG') THEN -C -C The generalized epsilon scheme -C -C M. N. Barber and C. J. Hamer, -C J. Austral. Math. Soc. (Series B) 23, 229 (1982) -C -C DELTA1 = 0 implies SN(N,K+1) = 0 (calculation not performed) -C DELTA2 = 0 implies SN(N+1,K) = SN(N,K) -C ---> convergence achieved for (N,K) -C - N_COUNT=0 -C - DO K=0,N_APP,2 -C - IF(I_VA.EQ.1) THEN - SI(K)=ALPHA - SI(K+1)=ALPHA - ELSEIF(I_VA.EQ.2) THEN - SI(K)=ZEROC - SI(K+1)=-ONEC - ENDIF -C - ENDDO -C - VAR=ONEC -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - NK(N,-1)=ZEROC - DELTA2=(SN(N+1,K)-SN(N,K)) - IF(CDABS(DELTA2).LT.EPSS) GOTO 23 - NK(N,K)=SI(K)*NK(N,K-1)+VAR/DELTA2 - DELTA1=(NK(N,K)-NK(N-1,K)) - IF(CDABS(DELTA1).LT.EPSS) GOTO 23 - SN(N,K+1)=SN(N,K)+VAR/DELTA1 - N_COUNT=N_COUNT+1 - 23 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'RHOA') THEN -C -C The rho scheme -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C -C I_VA = 1 : Osada's formulation (p 87 eq 11-2-1) -C I_VA = 2 : Standard rho (p 34 eq 6.2-2) -C I_VA = 3 : (p 37 eq 6.3-3) ??????? -C -C ALPHA is the decay parameter. It is -C generally advised to take it as -C an integer to ensure convergence -C - IF(DREAL(ALPHA).LT.0.0D0) THEN -C -C Drummond approximation to alpha : -C -C (non valid for the Taylor expansion -C of 1/1+x as in this case the denominator -C is always zero) -C - DO N=N_INI,N_APP -C - DELTA1=SN(N+3,1)-SN(N+2,1) - DELTA2=SN(N+2,1)-SN(N+1,1) - DELTA3=SN(N+1,1)-SN(N,1) -C - NUM=(DELTA2-DELTA3)*(DELTA1-DELTA2) - DEN=(DELTA1-DELTA2)*DELTA2-(DELTA2-DELTA3)*DELTA1 - ALPHA=NUM/DEN-ONEC -C - ENDDO -C - ENDIF -C - N_COUNT=0 -C - CALL INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - IF(I_VA.NE.3) THEN - DEN=(SN(N+1,K)-SN(N,K)) - ELSE - DELTA2=(SN(N+2,K)-SN(N+1,K)) - DELTA3=(SN(N+1,K)-SN(N,K)) - DEN=(XN(N+2*K)-XN(N+1))*DELTA3-(XN(N+2*K-1)- - & XN(N))*DELTA2 - ENDIF -C - IF(CDABS(DEN).LT.EPS) GOTO 30 - IF((I_VA.EQ.3).AND.(N.EQ.(N_APP-K-1))) GOTO 30 -C - IF(I_VA.EQ.1) THEN - NUM=(DFLOAT(K+1)+ALPHA) - ELSEIF(I_VA.EQ.2) THEN - NUM=XN(N+K+1)-XN(N) - ELSEIF(I_VA.EQ.3) THEN - NUM=(XN(N+2*K+1)-XN(N))*DELTA2*DELTA3 - ENDIF -C - SN(N,K+1)=SN(N+1,K-1)+NUM/DEN - N_COUNT=N_COUNT+1 -C - 30 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'THET') THEN -C -C The theta scheme -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C E. Calicetti et al, Phys. Rep. 446, 1 (2007) -C -C I_VA = 1 : standard formulation (eq 71 p 19) -C I_VA = 2 : generalized formulation (eq 10.1-10 p 74) -C - KMAX2=(N_APP-2)/2 -C - N_COUNT=0 - CALL INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) -C - DO K=0,KMAX2 -C - IF(I_VA.EQ.3) GOTO 51 -C - DO N=N_INI,N_APP-2*K-1 -C - DEN=SN(N+1,2*K)-SN(N,2*K) - IF(CDABS(DEN).LT.EPSS) GOTO 40 - IF(I_VA.EQ.1) THEN - NUM=1.D0 - ELSEIF(I_VA.EQ.2) THEN - NUM=XN(N+2*K+1)-XN(N) - ENDIF - SN(N,2*K+1)=SN(N+1,2*K-1)+NUM/DEN - N_COUNT=N_COUNT+1 -C - 40 CONTINUE -C - ENDDO -C - DO N=N_INI,N_APP-4*K-2 -C - DELTA2=SN(N+2,2*K)-SN(N+1,2*K) - DELTA3=SN(N+1,2*K+1)-SN(N,2*K+1) - DELTA4=SN(N+2,2*K+1)-SN(N+1,2*K+1) -C - IF(I_VA.EQ.1) THEN - NUM=DELTA2*DELTA4 - DEN=DELTA4-DELTA3 - ELSEIF(I_VA.EQ.2) THEN - NUM=(XN(N+2*K+2)-XN(N))*DELTA2*DELTA4 - DEN=(XN(N+2*K+2)-XN(N+1))*DELTA3-(XN(N+2*K+1)-XN( - & N))*DELTA4 - ENDIF - IF(CDABS(DEN).LT.EPSS) GOTO 50 -C - SN(N,2*K+2)=SN(N+1,2*K)+NUM/DEN -C - N_COUNT=N_COUNT+1 -C - 50 CONTINUE -C - ENDDO - GOTO 52 -C - 51 DO N=N_INI,N_APP-4*K-3 -C - DELTA1=SN(N+1,K)-SN(N,K) - DELTA2=SN(N+2,K)-SN(N+1,K) - DELTA3=SN(N+3,K)-SN(N+2,K) - DELTA4=DELTA3-DELTA2 - DELTA5=DELTA2-DELTA1 - DEN=DELTA3*DELTA5-DELTA1*DELTA4 - IF(CDABS(DEN).LT.EPSS) GOTO 53 -C - NUM=DELTA1*DELTA2*DELTA4 - SN(N,K+1)=SN(N+1,K)-NUM/DEN -C - N_COUNT=N_COUNT+1 -C - 53 CONTINUE -C - ENDDO -C - 52 CONTINUE -C - ENDDO -C - ELSEIF(METHOD.EQ.'LEGE') THEN -C -C The Legendre-Toeplitz scheme -C -C D. A. Smith and W. F. Ford, SIAM J. Numer. Anal. 16, 223 (1979) -C eq. (4.3), (4.4) p. 231 -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DELTA1=DFLOAT(K+K+1)/DFLOAT(K+1) - DELTA2=DFLOAT(K)/DFLOAT(K+1) -C - DO N=N_INI,N_APP-K-1 - SN(N,K+1)=DELTA1*(2.D0*SN(N+1,K)-SN(N,K))-DELTA2*SN( - & N,K-1) - N_COUNT=N_COUNT+1 -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'CHEB') THEN -C -C The Chebyshev-Toeplitz scheme -C -C D. A. Smith and W. F. Ford, SIAM J. Numer. Anal. 16, 223 (1979) -C eq. (4.4) - (4.9) p. 231 -C - SI(0)=ONEC - SI(1)=3.D0*ONEC -C - DO N=N_INI,N_APP-1 -C - TN(N,0)=SN(N,0) - TN(N,1)=SN(N,0)+2.D0*TN(N+1,0) - SN(N,0)=TN(N,0)/SI(0) - SN(N,1)=TN(N,1)/SI(1) -C - ENDDO -C - N_COUNT=0 -C - DO K=1,N_APP-1 -C - SI(K+1)=6.D0*SI(K)-SI(K-1) -C - DO N=N_INI,N_APP-K-1 -C - TN(N,K+1)=2.D0*TN(N,K)+4.D0*TN(N+1,K)-TN(N,K-1) - SN(N,K+1)=TN(N,K+1)/SI(K+1) - N_COUNT=N_COUNT+1 -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'OVER') THEN -C -C The Overholt scheme -C -C H. H. H. Homeier, J. Mol. Struc. (Theochem) 368, 81 (1996) -C (eq 36 p 85) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-2 -C - DELTA1=ONEC - DELTA2=ONEC - MUL=1.D0 -C - DO I=1,K+1 -C - DELTA1=DELTA1*(SN(N+K+1,0)-SN(N+K,0)) - DELTA2=DELTA2*(SN(N+K+2,0)-SN(N+K+1,0)) - MUL=MUL*EPS -C - ENDDO -C - IF(CDABS(DELTA1-DELTA2).LT.MUL) GOTO 60 -C - SN(N,K+1)=(DELTA1*SN(N+1,K)-DELTA2*SN(N,K))/(DELTA1- - & DELTA2) - N_COUNT=N_COUNT+1 -C - 60 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'DURB') THEN -C -C The Durbin scheme -C -C C. Brezinski and M. Redivo Zaglia, Comput. Appl. Math. 26, 171 (2007) -C (eq 25 p 185) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-4*K-4 -C - DELTA1=3.D0*SN(N+2,K)*SN(N+2,K)-4.D0*SN(N+1,K)*SN(N+ - & 3,K)+SN(N,K)*SN(N+4,K) - DELTA2=6.D0*SN(N+2,K)-4.D0*(SN(N+1,K)+SN(N+3,K))+SN( - & N,K)+SN(N+4,K) - IF(CDABS(DELTA2).LT.EPSS) GOTO 70 - SN(N,K+1)=DELTA1/DELTA2 - N_COUNT=N_COUNT+1 -C - 70 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD(2:4).EQ.'LEV') THEN -C -C The Levin schemes -C - CALL LEVIN(ALPHA,BETA,SN,I_WN,I_XN,N_APP,N_INI,L,N_COUNT, - & METHOD(1:1)) -C - ELSEIF(METHOD.EQ.'EULE') THEN -C -C The generalized Euler scheme -C -C D. A. Smith and W. F. Ford, SIAM J. Numer. Anal. 16, 223 (1979) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-K-1 -C - SN(N,K+1)=(SN(N+1,K)-BETA*SN(N,K))/(ONEC-BETA) - N_COUNT=N_COUNT+1 -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'VARI') THEN -C -C Various schemes -C -C I_VA = 1 : Weniger lambda scheme (p 87 eq 11-2-1) -C I_VA = 2 : Weniger sigma scheme (p 87 eq 11-2-2) -C I_VA = 3 : Weniger mu scheme (p 87 eq 11-2-3) -C I_VA = 4 : iterated rho scheme -C I_VA = 5 : Bjorstad, Dahlquist and Grosse scheme (p 77 eq 6-3-4) -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C - N_COUNT=0 -C - CALL INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-2*K-2 -C - IF(I_VA.EQ.1) THEN - A=BETA+DFLOAT(N) - B=A+ONEC - C=A - ELSEIF(I_VA.EQ.2) THEN - A=ALPHA+DFLOAT(N+K) - B=A+ONEC - C=A - ELSEIF(I_VA.EQ.3) THEN - A=ALPHA+DFLOAT(N-K) - B=A+ONEC - C=A - ELSEIF(I_VA.EQ.4) THEN - A=XN(N+2*K+2)-XN(N) - B=XN(N+2*K+1)-XN(N) - C=XN(N+2*K+2)-XN(N+1) - ELSEIF(I_VA.EQ.5) THEN - A=(DFLOAT(2*K+1)+ALPHA)/(DFLOAT(2*K)+ALPHA) - B=ONEC - C=ONEC - ENDIF -C - DELTA1=SN(N+1,K)-SN(N,K) - DELTA2=SN(N+2,K)-SN(N+1,K) - DEN=B*DELTA2-C*DELTA1 - IF(CDABS(DEN).LT.EPSS) GOTO 81 - SN(N,K+1)=SN(N+1,K)-A*DELTA1*DELTA2/DEN - N_COUNT=N_COUNT+1 -C - 81 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'ITHE') THEN -C -C Iterated theta schemes -C -C R. Thukral, Appl. Math. Comput. 187, 1502 (2007) <--- Lubkin -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C -C I_VA = 1 : Lubkin transform = iterated theta transform (p 79 eq 10-3-6) -C I_VA = 2 : iterated theta transform (p 84 eq 11-1-5) -C I_VA = 3 : 2nd modified iterated Aitken transform (p 85 eq 11-1-12) -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,N_APP-3*K-3 -C - DELTA1=SN(N+1,K)-SN(N,K) - DELTA2=SN(N+2,K)-SN(N+1,K) - DELTA3=SN(N+3,K)-SN(N+2,K) - IF(I_VA.EQ.1) THEN - NUM=-DELTA1*DELTA2*(DELTA3-DELTA2) - DEN=DELTA3*(DELTA2-DELTA1)-DELTA1*(DELTA3-DELTA2) - FIR=SN(N+1,K) - ELSEIF(I_VA.EQ.2) THEN - NUM=DELTA1*DELTA1*DELTA1*(DELTA3-DELTA2) - DEN=DELTA1*DELTA1*(DELTA3-DELTA2)-DELTA2*DELTA2*( - & DELTA2-DELTA1) - FIR=SN(N+1,K) - ELSEIF(I_VA.EQ.3) THEN - NUM=DELTA2*DELTA2*DELTA3*(DELTA3-DELTA2) - DEN=DELTA2*DELTA2*(DELTA3-DELTA2)-DELTA3*DELTA3*( - & DELTA2-DELTA1) - FIR=SN(N+2,K) - ENDIF -C - IF(CDABS(DEN).LT.EPSS) GOTO 90 - SN(N,K+1)=FIR+NUM/DEN - N_COUNT=N_COUNT+1 -C - 90 CONTINUE -C - ENDDO -C - ENDDO -C - ELSEIF(METHOD.EQ.'EALG') THEN -C -C E algorithm -C - N_COUNT=0 -C - CALL INTERP_POINTS(I_XN,2*N_APP,ALPHA,BETA,SN,XN) - CALL REMAIN_SERIES(I_GN,N_APP,SN,XN,GN) -C -C Computing the GN(N,K,I) -C - DO K=1,N_APP -C - DO N=N_INI,N_APP -C - DEN=GN(N+1,K-1,K)-GN(N,K-1,K) - IF(CDABS(DEN).LT.EPSS) GOTO 91 - DELTA1=GN(N,K-1,K)/DEN -C - DO I=0,N_APP -C - NUM=GN(N+1,K-1,I)-GN(N,K-1,I) - GN(N,K,I)=GN(N,K-1,I)-NUM*DELTA1 -C - ENDDO -C - 91 CONTINUE -C - ENDDO -C - ENDDO -C -C Computing SN(N,K) -C - DO K=1,N_APP -C - DO N=N_INI,N_APP-K -C - DEN=GN(N+1,K-1,K)-GN(N,K-1,K) - IF(CDABS(DEN).LT.EPSS) GOTO 92 - DELTA1=GN(N,K-1,K)/DEN - SN(N,K)=SN(N,K-1)-(SN(N+1,K-1)-SN(N,K-1))*DELTA1 -C - 92 CONTINUE -C - ENDDO -C - ENDDO -C - ENDIF -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE COEFFICIENTS(N_APP,COEF,BETA,I_SWITCH) -C -C Binomial coefficients ---> ISWITCH = 1 -C Power coefficients ---> ISWITCH = 2 -C Pochhammer coefficients ---> ISWITCH = 3 -C Pochhammer coefficients ---> ISWITCH = 4 (negative argument) -C -C -C Author : D. Sebilleau -C -C Last modified : 28 Sep 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - COMPLEX*16 COEF(-1:3*N_ORD_M,-1:N_ORD_M),BETA,ONEC -C - ONEC=(1.D0,0.D0) -C - COEF(0,0)=ONEC - COEF(1,0)=ONEC - COEF(1,1)=ONEC -C - IF(I_SWITCH.EQ.1) THEN -C - DO N=2,N_APP-1 -C - COEF(N,0)=ONEC - COEF(N,1)=DCMPLX(N) -C - DO K=2,N -C - COEF(N,K)=COEF(N-1,K-1)+COEF(N-1,K) -C - ENDDO -C - ENDDO -C - ELSEIF(I_SWITCH.EQ.2) THEN -C - DO N=0,N_APP -C - COEF(N,0)=ONEC - COEF(N,1)=BETA+DFLOAT(N) -C - DO K=1,N_APP -C - COEF(N,K)=COEF(N,K-1)*(BETA+DFLOAT(N)) -C - ENDDO -C - ENDDO -C - ELSEIF(I_SWITCH.EQ.3) THEN -C - DO N=0,N_APP -C - COEF(N,0)=ONEC -C - DO K=1,N_APP -C - COEF(N,K)=COEF(N,K-1)*(BETA+DFLOAT(N+K-1)) -C - ENDDO -C - ENDDO -C - ELSEIF(I_SWITCH.EQ.4) THEN -C - DO N=0,N_APP -C - COEF(N,0)=ONEC -C - DO K=1,N_APP -C - COEF(N,K)=COEF(N,K-1)*(-BETA-DFLOAT(N+K-1)) -C - ENDDO -C - ENDDO -C - ENDIF -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE LEVIN(ALPHA,BETA,SN,I_WN,I_XN,N_APP,N_INI,L,N_COUNT, - &OMEGA) -C -C This subroutine computes Levin-type transforms of -C a given sequence SN(N,1) -C -C E. J. Weniger, Comp. Phys. Rep. 10, 189 (1989) -C -C Author : D. Sebilleau -C -C Last modified : 28 Feb 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M),AN(-1:N_ORD_M),WN(-1: - &N_ORD_M) - COMPLEX*16 COEF(-1:3*N_ORD_M,-1:N_ORD_M),ZEROC,ONEC - COMPLEX*16 COEL(-1:3*N_ORD_M,-1:N_ORD_M),XN(-1:N_ORD_M) - COMPLEX*16 DEN(-1:N_ORD_M,-1:N_ORD_M),NUM(-1:N_ORD_M,-1:N_ORD_M) - COMPLEX*16 X,XK,XXK,ALPHA,BETA,INC,SUM_1,SUM_2 -C - CHARACTER*1 OMEGA -C - ZEROC=(0.D0,0.D0) - ONEC=(1.D0,0.D0) -C -C Recovery of the original series elements -C - DO N=N_INI,N_APP -C - AN(N)=SN(N,0)-SN(N-1,0) -C - ENDDO -C -C -C Choice of the remainder estimate omega_n -C -C Reference : H. H. H. Homeier, J. Comput. App. Math. 122, 81 (2000) -C -C - IF(OMEGA.EQ.'U') THEN -C - MAX_N=N_APP -C - DO N=N_INI,MAX_N - WN(N)=AN(N)*(DFLOAT(N)+BETA) - ENDDO -C - ELSEIF(OMEGA.EQ.'V') THEN -C - MAX_N=N_APP-1 -C - DO N=N_INI,MAX_N - WN(N)=AN(N)*AN(N+1)/(AN(N)-AN(N+1)) - ENDDO -C - ELSEIF(OMEGA.EQ.'T') THEN -C - MAX_N=N_APP -C - DO N=N_INI,MAX_N - WN(N)=AN(N) - ENDDO -C - ELSEIF(OMEGA.EQ.'D') THEN -C - MAX_N=N_APP-1 -C - DO N=N_INI,MAX_N - WN(N)=AN(N+1) - ENDDO -C - ELSEIF(OMEGA.EQ.'E') THEN -C - MAX_N=N_APP-2 -C - DO N=N_INI,MAX_N - WN(N)=AN(N)*AN(N+2)/AN(N+1) - ENDDO -C - ENDIF -C -C Choice of the Levin-type transform -C -C I_WN = 1 ---> L-type (eq 7.2-8 p 41) (eq 3-13b p 10) -C I_WN = 2 ---> S-type (eq 8.3-7 p 58) (eq 3-14b p 10) -C I_WN = 3 ---> M-type !! BETA > K-1 (eq 9.3-7 p 66) (eq 3-15b p 10) -C I_WN = 4 ---> C-type (eq 3-16c p 11) -C I_WN = 5 ---> Drummond-type (eq 9.5-5 p 70) -C I_WN = 6 ---> R-type (eq 7.14-12 p 47) -C -C References : E. J. Weniger, J. Math. Phys. 45 1209 (2004) -C E. J. Weniger, Comp. Phys. Rep. 10 189 (1989) -C -C -C Check for the possibility of L extension (the standard case is L = 0) -C - IF(L.GT.0) THEN - IF(I_WN.EQ.1) THEN - CALL COEFFICIENTS(MAX_N,COEL,BETA,2) - DO N=N_INI,MAX_N - WN(N)=WN(N)*COEL(N,L) - ENDDO - ELSEIF(I_WN.EQ.2) THEN - CALL COEFFICIENTS(N_APP,COEL,BETA,3) - DO N=N_INI,MAX_N - WN(N)=WN(N)*COEL(N,L) - ENDDO - ELSEIF(I_WN.EQ.3) THEN - CALL COEFFICIENTS(N_APP,COEL,BETA,4) - DO N=N_INI,MAX_N - WN(N)=WN(N)*COEL(N,L) - ENDDO - ENDIF - ENDIF -C - DO N=N_INI,MAX_N -C - X=BETA+DFLOAT(N) - COEF(N,0)=ONEC -C - DO K=1,N_APP -C - XK=DFLOAT(K) - XXK=X+XK -C - IF(I_WN.EQ.1) THEN - INC=XXK/(XXK+ONEC) - COEF(N,K)=X*(INC**K)/XXK - ELSEIF(I_WN.EQ.2) THEN - COEF(N,K)=(XXK-ONEC)*XXK/((XXK+XK-ONEC)*(XXK+XK)) - ELSEIF(I_WN.EQ.3) THEN - COEF(N,K)=(X-XK+ONEC)/(XXK+ONEC) - ELSEIF(I_WN.EQ.4) THEN - SUM_1=ZEROC - SUM_2=ZEROC - DO J=0,K - SUM_1=SUM_1+ALPHA*XXK - SUM_2=SUM_2+ALPHA*(XXK+ONEC) - ENDDO - COEF(N,K)=(ALPHA*X+XK-ONEC)*SUM_1/(ALPHA*XXK*SUM_2) - ELSEIF(I_WN.EQ.5) THEN - COEF(N,K)=ONEC - ELSEIF(I_WN.EQ.6) THEN - CALL INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) - ENDIF -C - ENDDO -C - ENDDO -C -C Computation of the recurrence relation for numerator -C and denominator -C - DO N=N_INI,MAX_N -C -C Starting values -C - NUM(N,0)=SN(N,0)/WN(N) - DEN(N,0)=ONEC/WN(N) - SN(N,0)=NUM(N,0)/DEN(N,0) -C - ENDDO -C - N_COUNT=0 -C - DO K=0,N_APP-1 -C - DO N=N_INI,MAX_N-K-1 -C - IF(I_WN.NE.6) THEN - NUM(N,K+1)=NUM(N+1,K)-COEF(N,K)*NUM(N,K) - DEN(N,K+1)=DEN(N+1,K)-COEF(N,K)*DEN(N,K) - ELSE - NUM(N,K+1)=(NUM(N+1,K)-NUM(N,K))/(XN(N+K+1)-XN(N)) - DEN(N,K+1)=(DEN(N+1,K)-DEN(N,K))/(XN(N+K+1)-XN(N)) - ENDIF - SN(N,K+1)=NUM(N,K+1)/DEN(N,K+1) - N_COUNT=N_COUNT+1 -C - ENDDO -C - ENDDO -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE INTERP_POINTS(I_XN,N_APP,ALPHA,BETA,SN,XN) -C -C This subroutine computes the interpolation -C points used in some algorithms -C -C -C Author : D. Sebilleau -C -C Last modified : 27 Feb 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M),XN(-1:N_ORD_M) - COMPLEX*16 ALPHA,BETA -C - DO N=0,N_APP -C - IF(I_XN.EQ.1) THEN - XN(N)=(SN(N+1,0)-SN(N,0)) - ELSEIF(I_XN.EQ.2) THEN - XN(N)=BETA+DFLOAT(N) - ELSEIF(I_XN.EQ.3) THEN - XN(N)=(BETA+DFLOAT(N))*(SN(N+1,0)-SN(N,0)) - ELSEIF(I_XN.EQ.4) THEN - XN(N)=(BETA+DFLOAT(N))**ALPHA - ELSEIF(I_XN.EQ.5) THEN - XN(N)=((2.D0**N)/BETA)**ALPHA - ENDIF -C - ENDDO -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE REMAIN_SERIES(I_GN,N_APP,SN,XN,GN) -C -C This subroutine computes the G_N series used to describe -C the correction on the remainder estimate -C -C I_GN = 1 : G_{i}(n) = G_{i-1}(n)*x_{n} -C I_GN = 2 : G_{i}(n) = x_{n+i-1} (G transformation) -C I_GN = 3 : G_{i}(n) = DELTA S_{n+i-1} (epsilon alg.) -C I_GN = 4 : G_{i}(n) = G_{i-1}(n)*S_{n}/x_{n} -C I_GN = 5 : G_{i}(n) = G_{i-1}(n)/x_{n} -C I_GN = 6 : G_{i}(n) = (DELTA S_{n})^i (Germain-Bonne) -C I_GN = 7 : G_{i}(n) = x_{n+i-2} (p process) -C -C -C Author : D. Sebilleau -C -C Last modified : 28 Feb 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M),XN(-1:N_ORD_M) - COMPLEX*16 GN(-1:N_ORD_M,-1:N_ORD_M,-1:N_ORD_M) - COMPLEX*16 CNK(-1:3*N_ORD_M,-1:N_ORD_M) - COMPLEX*16 SUM_J,MONE,ONE,ZEROC,ONEC -C - ZEROC=(0.D0,0.D0) - ONEC=(1.D0,0.D0) -C -C Computing binomial coefficients -C - IF(I_GN.EQ.6) CALL COEFFICIENTS(N_APP,CNK,ZEROC,1) -C - DO N=0,N_APP -C - GN(N,0,-1)=ZEROC - GN(N,0,0)=ONEC - MONE=-ONEC -C - DO I=1,N_APP -C - IF(I_GN.EQ.1) THEN - GN(N,0,I)=GN(N,0,I-1)*XN(N) - ELSEIF(I_GN.EQ.2) THEN - GN(N,0,I)=XN(N+I-1) - ELSEIF(I_GN.EQ.3) THEN - GN(N,0,I)=SN(N+I,0)-SN(N+I-1,0) - ELSEIF(I_GN.EQ.4) THEN - GN(N,0,I)=GN(N,0,I-1)*SN(N,0)/XN(N) - ELSEIF(I_GN.EQ.5) THEN - GN(N,0,I)=GN(N,0,I-1)/XN(N) - ELSEIF(I_GN.EQ.6) THEN - MONE=-MONE - SUM_J=ZEROC - ONE=-ONEC - DO J=0,I - ONE=-ONE - SUM_J=SUM_J+ONE*CNK(I,J)*SN(N+J,0) - ENDDO - GN(N,0,I)=MONE*SUM_J - ELSEIF(I_GN.EQ.7) THEN - IF(I.GT.1) THEN - GN(N,0,I)=SN(N+I-1,0)-SN(N+I-2,0) - ELSE - GN(N,0,I)=SN(N,0) - ENDIF - ENDIF -C - ENDDO -C - ENDDO -C - RETURN -C - END -C -C======================================================================= -C - SUBROUTINE CHECK_CONV(J_INI,J_FIN,N_TABLE,I_CONV,N_ACC,K_ACC, - &K_MAX,N_COEF,I_COL,ACC,RHO,SN) -C -C This subroutine checks the convergence of the acceleration process. -C For convergence, there must be a N_TABLE x N_TABLE square in the (n,k) table -C of SN(n,k)that contains equal values (SMALL being the accepted error) -C -C This way, both the column convergence (k fixed and n increasing) -C and the diagonal convergence (n=k increasing) are considered -C -C Important note: in some algorithms, such as the epsilon one, only -C odd columns are meaningful. If I_COL = 0, all -C columns are taken into account while for I_COL = 2 -C only even columns are considered -C -C Author : D. Sebilleau -C -C Last modified : 27 Feb 2011 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - COMPLEX*16 SN(-1:N_ORD_M,-1:N_ORD_M) -C - REAL*8 ACC,REF,COMP,RHO -C -C Step for the columns comparison and change of the size -C of the square when only even columns should be considered -C - IF(I_COL.EQ.0) THEN - I_STEP=1 - N_SQUARE=N_TABLE - N_START=N_TABLE - ELSEIF(I_COL.EQ.2) THEN - I_STEP=2 - N_SQUARE=N_TABLE+N_TABLE-1 - IF(MOD(N_SQUARE,2).EQ.0) THEN - N_START=N_SQUARE - ELSE - N_START=N_SQUARE+1 - ENDIF - ENDIF -C -C Positionning the bottom right corner of the square: (N,K), -C starting from (N_START,N_START). If only even columns -C should be considered, N_START must be even -C - K_FIN=J_FIN/N_COEF - IF(K_FIN.LT.N_START) THEN - I_CONV=0 - GOTO 10 - ENDIF -C - DO K=N_START,K_FIN+N_START-1,I_STEP - IF(K.GT.K_MAX) GOTO 20 - DO N=N_START,K_FIN-K,I_STEP -C -C Checking the values (N1,K1) in the selected square -C - REF=CDABS(SN(N,K)) - I_CONV=1 -C - DO N1=0,N_START,I_STEP - NN=N-N1 - DO K1=0,N_START,I_STEP -C - KK=K-K1 - COMP=CDABS(REF-SN(NN,KK)) - IF(COMP.GT.ACC) THEN - I_CONV=0 - ENDIF -C - ENDDO - ENDDO -C - IF(I_CONV.EQ.1) THEN -C -C All values in the square are equal (within ACC): -C Spectral radius taken as the left top corner value -C - N_ACC=N-N_START+I_STEP - K_ACC=K-N_START+I_STEP - RHO=REF - GOTO 10 - ENDIF -C - ENDDO - 20 CONTINUE - ENDDO -C - 10 RETURN -C - END -C -C======================================================================= -C - SUBROUTINE CONV_SERIES(ACC,N_TABLE,JORD,R_SN,*) -C -C This subroutine checks the convergence of the power method or of -C the Rayleigh quotient method at level JORD without convergence -C acceleration. For convergence, there must be N_TABLE equal values -C (ACC being the accepted error) of SN -C -C Author : D. Sebilleau -C -C Last modified : 27 Feb 2013 -C -C INCLUDE 'spec.inc' -C - USE DIM_MOD - REAL ACC,COMP,R_SN(N_ORD_M) -C - NCONV=0 - DO J=JORD,JORD-N_TABLE,-1 - COMP=ABS(R_SN(J)-R_SN(J-1)) - IF(COMP.LE.ACC) THEN - NCONV=NCONV+1 - ENDIF - ENDDO -C - IF(NCONV.EQ.N_TABLE-1) RETURN 1 -C - RETURN -C - END