From 183e704f34b9a7c2bd1d19d850dd1ea678407050 Mon Sep 17 00:00:00 2001 From: Sylvain Tricot Date: Thu, 13 Mar 2025 11:20:31 +0100 Subject: [PATCH] Boost the maximum scattering order up to 18. The maximum scattering order in the Rehr Albers series expansion was formerly limited to 10. The new limit is now 18. --- src/msspec/calculator.py | 6 +- src/msspec/parameters.py | 10 +- .../fortran/phd_se_noso_nosp_nosym/do_main.f | 2 +- .../phd_se_noso_nosp_nosym/findpaths6.f | 367 ++++++++++++++++++ .../phd_se_noso_nosp_nosym/findpaths7.f | 367 ++++++++++++++++++ .../phd_se_noso_nosp_nosym/findpaths8.f | 367 ++++++++++++++++++ .../phd_se_noso_nosp_nosym/findpaths9.f | 367 ++++++++++++++++++ 7 files changed, 1477 insertions(+), 9 deletions(-) create mode 100644 src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths6.f create mode 100644 src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths7.f create mode 100644 src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths8.f create mode 100644 src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths9.f diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 17f36ba..af2b6f4 100644 --- a/src/msspec/calculator.py +++ b/src/msspec/calculator.py @@ -17,8 +17,8 @@ # along with this msspec. If not, see . # # Source file : src/msspec/calculator.py -# Last modified: Tue, 25 Oct 2022 16:21:38 +0200 -# Committed by : Sylvain Tricot 1666707698 +0200 +# Last modified: Thu, 13 Mar 2025 11:20:31 +0100 +# Committed by : Sylvain Tricot """ @@ -379,7 +379,7 @@ class _MSCALCULATOR(Calculator): 'LI_M' : get_li(self.spec_parameters.extra_level) + 1, 'NEMET_M' : 1, 'NO_ST_M' : self.spec_parameters.calc_no, - 'NDIF_M' : 10, + 'NDIF_M' : 18, 'NSO_M' : 2, 'NTEMP_M' : 1, 'NODES_EX_M' : 3, diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index 8f1e8b3..c8e232b 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -19,7 +19,7 @@ # along with this msspec. If not, see . # # Source file : src/msspec/parameters.py -# Last modified: Thu, 27 Feb 2025 15:47:32 +0100 +# Last modified: Thu, 13 Mar 2025 11:20:31 +0100 # Committed by : Sylvain Tricot @@ -606,7 +606,7 @@ class SpecParameters(BaseParameters): Parameter('eigval_beta', types=float, default=1., fmt='.2f'), Parameter('calc_no', types=int, limits=[0, 8], default=1, fmt='d'), - Parameter('calc_ndif', types=int, limits=[1, 10], default=3, + Parameter('calc_ndif', types=int, limits=[1, 18], default=3, fmt='d'), Parameter('calc_ispher', types=int, limits=[0, 1], default=1, fmt='d'), @@ -640,7 +640,7 @@ class SpecParameters(BaseParameters): fmt='d'), Parameter('calc_ira', types=int, limits=[0, 1], default=0, fmt='d'), Parameter('calc_ipw', types=int, limits=[0, 1], default=0, fmt='d'), - Parameter('calc_ncut', types=int, limits=[0, 10], default=2, + Parameter('calc_ncut', types=int, limits=[0, 18], default=2, fmt='d'), Parameter('calc_pctint', types=float, limits=[1e-4, 999.9999], default=0.01, fmt='.4f'), @@ -1441,7 +1441,7 @@ class CalculationParameters(BaseParameters): It is only meaningful for the series expansion algorithm. Its value is limited to 8 but it is rarely necessary to go beyond 2 or 3."""), - Parameter('scattering_order', types=int, limits=(1, 10), default=3, + Parameter('scattering_order', types=int, limits=(1, 18), default=3, doc=""" The scattering order. Only meaningful for the 'expansion' algorithm. Its value is limited to 10."""), @@ -1509,7 +1509,7 @@ class CalculationParameters(BaseParameters): path is rejected and its contribution to the scattering path operator won’t be computed. """), - Parameter('scattering_order_cutoff', types=int, limits=(0, 10), + Parameter('scattering_order_cutoff', types=int, limits=(0, 18), default=2, doc=""" Used in conjunction with the ‘*plane_wave_normal*’ filter. It states to activate the plane wave approximation (which is fast but diff --git a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/do_main.f b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/do_main.f index d60e4cc..3ac8b23 100644 --- a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/do_main.f +++ b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/do_main.f @@ -1459,7 +1459,7 @@ CST 9 FORMAT(3X,F9.4,1X,F9.4,5X,E12.6,5X,E12.6) 95 FORMAT(////,31X,'AUGER LINE :',A6,//) 97 FORMAT(///,19X,'(PLANE WAVES MULTIPLE SCATTERING - ORDER ',I1,')') & - 98 FORMAT(///,17X,'(SPHERICAL WAVES MULTIPLE SCATTERING - ORDER ',I1, + 98 FORMAT(///,17X,'(SPHERICAL WAVES MULTIPLE SCATTERING - ORDER ',I2, &')') 100 FORMAT(///,8X,'<<<<<<<<<< WRONG NAME FOR THE INITIAL STATE',' >> &>>>>>>>>') diff --git a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths6.f b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths6.f new file mode 100644 index 0000000..bf8635c --- /dev/null +++ b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths6.f @@ -0,0 +1,367 @@ +C +C======================================================================= +C + SUBROUTINE FINDPATHS6(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI, + & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) +C +C This routine generates all the paths and filters them according to the +C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). +C It corresponds to the spin-independent case from a non spin-orbit +C resolved initial core state LI +C +C Last modified : 16 May 2007 +C + USE DIM_MOD +C + USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH + USE COOR_MOD + USE DEBWAL_MOD + USE INIT_L_MOD + USE PATH_MOD + USE ROT_MOD + USE TESTPA_MOD + USE TESTPB_MOD + USE TRANS_MOD + USE TLDW_MOD + USE VARIA_MOD +C + DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M) + DIMENSION JPOS(NDIF_M,3),R(NDIF_M) +C +C +C + COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK + COMPLEX IC,COMPL1,PW(0:NDIF_M) + COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M) + COMPLEX YLM1(0:NL_M,-NL_M:NL_M) + COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2 +C + DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/ +C + IC=(0.,1.) + IEULER=1 +C + IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP)) + IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP)) +C +C I_CP = 0 : all open paths generated +C I_CP = 1 : only closed paths generated +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO JTYP=1,N_TYP + IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP)) + IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP)) + ND=ND+1 +C +C I_ABS = 0 : the atom before the scatterer is not the absorber +C I_ABS = 1 : the atom before the scatterer is the absorber +C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only) +C + IF(ND.EQ.1) THEN + I_ABS=1 + ELSE + I_ABS=0 + ENDIF +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPJ=NATYP(JTYP) + ELSE + NBTYPJ=1 + ENDIF +C + DO JNUM=1,NBTYPJ + JATL=NCORR(JNUM,JTYP) + IF(JATL.EQ.IATL) GOTO 12 + XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL) + YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL) + ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + JPOS(ND,1)=JTYP + JPOS(ND,2)=JNUM + JPOS(ND,3)=JATL + NPATH(ND)=NPATH(ND)+1. + IF(ND.GT.1) THEN + COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/( + &R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(ITYP).EQ.0) THEN + IF(COSTHMIJ.LT.COSFWDI) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(ITYP).EQ.1) THEN + IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + ENDIF + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + RHOIJ=VK(JE)*R(ND) + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1.) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THIJ=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + CALL ARCSIN(COMPL1,CTROIS1,PHIIJ) + IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40 + ZSURFI=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(ITYP)=SIG2(R(ND-1),ITYP) + ENDIF + IF(ABS(ZSURFI).LE.SMALL) THEN + IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN + CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO + &STHMIJ) + ELSE + CSKZ2I=1. + ENDIF + UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.)) + ELSE + UII=UJ2(ITYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UI2=VK2(JE)*UII + CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED) + ENDIF + 40 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ)) + ENDIF + ENDIF + IF(ND.EQ.1) THEN + RHO01=RHOIJ + TH01=THIJ + PHI01=PHIIJ + CALL DJMN2(TH01,RLM01,LF2,2) + GOTO 30 + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + LMJ=LMAX(ITYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42 + CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ, + &BMIJ,CMIJ,RHOMI,RHOIJ) + 30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F + &REF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 42 + I_ABS=0 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO KTYP=1,N_TYP + ND=ND+1 + IF(ND.GT.NDIF) GOTO 20 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPK=NATYP(KTYP) + ELSE + NBTYPK=1 + ENDIF +C + DO KNUM=1,NBTYPK + KATL=NCORR(KNUM,KTYP) + IF(KATL.EQ.JATL) GOTO 22 + JPOS(ND,1)=KTYP + JPOS(ND,2)=KNUM + JPOS(ND,3)=KATL + XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL) + YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL) + ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF(IT(ND-1).EQ.1) GOTO 32 + RHOJK=R(ND)*VK(JE) + NPATH(ND)=NPATH(ND)+1. + COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1)) + &/(R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(JTYP).EQ.0) THEN + IF(COSTHIJK.LT.COSFWDJ) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(JTYP).EQ.1) THEN + IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH + &EN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + IF(IT(ND-1).EQ.1) GOTO 32 + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THJK=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + IF(ND-1.LT.NDIF) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50 + ZSURFJ=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(JTYP)=SIG2(R(ND-1),JTYP) + ENDIF + IF(ABS(ZSURFJ).LE.SMALL) THEN + IF(ABS(COSTHIJK-1.).GT.SMALL) THEN + CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2 + &.*COSTHIJK) + ELSE + CSKZ2J=1. + ENDIF + UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.)) + ELSE + UJJ=UJ2(JTYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UJ2=VK2(JE)*UJJ + CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED) + ENDIF + 50 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK)) + ENDIF + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + LMJ=LMAX(JTYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT)) + & IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 + IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + ENDIF + CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A + &IJK,BIJK,CIJK,RHOIJ,RHOJK) + CEX(ND)=CEXP(IC*RHOJK)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI + &JK,FREF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 32 + CALL FINDPATHS7(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK, + 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) + 32 DIJ=DIJ-R(ND) + 22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDDO + 20 CONTINUE + ND=ND-1 + ENDDO + 42 DIJ=DIJ-R(ND) + 12 IF(ND.GT.1) THEN + IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDIF + ENDDO + ND=ND-1 + ENDDO +C + RETURN +C + END diff --git a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths7.f b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths7.f new file mode 100644 index 0000000..efe07b1 --- /dev/null +++ b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths7.f @@ -0,0 +1,367 @@ +C +C======================================================================= +C + SUBROUTINE FINDPATHS7(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI, + & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) +C +C This routine generates all the paths and filters them according to the +C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). +C It corresponds to the spin-independent case from a non spin-orbit +C resolved initial core state LI +C +C Last modified : 16 May 2007 +C + USE DIM_MOD +C + USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH + USE COOR_MOD + USE DEBWAL_MOD + USE INIT_L_MOD + USE PATH_MOD + USE ROT_MOD + USE TESTPA_MOD + USE TESTPB_MOD + USE TRANS_MOD + USE TLDW_MOD + USE VARIA_MOD +C + DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M) + DIMENSION JPOS(NDIF_M,3),R(NDIF_M) +C +C +C + COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK + COMPLEX IC,COMPL1,PW(0:NDIF_M) + COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M) + COMPLEX YLM1(0:NL_M,-NL_M:NL_M) + COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2 +C + DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/ +C + IC=(0.,1.) + IEULER=1 +C + IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP)) + IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP)) +C +C I_CP = 0 : all open paths generated +C I_CP = 1 : only closed paths generated +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO JTYP=1,N_TYP + IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP)) + IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP)) + ND=ND+1 +C +C I_ABS = 0 : the atom before the scatterer is not the absorber +C I_ABS = 1 : the atom before the scatterer is the absorber +C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only) +C + IF(ND.EQ.1) THEN + I_ABS=1 + ELSE + I_ABS=0 + ENDIF +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPJ=NATYP(JTYP) + ELSE + NBTYPJ=1 + ENDIF +C + DO JNUM=1,NBTYPJ + JATL=NCORR(JNUM,JTYP) + IF(JATL.EQ.IATL) GOTO 12 + XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL) + YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL) + ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + JPOS(ND,1)=JTYP + JPOS(ND,2)=JNUM + JPOS(ND,3)=JATL + NPATH(ND)=NPATH(ND)+1. + IF(ND.GT.1) THEN + COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/( + &R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(ITYP).EQ.0) THEN + IF(COSTHMIJ.LT.COSFWDI) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(ITYP).EQ.1) THEN + IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + ENDIF + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + RHOIJ=VK(JE)*R(ND) + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1.) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THIJ=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + CALL ARCSIN(COMPL1,CTROIS1,PHIIJ) + IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40 + ZSURFI=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(ITYP)=SIG2(R(ND-1),ITYP) + ENDIF + IF(ABS(ZSURFI).LE.SMALL) THEN + IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN + CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO + &STHMIJ) + ELSE + CSKZ2I=1. + ENDIF + UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.)) + ELSE + UII=UJ2(ITYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UI2=VK2(JE)*UII + CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED) + ENDIF + 40 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ)) + ENDIF + ENDIF + IF(ND.EQ.1) THEN + RHO01=RHOIJ + TH01=THIJ + PHI01=PHIIJ + CALL DJMN2(TH01,RLM01,LF2,2) + GOTO 30 + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + LMJ=LMAX(ITYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42 + CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ, + &BMIJ,CMIJ,RHOMI,RHOIJ) + 30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F + &REF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 42 + I_ABS=0 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO KTYP=1,N_TYP + ND=ND+1 + IF(ND.GT.NDIF) GOTO 20 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPK=NATYP(KTYP) + ELSE + NBTYPK=1 + ENDIF +C + DO KNUM=1,NBTYPK + KATL=NCORR(KNUM,KTYP) + IF(KATL.EQ.JATL) GOTO 22 + JPOS(ND,1)=KTYP + JPOS(ND,2)=KNUM + JPOS(ND,3)=KATL + XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL) + YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL) + ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF(IT(ND-1).EQ.1) GOTO 32 + RHOJK=R(ND)*VK(JE) + NPATH(ND)=NPATH(ND)+1. + COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1)) + &/(R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(JTYP).EQ.0) THEN + IF(COSTHIJK.LT.COSFWDJ) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(JTYP).EQ.1) THEN + IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH + &EN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + IF(IT(ND-1).EQ.1) GOTO 32 + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THJK=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + IF(ND-1.LT.NDIF) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50 + ZSURFJ=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(JTYP)=SIG2(R(ND-1),JTYP) + ENDIF + IF(ABS(ZSURFJ).LE.SMALL) THEN + IF(ABS(COSTHIJK-1.).GT.SMALL) THEN + CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2 + &.*COSTHIJK) + ELSE + CSKZ2J=1. + ENDIF + UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.)) + ELSE + UJJ=UJ2(JTYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UJ2=VK2(JE)*UJJ + CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED) + ENDIF + 50 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK)) + ENDIF + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + LMJ=LMAX(JTYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT)) + & IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 + IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + ENDIF + CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A + &IJK,BIJK,CIJK,RHOIJ,RHOJK) + CEX(ND)=CEXP(IC*RHOJK)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI + &JK,FREF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 32 + CALL FINDPATHS8(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK, + 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) + 32 DIJ=DIJ-R(ND) + 22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDDO + 20 CONTINUE + ND=ND-1 + ENDDO + 42 DIJ=DIJ-R(ND) + 12 IF(ND.GT.1) THEN + IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDIF + ENDDO + ND=ND-1 + ENDDO +C + RETURN +C + END diff --git a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths8.f b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths8.f new file mode 100644 index 0000000..0f0e50e --- /dev/null +++ b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths8.f @@ -0,0 +1,367 @@ +C +C======================================================================= +C + SUBROUTINE FINDPATHS8(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI, + & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) +C +C This routine generates all the paths and filters them according to the +C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). +C It corresponds to the spin-independent case from a non spin-orbit +C resolved initial core state LI +C +C Last modified : 16 May 2007 +C + USE DIM_MOD +C + USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH + USE COOR_MOD + USE DEBWAL_MOD + USE INIT_L_MOD + USE PATH_MOD + USE ROT_MOD + USE TESTPA_MOD + USE TESTPB_MOD + USE TRANS_MOD + USE TLDW_MOD + USE VARIA_MOD +C + DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M) + DIMENSION JPOS(NDIF_M,3),R(NDIF_M) +C +C +C + COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK + COMPLEX IC,COMPL1,PW(0:NDIF_M) + COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M) + COMPLEX YLM1(0:NL_M,-NL_M:NL_M) + COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2 +C + DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/ +C + IC=(0.,1.) + IEULER=1 +C + IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP)) + IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP)) +C +C I_CP = 0 : all open paths generated +C I_CP = 1 : only closed paths generated +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO JTYP=1,N_TYP + IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP)) + IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP)) + ND=ND+1 +C +C I_ABS = 0 : the atom before the scatterer is not the absorber +C I_ABS = 1 : the atom before the scatterer is the absorber +C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only) +C + IF(ND.EQ.1) THEN + I_ABS=1 + ELSE + I_ABS=0 + ENDIF +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPJ=NATYP(JTYP) + ELSE + NBTYPJ=1 + ENDIF +C + DO JNUM=1,NBTYPJ + JATL=NCORR(JNUM,JTYP) + IF(JATL.EQ.IATL) GOTO 12 + XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL) + YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL) + ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + JPOS(ND,1)=JTYP + JPOS(ND,2)=JNUM + JPOS(ND,3)=JATL + NPATH(ND)=NPATH(ND)+1. + IF(ND.GT.1) THEN + COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/( + &R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(ITYP).EQ.0) THEN + IF(COSTHMIJ.LT.COSFWDI) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(ITYP).EQ.1) THEN + IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + ENDIF + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + RHOIJ=VK(JE)*R(ND) + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1.) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THIJ=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + CALL ARCSIN(COMPL1,CTROIS1,PHIIJ) + IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40 + ZSURFI=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(ITYP)=SIG2(R(ND-1),ITYP) + ENDIF + IF(ABS(ZSURFI).LE.SMALL) THEN + IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN + CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO + &STHMIJ) + ELSE + CSKZ2I=1. + ENDIF + UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.)) + ELSE + UII=UJ2(ITYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UI2=VK2(JE)*UII + CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED) + ENDIF + 40 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ)) + ENDIF + ENDIF + IF(ND.EQ.1) THEN + RHO01=RHOIJ + TH01=THIJ + PHI01=PHIIJ + CALL DJMN2(TH01,RLM01,LF2,2) + GOTO 30 + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + LMJ=LMAX(ITYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42 + CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ, + &BMIJ,CMIJ,RHOMI,RHOIJ) + 30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F + &REF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 42 + I_ABS=0 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO KTYP=1,N_TYP + ND=ND+1 + IF(ND.GT.NDIF) GOTO 20 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPK=NATYP(KTYP) + ELSE + NBTYPK=1 + ENDIF +C + DO KNUM=1,NBTYPK + KATL=NCORR(KNUM,KTYP) + IF(KATL.EQ.JATL) GOTO 22 + JPOS(ND,1)=KTYP + JPOS(ND,2)=KNUM + JPOS(ND,3)=KATL + XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL) + YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL) + ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF(IT(ND-1).EQ.1) GOTO 32 + RHOJK=R(ND)*VK(JE) + NPATH(ND)=NPATH(ND)+1. + COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1)) + &/(R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(JTYP).EQ.0) THEN + IF(COSTHIJK.LT.COSFWDJ) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(JTYP).EQ.1) THEN + IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH + &EN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + IF(IT(ND-1).EQ.1) GOTO 32 + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THJK=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + IF(ND-1.LT.NDIF) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50 + ZSURFJ=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(JTYP)=SIG2(R(ND-1),JTYP) + ENDIF + IF(ABS(ZSURFJ).LE.SMALL) THEN + IF(ABS(COSTHIJK-1.).GT.SMALL) THEN + CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2 + &.*COSTHIJK) + ELSE + CSKZ2J=1. + ENDIF + UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.)) + ELSE + UJJ=UJ2(JTYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UJ2=VK2(JE)*UJJ + CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED) + ENDIF + 50 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK)) + ENDIF + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + LMJ=LMAX(JTYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT)) + & IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 + IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + ENDIF + CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A + &IJK,BIJK,CIJK,RHOIJ,RHOJK) + CEX(ND)=CEXP(IC*RHOJK)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI + &JK,FREF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 32 + CALL FINDPATHS9(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK, + 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) + 32 DIJ=DIJ-R(ND) + 22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDDO + 20 CONTINUE + ND=ND-1 + ENDDO + 42 DIJ=DIJ-R(ND) + 12 IF(ND.GT.1) THEN + IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDIF + ENDDO + ND=ND-1 + ENDDO +C + RETURN +C + END diff --git a/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths9.f b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths9.f new file mode 100644 index 0000000..0053b14 --- /dev/null +++ b/src/msspec/spec/fortran/phd_se_noso_nosp_nosym/findpaths9.f @@ -0,0 +1,367 @@ +C +C======================================================================= +C + SUBROUTINE FINDPATHS9(ND,ITYP,IATL,I_CP,R,XR,YR,ZR,RHOMI,THMI, + & PHIMI,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) +C +C This routine generates all the paths and filters them according to the +C criteria given in the input data file (IFSPH,IFWD,IPW,ILENGTH). +C It corresponds to the spin-independent case from a non spin-orbit +C resolved initial core state LI +C +C Last modified : 16 May 2007 +C + USE DIM_MOD +C + USE APPROX_MOD , ILE => ILENGTH, RLE => RLENGTH + USE COOR_MOD + USE DEBWAL_MOD + USE INIT_L_MOD + USE PATH_MOD + USE ROT_MOD + USE TESTPA_MOD + USE TESTPB_MOD + USE TRANS_MOD + USE TLDW_MOD + USE VARIA_MOD +C + DIMENSION XR(NDIF_M),YR(NDIF_M),ZR(NDIF_M) + DIMENSION JPOS(NDIF_M,3),R(NDIF_M) +C +C +C + COMPLEX PW1,PWI,FTHETA,RHOMI,RHOIJ,RHOJK + COMPLEX IC,COMPL1,PW(0:NDIF_M) + COMPLEX TAU(LINMAX,LINFMAX,NATCLU_M) + COMPLEX YLM1(0:NL_M,-NL_M:NL_M) + COMPLEX YLM2(0:NL_M,-NL_M:NL_M),CTL,CTL2 +C + DATA XCOMP,PI4,SMALL /1.E-10,12.566371,0.0001/ +C + IC=(0.,1.) + IEULER=1 +C + IF(IFWD.EQ.1) COSFWDI=COS(RTHFWD(ITYP)) + IF(IBWD(ITYP).EQ.1) COSBWDI=COS(RTHBWD(ITYP)) +C +C I_CP = 0 : all open paths generated +C I_CP = 1 : only closed paths generated +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO JTYP=1,N_TYP + IF(IFWD.EQ.1) COSFWDJ=COS(RTHFWD(JTYP)) + IF(IBWD(JTYP).EQ.1) COSBWDJ=COS(RTHBWD(JTYP)) + ND=ND+1 +C +C I_ABS = 0 : the atom before the scatterer is not the absorber +C I_ABS = 1 : the atom before the scatterer is the absorber +C I_ABS = 2 : the atom after the scatterer is the absorber (XAS only) +C + IF(ND.EQ.1) THEN + I_ABS=1 + ELSE + I_ABS=0 + ENDIF +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPJ=NATYP(JTYP) + ELSE + NBTYPJ=1 + ENDIF +C + DO JNUM=1,NBTYPJ + JATL=NCORR(JNUM,JTYP) + IF(JATL.EQ.IATL) GOTO 12 + XR(ND)=SYM_AT(1,JATL)-SYM_AT(1,IATL) + YR(ND)=SYM_AT(2,JATL)-SYM_AT(2,IATL) + ZR(ND)=SYM_AT(3,JATL)-SYM_AT(3,IATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + JPOS(ND,1)=JTYP + JPOS(ND,2)=JNUM + JPOS(ND,3)=JATL + NPATH(ND)=NPATH(ND)+1. + IF(ND.GT.1) THEN + COSTHMIJ=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1))/( + &R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(ITYP).EQ.0) THEN + IF(COSTHMIJ.LT.COSFWDI) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(ITYP).EQ.1) THEN + IF((COSTHMIJ.GT.COSBWDI).AND.(COSTHMIJ.LT.-SMALL)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHMIJ.LT.COSFWDI).AND.(COSTHMIJ.GE.0.)) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + ENDIF + IF((IT(ND-1).EQ.1).AND.(ND.GT.1)) GOTO 42 + RHOIJ=VK(JE)*R(ND) + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1.) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THIJ=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + CALL ARCSIN(COMPL1,CTROIS1,PHIIJ) + IF((ND.GT.1).AND.((ND-1).LT.NDIF)) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 40 + ZSURFI=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(ITYP)=SIG2(R(ND-1),ITYP) + ENDIF + IF(ABS(ZSURFI).LE.SMALL) THEN + IF(ABS(COSTHMIJ-1.).GT.SMALL) THEN + CSKZ2I=(CTROIS1-COS(THMI))*(CTROIS1-COS(THMI))/(2.-2.*CO + &STHMIJ) + ELSE + CSKZ2I=1. + ENDIF + UII=UJ2(ITYP)*(1.+CSKZ2I*(RSJ-1.)) + ELSE + UII=UJ2(ITYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UI2=VK2(JE)*UII + CALL DWSPH(ITYP,JE,XK2UI2,TLT,ISPEED) + ENDIF + 40 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UII*(1.-COSTHMIJ)) + ENDIF + ENDIF + IF(ND.EQ.1) THEN + RHO01=RHOIJ + TH01=THIJ + PHI01=PHIIJ + CALL DJMN2(TH01,RLM01,LF2,2) + GOTO 30 + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHMIJ,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + LMJ=LMAX(ITYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THIJ,PHIIJ,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,JTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT))IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 42 + CALL EULER(THIJ,PHIIJ,THMI,PHIMI,AMIJ,BMIJ,CMIJ,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,ITYP,JTYP,JE,I_ABS,ISPEED,ISPHER,AMIJ, + &BMIJ,CMIJ,RHOMI,RHOIJ) + 30 CEX(ND)=CEXP(IC*RHOIJ)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(JATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOIJ,THIJ,PHIIJ,F + &REF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 42 + I_ABS=0 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF-1)) THEN + N_TYP=N_PROT + ELSE + N_TYP=1 + ENDIF +C + DO KTYP=1,N_TYP + ND=ND+1 + IF(ND.GT.NDIF) GOTO 20 +C + IF((I_CP.EQ.0).OR.(ND.NE.NDIF)) THEN + NBTYPK=NATYP(KTYP) + ELSE + NBTYPK=1 + ENDIF +C + DO KNUM=1,NBTYPK + KATL=NCORR(KNUM,KTYP) + IF(KATL.EQ.JATL) GOTO 22 + JPOS(ND,1)=KTYP + JPOS(ND,2)=KNUM + JPOS(ND,3)=KATL + XR(ND)=SYM_AT(1,KATL)-SYM_AT(1,JATL) + YR(ND)=SYM_AT(2,KATL)-SYM_AT(2,JATL) + ZR(ND)=SYM_AT(3,KATL)-SYM_AT(3,JATL) + R(ND)=SQRT(XR(ND)*XR(ND)+YR(ND)*YR(ND)+ZR(ND)*ZR(ND)) + DIJ=DIJ+R(ND) + IF((ILE.EQ.1).AND.(DIJ.GT.RLE)) IT(ND-1)=1 + IF(IT(ND-1).EQ.1) GOTO 32 + RHOJK=R(ND)*VK(JE) + NPATH(ND)=NPATH(ND)+1. + COSTHIJK=(XR(ND)*XR(ND-1)+YR(ND)*YR(ND-1)+ZR(ND)*ZR(ND-1)) + &/(R(ND)*R(ND-1)) + IF(IFWD.EQ.1) THEN + IF(IBWD(JTYP).EQ.0) THEN + IF(COSTHIJK.LT.COSFWDJ) THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ELSEIF(IBWD(JTYP).EQ.1) THEN + IF((COSTHIJK.GT.COSBWDJ).AND. (COSTHIJK.LT.-SMALL)) TH + &EN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + IF((COSTHIJK.LT.COSFWDJ).AND.(COSTHIJK.GE.0.))THEN + NTHOF=NTHOF+1 + IN(ND-1)=1 + IF(NTHOF.GT.NTHOUT) THEN + IT(ND-1)=1 + ENDIF + ENDIF + ENDIF + ENDIF + IF(IT(ND-1).EQ.1) GOTO 32 + CTROIS1=ZR(ND)/R(ND) + IF(CTROIS1.GT.1) THEN + CTROIS1=1. + ELSEIF(CTROIS1.LT.-1.) THEN + CTROIS1=-1. + ENDIF + THJK=ACOS(CTROIS1) + COMPL1= XR(ND)+IC*YR(ND) + IF(ND-1.LT.NDIF) THEN + IF((IDWSPH.EQ.1).AND.(ISPEED.EQ.1)) GOTO 50 + ZSURFJ=ZSURF-ZR(ND-1) + IF(IDCM.EQ.1) THEN + UJ2(JTYP)=SIG2(R(ND-1),JTYP) + ENDIF + IF(ABS(ZSURFJ).LE.SMALL) THEN + IF(ABS(COSTHIJK-1.).GT.SMALL) THEN + CSKZ2J=(CTROIS1-COS(THIJ))*(CTROIS1-COS(THIJ))/(2.-2 + &.*COSTHIJK) + ELSE + CSKZ2J=1. + ENDIF + UJJ=UJ2(JTYP)*(1.+CSKZ2J*(RSJ-1.)) + ELSE + UJJ=UJ2(JTYP) + ENDIF + IF((ISPEED.EQ.0).AND.(IDWSPH.EQ.1)) THEN + XK2UJ2=VK2(JE)*UJJ + CALL DWSPH(JTYP,JE,XK2UJ2,TLT,ISPEED) + ENDIF + 50 IF(IDWSPH.EQ.1) THEN + DW(ND-1)=1. + ELSE + DW(ND-1)=EXP(-VK2(JE)*UJJ*(1.-COSTHIJK)) + ENDIF + ENDIF + IF(IPW.EQ.1) THEN + CALL FACDIF(COSTHIJK,JPOS(ND-1,1),JE,FTHETA) + PWI=FTHETA*DW(ND-1)/R(ND) + PW(ND)=PW(ND-1)*PWI + CTL2=PI4*PW(ND)*CEX(1)/VK(JE) + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + LMJ=LMAX(JTYP,JE) + IF(ND.GT.NCUT) THEN + IT(ND)=1 + ELSE + IT(ND)=0 + ENDIF + CALL HARSPH2(NL_M,TH01,PHI01,YLM1,LF2) + CALL HARSPH2(NL_M,THJK,PHIJK,YLM2,LMJ) + XMAXT=0. + DO LJ=0,LMJ + CTL=CTL2*TL(LJ,1,KTYP,JE)*YLM2(LJ,0) + DO LF=LF1,LF2,ISTEP_LF + PW1=CTL*YLM1(LF,0)*TL(LF,1,1,JE) + XMAXT=AMAX1(XMAXT,CABS(PW1)) + ENDDO + ENDDO + IF((PCTINT*FREF-XMAXT.LT.-XCOMP).AND.(ND.GT.NCUT)) + & IT(ND)=0 + ENDIF + IF((IT(ND-1).EQ.1).OR.(IT(ND).EQ.1)) GOTO 32 + IF((ND.LT.NDIF).OR.(IPW.EQ.0)) THEN + CALL ARCSIN(COMPL1,CTROIS1,PHIJK) + ENDIF + CALL EULER(THJK,PHIJK,THIJ,PHIIJ,AIJK,BIJK,CIJK,IEULER) + IF((I_CP.EQ.1).AND.(ND.EQ.NDIF)) I_ABS=2 + CALL MATDIF(NO,ND-1,LF2,JTYP,KTYP,JE,I_ABS,ISPEED,ISPHER,A + &IJK,BIJK,CIJK,RHOIJ,RHOJK) + CEX(ND)=CEXP(IC*RHOJK)/R(ND) + CEXDW(ND)=CEX(ND)*DW(ND-1) + IF((IJ.EQ.1).OR.(ND.EQ.NCUT)) THEN + IF((I_CP.EQ.0).OR.(KATL.EQ.1)) THEN + CALL PATHOP(JPOS,ND,JE,I_CP,RHO01,PHI01,RHOJK,THJK,PHI + &JK,FREF,IJ,DIJ,TAU) + NPATH2(ND)=NPATH2(ND)+1. + ENDIF + ENDIF + IF(ND.EQ.NDIF) GOTO 32 +c CALL FINDPATHS(ND,KTYP,KATL,I_CP,R,XR,YR,ZR,RHOJK, +c 1 THJK,PHIJK,ZSURF,JPOS,PW,JE,FREF,DIJ,TAU) + 32 DIJ=DIJ-R(ND) + 22 IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDDO + 20 CONTINUE + ND=ND-1 + ENDDO + 42 DIJ=DIJ-R(ND) + 12 IF(ND.GT.1) THEN + IF(IN(ND-1).EQ.1) NTHOF=NTHOF-1 + IT(ND-1)=0 + IN(ND-1)=0 + ENDIF + ENDDO + ND=ND-1 + ENDDO +C + RETURN +C + END