diff --git a/src/install.sh b/src/install.sh index 169f653..7d18681 100755 --- a/src/install.sh +++ b/src/install.sh @@ -69,7 +69,8 @@ error_exit () { printf " for more informations.\n" printf " Below is an excerpt of the last 10 lines:\n" printf "\n" - tail -v -n 10 "${LOGFILE}" +# tail -v -n 10 "${LOGFILE}" + tail -n 10 "${LOGFILE}" # uninstall if needed pip uninstall -y msspec @@ -91,7 +92,7 @@ get_os () { break fi done - + } @@ -105,7 +106,7 @@ read_yes_no () { VALID=1 while test $VALID -ne 0; do printf "${PROMPT}" - if [ x"$BYPASS" = "xn" ]; then + if [ x"$BYPASS" = "xn" ]; then read "RESPONSE" else printf "\n" @@ -128,7 +129,7 @@ read_folder () { ANSWER="" unset FOLDER printf "${PROMPT}" - if [ x"$BYPASS" = "xn" ]; then + if [ x"$BYPASS" = "xn" ]; then read "ANSWER" else printf "\n" @@ -170,7 +171,7 @@ wrap () { if [ "$DEBUG" = "y" ]; then eval "($1) 2>&1" | tee -a ${LOGFILE} echo ${PIPESTATUS[0]} >${GETOUTFILE} - + rc=$(cat $GETOUTFILE) if [ $rc != 0 ]; then error_exit @@ -194,7 +195,7 @@ update_path () { ${ECHO} "" ${ECHO} "Please source again your shell configuration file by typing:" ${ECHO} "source $SHELLRC" - fi + fi } @@ -230,7 +231,7 @@ patience () { ${ECHO} "\b\b Done." } -show_help () { +show_help () { echo "Usage: $SCRIPT_NAME [OPTIONS]" echo "List of possible options:" echo " -p FOLDER Default FOLDER for installation." @@ -252,7 +253,8 @@ check_dependencies () { gcc_ver_min=7 log_message "Checking if gfortran >= $gcc_ver_min is installed..." command -V gfortran || return 1 - gcc_ver=$(gcc -dumpversion) +# gcc_ver=$(gcc -dumpversion) + gcc_ver=$(echo $(gcc -dumpversion)| cut -d. -f1) echo "You have gfortran version $gcc_ver installed" test $gcc_ver -ge $gcc_ver_min || return 1 log_message "Ok\n" @@ -278,9 +280,9 @@ check_dependencies () { # we need virtualenv log_message "Checking if virtualenv is installed..." - command -V virtualenv || return 1 + #command -V virtualenv || return 1 log_message "Ok\n" - + } local_install () { @@ -322,7 +324,7 @@ local_install () { y) # Create the destination folder read_folder "Please, type in the base installation folder" "${DEFAULT_DEST}" DEST_FOLDER=${FOLDER}/MsSpec-${VERSION} - if [ -d $DEST_FOLDER ]; then + if [ -d $DEST_FOLDER ]; then export DEST_ALREADY_EXISTS="y" # to avoid cleaning the existing DEST_FOLDER in case of any problem read_yes_no "The folder $DEST_FOLDER already exists. Overwrite" "y" case ${ANSWER} in @@ -348,7 +350,7 @@ local_install () { # build and run setuptools to create a source distribution wrap "make sdist" \ "Building MsSpec python package" - + # install the package with pip wrap "pip install $pip_opt dist/msspec-*.tar.gz" \ "Installing pymsspec python package and its dependencies" @@ -377,7 +379,7 @@ local_install () { if [ x"${ANSWER}" = "xn" ]; then rm -f "${LOGFILE}" fi - + # self-explanatory success_message } diff --git a/src/msspec/phagen/fortran/phagen_scf.f b/src/msspec/phagen/fortran/phagen_scf.f index 3abcb83..3e20fc2 100644 --- a/src/msspec/phagen/fortran/phagen_scf.f +++ b/src/msspec/phagen/fortran/phagen_scf.f @@ -27,8 +27,8 @@ c .. INCOMING WAVE BOUNDARY CONDITIONS c C .................................... C -C bug corrected in subroutine -C GET_CORE_STATE +C bug corrected in subroutine +C GET_CORE_STATE C (FDP 18th May 2006) C C bug corrected in subroutine @@ -60,27 +60,27 @@ C C modified to impose lmaxt externally C (CRN : july 2009) C -C modified to include quadrupole +C modified to include quadrupole C radial matrix elements C (CRN : june 2012) C C File formats for radial integrals C modified (DS 8th january 2013) C -C modified to introduce t-matrix +C modified to introduce t-matrix C calculation in the eikonal approximation C (CRN : march 2013) C C bug corrected in routine linlogmesh: rhon ---> r_sub C (CRN : april 2013) C -C modified to calculate tmatrix, radial integrals +C modified to calculate tmatrix, radial integrals C and atomic cross sections on linearlog mesh C (CRN: september 2012 and april 2013) C C bug corrected in routine pgenll2: complex*16 dnm. C v potential converted to complex*16 in routines -C pgenll1m and pgenll2 +C pgenll1m and pgenll2 C (CRN: april 2013) C C bug corrected in the calculation of the total mfp = amfpt @@ -101,7 +101,7 @@ C statement 13824 C C bug corrected in subroutine calc_edge (xion = 0 for ground state) C (CNR: June 2017) - + implicit real*8 (a-h,o-z) c include 'msxas3.inc' @@ -140,7 +140,7 @@ c if (calctype.eq.'xpd') then call system('mkdir -p tl') C!!!! INPUT FILE TO LOAD open (idat,file='../input/input.ms',status='old') -C!!!! +C!!!! open (iphas,file='div/phases.dat',form='formatted', 1 status='unknown') open (iedl0,file='div/exdl0.dat',form='unformatted', @@ -296,7 +296,7 @@ c vinput = .false. potype='hedin' potgen='in' - cip=0.0 + cip=0.0 relc='nr' eikappr=' no' coor='angs' @@ -340,7 +340,7 @@ c c.....definition of lmax_mode: c..... lmax_mode = 0: lmaxn(na)=lmax_, independent of energy and atom number c..... lmax_mode = 1: lmaxn(na)= km*rs(na)+1, where km=(emax)^{1/2} -c..... lmax_mode = 2: lmaxn(na)= ke*rs(na)+1, where ke=(e)^{1/2}, where +c..... lmax_mode = 2: lmaxn(na)= ke*rs(na)+1, where ke=(e)^{1/2}, where c..... e is the running energy c c.. read control cards in namelist &job @@ -348,15 +348,15 @@ c read(idat,job) read(idat,*) c -c.....convert lengths in au if coor='angs'. Coordinates will be converted +c.....convert lengths in au if coor='angs'. Coordinates will be converted c in subroutine inoor - if(coor.eq.'angs'.and.lambda.ne.0) then + if(coor.eq.'angs'.and.lambda.ne.0) then lambda = lambda/real(antoau) else lambda = 20.0 ! in au corresponding to kappa = 0.05 (see subroutine cont) endif c - if(coor.eq.'angs'.and.rsh.ne.0) then + if(coor.eq.'angs'.and.rsh.ne.0) then rsh = rsh/antoau else rsh = 1.0d0 ! in au @@ -381,12 +381,12 @@ c call exit endif c - if(calctype.eq.'els') then + if(calctype.eq.'els') then lmax_mode = 2 einl = dble(einc - esct - cip) if(cip.ne.0.0.and.einl.lt.0.0d0) then write(6,*)' unable to excite chosen edge:', - & ' einc - esct - cip less than zero =', einl + & ' einc - esct - cip less than zero =', einl call exit endif endif @@ -1239,7 +1239,7 @@ c einl = dble(einc - esct - cip) if(einl.lt.0.0d0) then write(6,*)' unable to excite chosen edge:', - & ' einc - esct - cip less than zero =', einl + & ' einc - esct - cip less than zero =', einl call exit endif endif @@ -2352,7 +2352,8 @@ c dgc contains large component radial functions c common /deux/ dvn(251), dvf(251), d(251), dc(251), dgc(251,30) c passc and rho0 contain 4*pi*r^2*rho(r) c - dimension r(mp),r_hs(441),rho0_hs(441) +cKMD dimension r(mp),r_hs(441),rho0_hs(441) + dimension r(mp),r_hs(*),rho0_hs(*) C dimension dum1(mp), dum2(mp) dimension vcoul(mp), rho0(mp), enp(ms) @@ -3456,15 +3457,15 @@ C 1020 KC=KC+1 1000 CONTINUE C WRITE(7,*) 'CHECKING MUFFIN-TIN RADII' - IF(OPTRSH.EQ.'y') THEN + IF(OPTRSH.EQ.'y') THEN WRITE(6,*) ' MT radii for Hydrogen atoms set to rsh' WRITE(7,*) ' MT radii for Hydrogen atoms set to rsh =', RSH ELSE WRITE(6,*) ' MT radii for Hydrogen atoms determined by stdcrm', - & ' unless other options are specified' + & ' unless other options are specified' WRITE(7,*) ' MT radii for Hydrogen atoms determined by stdcrm', - & ' unless other options are specified' - ENDIF + & ' unless other options are specified' + ENDIF WRITE(7,*) ' M, Z(M), MN, Z(MN), AN(MN,M),', & ' RSC(M), RSC(MN), RS(M), RS(MN)' C @@ -3482,7 +3483,7 @@ C C C IF OPTRSH='y' MT RADIUS FOR H ATOMs SET TO RSH IN INPUT ! Added 16 Jul 2013 C - IF(NZM(M).EQ.1.AND.OPTRSH.EQ.'y') THEN + IF(NZM(M).EQ.1.AND.OPTRSH.EQ.'y') THEN WRITE(6,*) ' MT radius', RS(M),' for H atom', M, & ' set to', RSH RS(M) = RSH @@ -3685,7 +3686,7 @@ C common/options/rsh,ovlpfac,vc0,rs0,vinput,absorber,hole,mode, - & ionzst,potype,norman,coor,charelx,edge,potgen + & ionzst,potype,norman,coor,charelx,edge,potgen C C**** CONT_SUB DIMENSIONING VARIABLES @@ -5212,7 +5213,7 @@ c character*3 eikappr character*5 potype c - logical do_r_in + logical do_r_in c c write(6,11) jat,jd,jf,jlmax,jn,jrd,jsd,j1d c @@ -5232,7 +5233,7 @@ ctn write(30,13) ctn 13 format(2x,' e xe natom l ' ctn $ ' atmat ') c -C WARNING: COMMONS /FCNR/ AND /PARAM/ ARE AVAILABLE ONLY AFTER SUBROUTINE +C WARNING: COMMONS /FCNR/ AND /PARAM/ ARE AVAILABLE ONLY AFTER SUBROUTINE C INPUT_CONT IS CALLED c c do not change in this version! @@ -5280,7 +5281,7 @@ c 770 format( 1x,' generating t_l (for030) and', &' atomic cross section (for050)') c -c construct log-linear x mesh +c construct log-linear x mesh c call llmesh c @@ -5322,7 +5323,7 @@ c C COMMON /LLM/ ALPHA, BETA c - COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), + COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), & PAX(RDX_,F_), RAMFNR(N_), RAMFSR(N_), RAMFSOP(N_), & RAMFSOA(N_) COMPLEX PX, PX0, PPX, PAX, RAMFNR, RAMFSR, RAMFSOP, RAMFSOA @@ -5655,13 +5656,13 @@ c write(6,*) ' ' c c nels = 1 - if(calctype.eq.'els'.or.calctype.eq.'e2e') then + if(calctype.eq.'els'.or.calctype.eq.'e2e') then c nels = 3 c c calculate cluster size for effective integration of eels tme c - kappa = 1.d0/dble(lambda) ! to account for thomas-fermi screening - ! length = 2.9*0.529/(r_s)^(1/2) + kappa = 1.d0/dble(lambda) ! to account for thomas-fermi screening + ! length = 2.9*0.529/(r_s)^(1/2) ! default = 1/20 = 0.05 (au)^{-1} c do i = 1, ndat @@ -5704,7 +5705,7 @@ c.....Calculation of tl and rme for xpd, xas and rexs c c if (calctype.eq.'xpd'.or.calctype.eq.'xas'.or. - 1 calctype.eq.'rex' .or. calctype.eq.'aed'.or. + 1 calctype.eq.'rex' .or. calctype.eq.'aed'.or. 2 calctype.eq.'led') then c nks = 1 !ficticious: in this section only for writing purposes @@ -5717,8 +5718,8 @@ c c if(calctype.eq.'xpd'.or.calctype.eq.'xas'.or. 1 calctype.eq.'rex') then - write(55,830) - write(55,840) + write(55,830) + write(55,840) write(55,850) write(55,840) endif @@ -5731,19 +5732,19 @@ c c calculate energy dependent potential: c if( irho .ne. 0 ) then - if(ne.eq.1) write(6,*) ' irho =', irho, + if(ne.eq.1) write(6,*) ' irho =', irho, & ' entering vxc to calculate energy', & ' dependent exchange' call vxc ( doit ) else - if(ne.eq.1.and.nks.eq.1) then - write(6,*) ' irho =', irho, ' energy independent potential' + if(ne.eq.1.and.nks.eq.1) then + write(6,*) ' irho =', irho, ' energy independent potential' write(6,*)' constant interstitial potential vcon =', vcon endif endif - ev=e-vcon - write(6,*) ' energy dependent vcon = ', vcon,' at energy', e -C + ev=e-vcon + write(6,*) ' energy dependent vcon = ', vcon,' at energy', e +C C CONSTRUCT RELATIVISTIC POTENTIAL ON LINEAR-LOG MESH C CALL VREL @@ -5843,15 +5844,15 @@ c c c calculate atomic t-matrix elements atm(n) C -c if(ne.eq.1.and.nks.eq.1) write(6,*) - if(ne.eq.1) write(6,*) +c if(ne.eq.1.and.nks.eq.1) write(6,*) + if(ne.eq.1) write(6,*) & ' calculating atomic t-matrix elements atm(n)' c call smtx(ne,lmax_mode) c c calculate the radial integrals of transition matrix elements: -c - if(calctype.ne.'led') then +c + if(calctype.ne.'led') then call radial(doit,imvhl) endif @@ -5871,10 +5872,10 @@ c if(lxp.gt.f_) lxp=f_ - 1 call writewf(lxp) c -c energy dependent factors for dipole and quadrupole absoprtion; +c energy dependent factors for dipole and quadrupole absoprtion; c factor 1/3 for unpolarized absorption c - if(ne.eq.1) + if(ne.eq.1) & write(6,*) ' check ionization potential:', cip edfct= facts*(cip+e)*2./3.0 edfctq = 2.0/5.0*3.0/16.0*edfct*((cip+e)*fsc)**2 @@ -5886,7 +5887,7 @@ c write(6,44) 44 format(' --------------------------------------------------', 1 '---------------') - if(gamma.ne.0.0.and.ne.eq.1.and.nks.eq.1) then + if(gamma.ne.0.0.and.ne.eq.1.and.nks.eq.1) then amfph = 0.529/gamma/2 write(6,43) amfph,e 43 format(' average mean free path due to finite gamma: mfp =' @@ -5935,13 +5936,13 @@ c * ,f10.5,' angstrom at energy ', f10.5 ,/) endif 802 continue - if(gamma.ne.0.0.and.ne.eq.1) then + if(gamma.ne.0.0.and.ne.eq.1) then amfpt = 0.529/(amfpi + gamma)/2.0 write(6,46) amfpt, e endif 46 format(' total mean free path due to Im V and gamma: mfp =' * ,f10.5,' angstrom at energy ', f10.5) - if(ne.eq.1.and.amfpt.eq.0.0.and.nks.eq.1) write(6,*) + if(ne.eq.1.and.amfpt.eq.0.0.and.nks.eq.1) write(6,*) & ' infinite mean free path for gamma: mfp = 0.0 and Im V = 0.0 ' write(6,44) write(6,*) ' ' @@ -5952,7 +5953,7 @@ c write(50,*)' &&&&&&&&&&&&&&&&&&&&&&&&& ' write(50,*)' ------------------------- ' c - if (xasxpd) then + if (xasxpd) then write(50,*) ' dipole atomic cross section' else write(50,*) ' dipole rexs matrix elements' @@ -5963,8 +5964,8 @@ c do 800 i=1,2 if((l0i.eq.0).and.(i.eq.1)) goto 800 np= l0i + (-1)**i - amem = dmx(i) - amem1 = dmx1(i) + amem = dmx(i) + amem1 = dmx1(i) pamel = amem1*cmplx(atm(nstart+np))*edfct c write(50,*)'nr ', amem1*xe/pai/(l0i - 1 + i) cofct(ne,i) = amem*cmplx(atm(nstart+np))**2*edfct*xe/pai @@ -5976,10 +5977,10 @@ c write(50,*)'nr ', amem1*xe/pai/(l0i - 1 + i) rexssme = dmx1(i)/(l0i-1+i) c cofct(ne,i) = cofct(ne,i)/sigma0 c write(6,*) sigma0,sigma0r - if (calctype.eq.'xas') then + if (calctype.eq.'xas') then write(50,805) e,vcon,amfpt,sigma0,rexsrme,rexssme else - write(50,806) e,vcon,amfpt,rexsrme,rexssme + write(50,806) e,vcon,amfpt,rexsrme,rexssme endif c if(i.eq.2) write(98,*) e*13.6, sigma0 @@ -6004,8 +6005,8 @@ c lf = l0i + i if(lf.le.0) go to 900 np = l0i + i - amem = qmx(n) - amem1 = qmx1(n) + amem = qmx(n) + amem1 = qmx1(n) pamel = amem1*cmplx(atm(nstart+np))*edfctq qcofct(ne,n) = amem*cmplx(atm(nstart+np))**2*edfctq*xe/pai pamel0 = qcofct(ne,n)/cmplx(atm(nstart+np)) @@ -6020,17 +6021,17 @@ c write(6,*) sigma0,sigma0r write(50,805) e,vcon,amfpt,sigma0,rexsrme,rexssme else write(50,806) e,vcon,amfpt,rexsrme,rexssme - endif + endif 900 continue c if (xasxpd) then write(50,*)' ------------------------- ' - write(50,*) ' dipole and quadrupole cross section with ', + write(50,*) ' dipole and quadrupole cross section with ', & 'relativistic corrections of type: ', relc write(50,*)' ------------------------- ' else write(50,*)' ------------------------- ' - write(50,*) ' dipole and quadrupole rexs matrix elements', + write(50,*) ' dipole and quadrupole rexs matrix elements', & ' with relativistic corrections of type: ', relc write(50,*)' ------------------------- ' endif @@ -6047,7 +6048,7 @@ c do 910 i=1,2 if((l0i.eq.0).and.(i.eq.1)) goto 910 np= l0i + (-1)**i - amem = dmxx(i) + amem = dmxx(i) amem1 = dmxx1(i) if(relc.eq.'nr') then atmd = atmnr(nstart+np) @@ -6067,7 +6068,7 @@ c write(50,*)'nr-rc ', amem1*xe/pai/(l0i - 1 + i) rexssme = dmxx1(i)/(l0i-1+i) c cofct(ne,i) = cofct(ne,i)/sigma0 c write(6,*) sigma0,sigma0r - if (calctype.eq.'xas') then + if (calctype.eq.'xas') then write(50,805) e,vcon,amfpt,sigma0,rexsrme,rexssme else write(50,806) e,vcon,amfpt,rexsrme,rexssme @@ -6095,7 +6096,7 @@ c lf = l0i + i if(lf.le.0) go to 920 np = l0i + i - amem = qmxx(n) + amem = qmxx(n) amem1 = qmxx1(n) if(relc.eq.'nr') then atmd = atmnr(nstart+np) @@ -6103,7 +6104,7 @@ c atmd = atmsr(nstart+np) else atmd = atmsop(nstart+np) - endif + endif pamel = amem1*atmd*edfctq qcofct(ne,n) = amem*atmd**2*edfctq*xe/pai pamel0 = qcofct(ne,n)/atmd @@ -6133,7 +6134,7 @@ c do 930 i=1,2 if((l0i.eq.0).and.(i.eq.1)) goto 930 np= l0i + (-1)**i - amem = dmxxa(i) + amem = dmxxa(i) amem1 = dmxxa1(i) atmd = atmsoa(nstart+np) pamel = amem1*atmd*edfct @@ -6175,7 +6176,7 @@ c lf = l0i + i if(lf.le.0) go to 940 np = l0i + i - amem = qmxxa(n) + amem = qmxxa(n) amem1 = qmxxa1(n) atmd = atmsoa(nstart+np) pamel = amem1*atmd*edfctq @@ -6210,15 +6211,15 @@ c 2 0.0,0.0, c 3 0.0,0.0, c 4 csqrt(qmx(3)*xe/pai) C - elseif(l0i.eq.1) then + elseif(l0i.eq.1) then C c write(55,860) csqrt(dmx(1)*xe/pai/l0i), c 1 csqrt(dmx(2)*xe/pai/(l0i+1)), -c 2 0.0,0.0, -c 3 csqrt(qmx(2)*xe/pai), +c 2 0.0,0.0, +c 3 csqrt(qmx(2)*xe/pai), c 4 csqrt(qmx(3)*xe/pai) C - else + else C c write(55,860) csqrt(dmx(1)*xe/pai/l0i), c 1 csqrt(dmx(2)*xe/pai/(l0i+1)), @@ -6238,15 +6239,15 @@ C 3 0.0,0.0, 4 csqrt(qmxx(3)*xe/pai),reg_type C - elseif(l0i.eq.1) then + elseif(l0i.eq.1) then C write(55,860) csqrt(dmxx(1)*xe/pai/l0i), 1 csqrt(dmxx(2)*xe/pai/(l0i+1)), - 2 0.0,0.0, - 3 csqrt(qmxx(2)*xe/pai), + 2 0.0,0.0, + 3 csqrt(qmxx(2)*xe/pai), 4 csqrt(qmxx(3)*xe/pai),reg_type C - else + else C write(55,860) csqrt(dmxx(1)*xe/pai/l0i), 1 csqrt(dmxx(2)*xe/pai/(l0i+1)), @@ -6267,15 +6268,15 @@ C 3 0.0,0.0, 4 csqrt(qmxxa(3)*xe/pai) C - elseif(l0i.eq.1) then + elseif(l0i.eq.1) then C write(55,860) csqrt(dmxxa(1)*xe/pai/l0i), 1 csqrt(dmxxa(2)*xe/pai/(l0i+1)), - 2 0.0,0.0, - 3 csqrt(qmxxa(2)*xe/pai), + 2 0.0,0.0, + 3 csqrt(qmxxa(2)*xe/pai), 4 csqrt(qmxxa(3)*xe/pai) C - else + else C write(55,860) csqrt(dmxxa(1)*xe/pai/l0i), 1 csqrt(dmxxa(2)*xe/pai/(l0i+1)), @@ -6287,7 +6288,7 @@ C c endif c - if(calctype.ne.'xpd') then + if(calctype.ne.'xpd') then if(l0i.eq.0) then c write(55,*) '========dq irregular me: hs mesh===============' C @@ -6365,7 +6366,7 @@ c if(expmode.eq.'cis') then if(nks.eq.1) e = einc if(nks.eq.2) e = einc - cip - emin - deltae - if(nks.eq.3) e = emin + deltae + if(nks.eq.3) e = emin + deltae elseif(expmode.eq.'cfs') then if(nks.eq.1) e = esct + cip + emin + deltae if(nks.eq.2) e = esct @@ -6373,7 +6374,7 @@ c elseif(expmode.eq.'cel') then if(nks.eq.1) e = einc + deltae if(nks.eq.2) e = einc - cip - emin + deltae - if(nks.eq.3) e = emin + if(nks.eq.3) e = emin endif c ev=e-vcon @@ -6386,24 +6387,24 @@ c c calculate energy dependent potential: c if( irho .ne. 0 ) then - if(ne.eq.1) write(6,*) ' irho =', irho, + if(ne.eq.1) write(6,*) ' irho =', irho, & ' entering vxc to calculate energy', & ' dependent exchange' call vxc ( doit ) else - if(ne.eq.1.and.nks.eq.1) then + if(ne.eq.1.and.nks.eq.1) then write(6,*) ' irho =', irho, ' energy independent', - 1 ' potential' + 1 ' potential' write(6,*)' constant interstitial potential vcon =', 1 vcon endif endif ev=e-vcon - if( irho .ne. 0 ) + if( irho .ne. 0 ) & write(6,*) ' energy dependent vcon = ', vcon, 1 ' at energy', e,' Ryd' - -C + +C C CONSTRUCT RELATIVISTIC POTENTIAL ON LINEAR-LOG MESH C CALL VREL @@ -6514,7 +6515,7 @@ c c c and corresponding radial integrals of transition matrix elements: c - if(nks.eq.3) then + if(nks.eq.3) then write(55,823) ne ! energy point call radialx_eels(neff) call writeelswf @@ -6763,7 +6764,7 @@ c iout = 5 id=1 n = nas c -c kx = kmax(n) ! value used in older versions (contains the 3 points +c kx = kmax(n) ! value used in older versions (contains the 3 points C outside the muffin-tin radius that were used for interpolation) c kx = kmax(n) - 3 @@ -7121,7 +7122,7 @@ c $n_=ltot_*ua_,rd_=440,sd_=ua_-1) C c.....this subroutine calculates the radial matrix elements -c.....necessary for eels cross-section +c.....necessary for eels cross-section c.....using a linear-log mesh c common/mtxele/ nstart,nlast @@ -7152,7 +7153,7 @@ C C COMMON /PDQX/ PX(RDX_,F_),DPX(RDX_,F_),PSX(F_),DPSX(F_),RAMFX(N_) C COMPLEX PX,DPX,PSX,DPSX,RAMFX c - COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), + COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), & PAX(RDX_,F_), RAMFNR(N_), RAMFSR(N_), RAMFSOP(N_), & RAMFSOA(N_) COMPLEX PX, PX0, PPX, PAX, RAMFNR, RAMFSR, RAMFSOP, RAMFSOA @@ -7202,7 +7203,7 @@ C C C c*************************************************************************** -c note that here rpix(k) = r**3*pi(k). +c note that here rpix(k) = r**3*pi(k). c wf rpix(k) is already normalized c (see subroutine corewf) c*************************************************************************** @@ -7212,15 +7213,15 @@ c id = 1 na = nas c -c.....calculate direct and exchange Coulomb integral on absorber and different -c.....spheres +c.....calculate direct and exchange Coulomb integral on absorber and different +c.....spheres c - nt0a=n0(na) + nt0a=n0(na) ntxa=nt0a+nterms(na)-1 dxa = hx(na) nstart = nt0a nlast = ntxa -c write(6,*) 'in radialx_eels', nt0a, ntxa +c write(6,*) 'in radialx_eels', nt0a, ntxa c write(6,*) ' ' write(6,*)' writing eels (e2e) regular direct terms' @@ -7236,8 +7237,8 @@ c c do 30 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 40 n2 = nt0b, ntxb @@ -7260,12 +7261,12 @@ c l1 = lc + 1 if(l1.gt.lexp_) cycle call coulss(rid1,rid2,il(1,l1,na), - & kl(1,l1,na),kmx(na),dxa,pi,vc) + & kl(1,l1,na),kmx(na),dxa,pi,vc) write(55,10) na, l, lp, ls, lc, vc/ramfprd !, vc enddo endif c - 50 continue + 50 continue c 40 continue c @@ -7286,8 +7287,8 @@ c c do 130 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 140 n2 = nt0b, ntxb @@ -7318,10 +7319,10 @@ c if(abs(vcdr).lt.1.e-9) cycle write(55,11) na, nb, l, lp, ls, lc, lcp, vcdr enddo - enddo + enddo endif c - 150 continue + 150 continue c 140 continue c @@ -7343,8 +7344,8 @@ c c do 31 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 41 n2 = nt0b, ntxb @@ -7365,14 +7366,14 @@ c if(na.eq.nb) then do lc = lc_min, lc_max, 2 l1 = lc + 1 - if(l1.gt.lexp_) cycle + if(l1.gt.lexp_) cycle call coulss(rid3,rid4,il(1,l1,na), & kl(1,l1,na),kmx(na),dxa,pi,vcx) write(55,10) na, l, lp, ls, lc, vcx/ramfprx enddo endif c - 51 continue + 51 continue c 41 continue c @@ -7393,8 +7394,8 @@ C c do 131 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 141 n2 = nt0b, ntxb @@ -7428,7 +7429,7 @@ c enddo endif c - 151 continue + 151 continue c 141 continue c @@ -7462,8 +7463,8 @@ c c do 32 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 42 n2 = nt0b, ntxb @@ -7491,12 +7492,12 @@ c call sstrop(rid2,il(1,l1,na), & kl(1,l1,na),kmx(na),dxa,pi,trop) do k = 1, kmx(na) - rid4(k) = rid1(k)*trop(k) - rid3(k) = rid(k)*trop(k) + rid4(k) = rid1(k)*trop(k) + rid3(k) = rid(k)*trop(k) enddo call irregint1(rid3,rid4,kmx(na),dxa,vc) if(abs(vc/ramfsr3(l+1,na)).lt.1.e-10) cycle - write(55,10) na, l, lp, ls, lc, vc/ramfsr3(l+1,na) + write(55,10) na, l, lp, ls, lc, vc/ramfsr3(l+1,na) enddo else do lc=abs(l-l0i), l+l0i, 2 @@ -7516,10 +7517,10 @@ c if(abs(vcdr).lt.1.e-10) cycle write(55,11) na, nb, l, lp, ls, lc, lcp, vcdr enddo - enddo + enddo endif c - 52 continue + 52 continue c 42 continue c @@ -7547,8 +7548,8 @@ c c do 33 nat2 = 1, neff nb = nat2 - if(neq(nat2).ne.0) nb = neq(nat2) - nt0b=n0(nb) + if(neq(nat2).ne.0) nb = neq(nat2) + nt0b=n0(nb) ntxb=nt0b+nterms(nb)-1 dxb = hx(nb) do 43 n2 = nt0b, ntxb @@ -7576,12 +7577,12 @@ c call sstrop(rid2,il(1,l1,na), & kl(1,l1,na),kmx(na),dxa,pi,trop) do k = 1, kmx(na) - rid4(k) = rid1(k)*trop(k) - rid3(k) = rid(k)*trop(k) + rid4(k) = rid1(k)*trop(k) + rid3(k) = rid(k)*trop(k) enddo call irregint1(rid3,rid4,kmx(na),dxa,vc) if(abs(vc/ramfsr2(l+1,na)).lt.1.e-10) cycle - write(55,10) na, l, lp, ls, lc, vc/ramfsr2(l+1,na) + write(55,10) na, l, lp, ls, lc, vc/ramfsr2(l+1,na) enddo else do lc=abs(l-l0i), l+l0i, 2 @@ -7601,10 +7602,10 @@ c if(abs(vcdr).lt.1.e-10) cycle write(55,11) na, nb, l, lp, ls, lc, lcp, vcdr enddo - enddo + enddo endif c - 53 continue + 53 continue c 43 continue c @@ -7635,13 +7636,13 @@ c dimension rho1(kmx), rho2(kmx), il(kmx), kl(kmx) dimension rid(rdx_), a(rdx_), p(rdx_) complex rho1, rho2, vc, vc1, vc2 - complex*16 rid, a, p + complex*16 rid, a, p real*8 il, kl c id = 1 do k = 1, kmx rid(k) = il(k)*dcmplx(rho2(k)) - enddo + enddo call integrcmdp(rid,dx,kmx,a,id) do k = 1, kmx rid(k) = kl(k)*dcmplx(rho2(k)) @@ -7681,7 +7682,7 @@ c id = 1 do k = 1, kmx1 rid(k) = rho1(k)*real(ila(k)) - enddo + enddo call integrcm(rid,dx1,kmx1,a1,id) c call interp(r1(kpl1-3),a1(kpl1-3),7,rs1,vc1,dummy,.false.) vc1 = a1(kmx1) @@ -7694,7 +7695,7 @@ c c call interp(r2(kpl2-3),a2(kpl2-3),7,rs2,vc2,dummy,.false.) vc2 = a2(kmx2) c - vc = vc1*vc2*8.0*pi + vc = vc1*vc2*8.0*pi return end c @@ -7706,13 +7707,13 @@ c dimension rho2(kmx), il(kmx), kl(kmx), trop(kmx) dimension rid(rdx_), a(rdx_), p(rdx_) complex rho2 - complex*16 rid, a, p, trop + complex*16 rid, a, p, trop real*8 il, kl c id = 1 do k = 1, kmx rid(k) = il(k)*dcmplx(rho2(k)) - enddo + enddo call integrcmdp(rid,dx,kmx,a,id) do k = 1, kmx rid(k) = kl(k)*dcmplx(rho2(k)) @@ -7750,7 +7751,7 @@ c call interp(r2(kpl2-3),a2(kpl2-3),7,rs2,vc2,dummy,.false.) do k = 1, kmx1 rid(k) = ila(k)*a2(kmx2)*8.0*pi enddo -c +c return end c @@ -7762,12 +7763,12 @@ c dimension rho1(kmx), rho2(kmx), il(kmx), kl(kmx) dimension rid(rdx_), a(rdx_), p(rdx_) complex rho1, rho2, vc, vc1, vc2 - complex rid, a, p, rl, hl + complex rid, a, p, rl, hl c id = 1 do k = 1, kmx rid(k) = rl(k)*dcmplx(rho2(k)) - enddo + enddo call integrcm(rid,dx,kmx,a,id) do k = 1, kmx rid(k) = hl(k)*dcmplx(rho2(k)) @@ -7806,7 +7807,7 @@ c id = 1 do k = 1, kmx rid(k) = dcmplx(rho2(k)) - enddo + enddo call integrcm(rid,dx,kmx,a,id) do k = 1, kmx rid(k) = dcmplx(rho1(k)) @@ -8343,10 +8344,10 @@ c write(6,*) ' impinging electron wave vector kappa =', real(xe) c do na=1,nuatom - write(45,*)'atom number ', na,'(z =', z(na),')' - write(35,*)'atom number ', na,'(z =', z(na),')' + write(45,*)'atom number ', na,'(z =', z(na),')' + write(35,*)'atom number ', na,'(z =', z(na),')' c write(6,*)' atom number ', na,'(z =', z(na),')' - z0 = z(na) + z0 = z(na) call tbmat(db,rs(na),kplace(na),z0,r(1,na),v(1,na),real(xe)) enddo c @@ -8366,7 +8367,7 @@ c c dimension v(kmax),r(kmax), z(rd_) complex v, z -c +c dimension x(nt_), rx(nt_), rid(nt_), rid1(nt_) c complex cu, tb, zb, z1, zx, dzx, d2zx, rid, rid1, dbf, dbs @@ -8393,7 +8394,7 @@ c do ib = 1, nb c b = exp((ib-1)*db + b0) nb = nint(rs/db) c write(6,*) 'nb =', nb - do ib = 1, nb - 1 + do ib = 1, nb - 1 b = (ib-1)*db + db c dx = 0.005 @@ -8939,7 +8940,7 @@ C C IZ=NZ2(I_ABSORBER+I_OUTER_SPHERE) c open(unit=697,file='get1.dat',status='unknown') - if(iz.eq.0) then + if(iz.eq.0) then iz=1 ! in case an empty sphere is the first atom write(6,*) ' warning check! empty sphere is the first atom ' endif @@ -9020,7 +9021,7 @@ C RETURN END c -C +C SUBROUTINE COREWF(NAS,IZC,HOLE) C INCLUDE 'msxas3.inc' @@ -9065,12 +9066,12 @@ C XION = 0.D0 CALL GET_INTRP_CORE(IZ,ITYHOLE,ITYRADIAL,XION,CWFX,RXD,KMXN) C -C*** NOTE THAT CWFX=P_NK (UPPER COMPONENT OF DIRAC EQU.) IS NOT NORMALIZED +C*** NOTE THAT CWFX=P_NK (UPPER COMPONENT OF DIRAC EQU.) IS NOT NORMALIZED C SINCE ALSO Q_NK (LOWER COMPONENT) MUST ALSO BE CONSIDERED. -C ALSO NOTE THE RELATION TO THE SCHRODINGER RADIAL FUNCTION R*R_L = P_NK. +C ALSO NOTE THE RELATION TO THE SCHRODINGER RADIAL FUNCTION R*R_L = P_NK. C THIS RELATION HOLDS IN THE LIMIT C --> INFINITY. c -c.....Find normalization constant in ll-mesh. +c.....Find normalization constant in ll-mesh. c do i = 1, kmxn xi = sngl(cwfx(i)) @@ -9099,7 +9100,7 @@ c C C C*********************************************************************** -C +C subroutine get_intrp_core(iz,ihole,i_radial,xion,cwfx,rx,kmxn) c c @@ -9144,7 +9145,7 @@ c xion=0.d0 c c compute relativistic Hartrer-Fock-Slater charge density (on log mesh) -c +c call scfdat (title, ifr, iz, ihole, xion, amass, beta, iprint, 1 vcoul, rho0, dum1, dum2, enp, eatom) c @@ -9305,7 +9306,7 @@ c kxe = nint(alog(emax_exc - emin_exc + 1.)/de_exc + 1.) kxe = nint((xemax-xemin)/xdelta + 1.) if(kxe.gt.nep_)then c write(lunout,730) kxe - write(6,730) kxe + write(6,730) kxe 730 format(//, & ' increase the dummy dimensioning variable, nep_. ', & /,'it should be at least equal to: ', i5,/) @@ -9413,14 +9414,14 @@ c c write(lunout,741) write(6,742) emax 742 format(/,1x,' lmax assignment based on', - & ' lmax = r_mt * k_max + 2',/, + & ' lmax = r_mt * k_max + 2',/, & ' at energy emax =',f12.6) c else c write(lunout,741) - write(6,743) + write(6,743) 743 format(/,1x,' lmax assignment based on', - & ' l_max = r_mt * k_e + 2',/, + & ' l_max = r_mt * k_e + 2',/, & ' where e is the running energy') c endif @@ -9460,9 +9461,9 @@ c lmax2(i) = nint(rs(i)*sqrt(emax)) + 2 if(l_max.lt.lmax2(i)) l_max=lmax2(i) if(i.eq.ndat) then -c write(lunout,844) +c write(lunout,844) write(6,844) - endif + endif 844 format(1x,' optimal lmax chosen according to the running', & ' energy e for each atom') enddo @@ -11608,7 +11609,7 @@ c iz=29 1 (ihole.eq.ilast .and. iocc(index,ihole)-real(delion).lt.1) ) then c call wlog(' Cannot remove an electron from this level',1) write(6,*)' Cannot remove an electron from level =', ihole - write(6,*) ' stop in getorb ' + write(6,*) ' stop in getorb ' stop 'GETORB-1' endif 11 continue @@ -13096,7 +13097,7 @@ c il = istrln (string) if (il .eq. 0) then if(iprint.eq.1) print 10 - write(11,10) + write(11,10) else if(iprint.eq.1) print 10, string(1:il) write(11,10) string(1:il) @@ -13296,7 +13297,7 @@ C X0 = 10.0 + LOG(ZAT) XMAX = ALPHA*R_SUB + BETA*LOG(R_SUB) KMX(N) = NINT ( (XMAX + X0 + DPAS) / DPAS ) IF(KMX(N).GT.RDX_) THEN - WRITE(6,*) + WRITE(6,*) & 'INCREASE PARAMETER RDX_. IT SHOULD BE AT LEAST ', KMX(N) CALL EXIT ENDIF @@ -13307,8 +13308,8 @@ C CHECK IN LLMESH c write(6,'(2i5,4e15.6)') n,kmx(n),rkmx,r_sub,xmax,rho_1 c flush(6) C - CALL LINLOGMESH ( I_END, HX(N), X(1,N), RX(1,N), DO_R_IN, - & KMX(N), KPLX(N), NR, RHO_1, R_SUB, R_IN, + CALL LINLOGMESH ( I_END, HX(N), X(1,N), RX(1,N), DO_R_IN, + & KMX(N), KPLX(N), NR, RHO_1, R_SUB, R_IN, & ALPHA, BETA ) c c if(n.eq.ndat) then @@ -13329,14 +13330,14 @@ c return end c - subroutine linlogmesh ( i_end, drho, rho, r_real, do_r_in, - & kmax, kplace, nr, rho_1, r_sub, r_in, + subroutine linlogmesh ( i_end, drho, rho, r_real, do_r_in, + & kmax, kplace, nr, rho_1, r_sub, r_in, & alpha, beta ) ! ! Set up log + linear radial mesh. ! ! rho = alpha * r_real + beta * log ( r_real ) -! +! ! rho_i = rho_{i-1} + drho ! ! @@ -13351,7 +13352,7 @@ c ! rho_1 : the first point in loglinear space ! r_sub : radius of bounding sphere in loglinear space, r_sub => rho(kplace) ! r_in : -! alpha : parameter for linear part +! alpha : parameter for linear part ! beta : parameter for log part c implicit double precision (a-h,o-z) @@ -13411,7 +13412,7 @@ c check = .true. ! rn = ( rhon - beta * log ( rhon ) ) / alpha ! correction 2nd April 2013 rn = ( rhon - beta * log ( r_sub ) ) / alpha ! - do i = kplace - 1, 1, - 1 + do i = kplace - 1, 1, - 1 k = 0 ! @@ -13434,7 +13435,7 @@ c check = .true. ! if ( myrank .eq. 0 ) then - write ( unit = 6, fmt = * ) "Error occurred at radialmesh!", + write ( unit = 6, fmt = * ) "Error occurred at radialmesh!", & "rn = 0" end if @@ -13446,7 +13447,7 @@ c check = .true. end if ! - epsilon = ( alpha * rn + beta * log ( rn ) - rho ( i ) ) / + epsilon = ( alpha * rn + beta * log ( rn ) - rho ( i ) ) / & ( alpha * rn + beta ) ! ! MPI @@ -13507,7 +13508,7 @@ c check = .true. ! MPI ! - epsilon = ( alpha * rn + beta * log ( rn ) - rho ( i ) ) / + epsilon = ( alpha * rn + beta * log ( rn ) - rho ( i ) ) / & ( alpha * rn + beta ) ! ! MPI @@ -13545,18 +13546,18 @@ c check = .true. ! if ( check .and. myrank .eq. 0 ) then - write ( unit = 99, fmt = * ) '# i rho r rho ( r )', + write ( unit = 99, fmt = * ) '# i rho r rho ( r )', & ' dr' i = 1 - write ( unit = 99, fmt = "( i4, 4es20.10 )" ) i, rho ( i ), - & r_real ( i ), + write ( unit = 99, fmt = "( i4, 4es20.10 )" ) i, rho ( i ), + & r_real ( i ), & alpha * r_real ( i ) + beta * log ( r_real ( i ) ) ! do i = 2, nr write ( unit = 99, fmt = "( i4, 4es20.10 )" ) i,rho ( i ), - & r_real ( i ), - & alpha * r_real ( i ) + beta * log ( r_real ( i ) ), + & r_real ( i ), + & alpha * r_real ( i ) + beta * log ( r_real ( i ) ), & r_real ( i ) - r_real ( i - 1 ) end do @@ -13620,7 +13621,7 @@ C COMMON /FCNR/KXE,H(D_),VCONS(2), 1 R(RD_,D_),V(RD_,SD_),ICHG(10,D_),KPLACE(AT_),KMAX(AT_) COMPLEX VCONS,V -C +C COMMON /FCNRLM/X(RDX_,D_), RX(RDX_,D_), HX(D_), VX(RDX_,SD_), & VXR(RDX_,SD_), DVX(RDX_,SD_), BX(RDX_,SD_), & VXSO(RDX_,SD_), KMX(AT_), KPLX(AT_) @@ -13638,9 +13639,9 @@ c COMPLEX ZTMP(0:RD_), ZX, DZX, D2ZX REAL*4 RTMP(0:RD_) C - DATA FSC,FSCS4 /7.29735E-3,1.331283E-5/ + DATA FSC,FSCS4 /7.29735E-3,1.331283E-5/ C -C INTERPOLATE POTENTIAL ON THE LOG-LINEAR MESH +C INTERPOLATE POTENTIAL ON THE LOG-LINEAR MESH C AND ADD RELATIVISTIC CORRECTIONS, INCLUDING SPIN-ORBIT INTERACTION C C WRITE(7,*) ' I RX(I), VX(I), VXSR(I), VXSO(I), BX(I) ' @@ -13657,7 +13658,7 @@ C ENDDO C NS = N - DO IS=1,NSPINS + DO IS=1,NSPINS DO I = 1, KMAX(N) ZTMP(I) = -V(I,NS) * RTMP(I) C WRITE(6,*) N, IS, I, RTMP(I), ZTMP(I) @@ -13684,19 +13685,19 @@ C VX(I,NS) = -ZX/RX(I,N) BX(I,NS) = FSCS4/(1.0 + FSCS4*(E - VX(I,NS))) DVX(I,NS) = -(DZX/RX(I,N) - ZX/RX(I,N)**2) - VXR(I,NS) = VX(I,NS) - FSCS4*(E - VX(I,NS))**2 + - & 0.5*BX(I,NS)*( -D2ZX/RX(I,N) + + VXR(I,NS) = VX(I,NS) - FSCS4*(E - VX(I,NS))**2 + + & 0.5*BX(I,NS)*( -D2ZX/RX(I,N) + & 1.5*BX(I,NS)*(DVX(I,NS))**2 ) VXSO(I,NS) = BX(I,NS)*DVX(I,NS)/RX(I,N) -C WRITE(15,1) I, RX(I,N), VX(I,NS), VXR(I,NS), -C & VXSO(I,NS), BX(I,NS) -1 FORMAT(I5,9E15.6) +C WRITE(15,1) I, RX(I,N), VX(I,NS), VXR(I,NS), +C & VXSO(I,NS), BX(I,NS) +1 FORMAT(I5,9E15.6) ENDDO NS=NS+NDAT ENDDO C ENDDO -C +C RETURN C END @@ -13822,15 +13823,15 @@ C C COMMON /FCNR/KXE, H(D_),VCONS(2), 1 R(RD_,D_),V(RD_,SD_),ICHG(10,D_),KPLACE(AT_),KMAX(AT_) - COMPLEX VCONS,V -C + COMPLEX VCONS,V +C COMMON /FCNRLM/X(RDX_,D_), RX(RDX_,D_), HX(D_), VX(RDX_,SD_), & VXR(RDX_,SD_), DVX(RDX_,SD_), BX(RDX_,SD_), & VXSO(RDX_,SD_), KMX(AT_), KPLX(AT_) COMPLEX VX, VXR, DVX, BX, VXSO C COMPLEX VXP(RDX_), VXA(RDX_), BD(RDX_) -C +C COMPLEX PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), PAX(RDX_,F_) COMPLEX PSX(N_), DPSX(N_), STMAT, RAMFX(N_) COMPLEX PS0(N_), DPS0(N_), STMAT0, RAMF0(N_) @@ -13845,7 +13846,7 @@ c common /flag/ inmsh,inv,inrho,insym,iovrho,iosym, 1 imvhl,nedhlp c - complex pss(6),dpss(6), + complex pss(6),dpss(6), & ramfnr(n_), ramfsr(n_), ramfsop(n_), ramfsoa(n_) c character*8 name0 ,nsymbl !added 29/3/2013 @@ -13899,7 +13900,7 @@ C atmnr(j)=(0.00,0.00) atmsr(j)=(0.00,0.00) atmsop(j)=(0.00,0.00) - 5 atmsoa(j)=(0.00,0.00) + 5 atmsoa(j)=(0.00,0.00) c c write(70,*) ' non relativistic stmat and phase shifts ' c write(80,*) ' scalar relativistic stmat and phase shifts ' @@ -13908,9 +13909,9 @@ c c calculate t-matrix elements: c stmat: inverse t-m elements (atomic spheres) c ramf: for normalization of ps(k) functions -c +c c write(19,18) e, xe - write(81,*) ' e, vcon, xe, relc =', e, real(vcon), + write(81,*) ' e, vcon, xe, relc =', e, real(vcon), & real(xe), relc c write(84,*) ' e, vcon, xe =', e, vcon, xe c 18 FORMAT(' E =', F10.5,5X,' XE =',2F10.5,' GAMMA =',F10.5) @@ -13921,7 +13922,7 @@ c write(80,77) na write(90,77) na ns=ns+1 - 25 nt0a=n0(na) + 25 nt0a=n0(na) ntxa=nt0a+nterms(na)-1 if (na.eq.nas) then nstart=nt0a @@ -13945,14 +13946,14 @@ C if(lmax_mode.eq.2.and.l.gt.lmxne(na,ne)) goto 45 np=npabs C -c if(relc.eq.'nr') then +c if(relc.eq.'nr') then c rx1 = dble(rx(1,na)) - rx2 = dble(rx(2,na)) + rx2 = dble(rx(2,na)) y0(l) = dcmplx(rx1**(l+1),0.d0) y1(l) = dcmplx(rx2**(l+1),0.d0) c - call pgenll1m(l, e, hx(na), rx(1,na), vx(1,ns), bd, + call pgenll1m(l, e, hx(na), rx(1,na), vx(1,ns), bd, & kmx(na), kplx(na), rs(na), px(1,np), psx(nn), & dpsx(nn), ramf00, stmat, y0(l),y1(l)) c @@ -13974,13 +13975,13 @@ c 1000 format(3x,f9.4,1x,f9.4,5x,f12.9,5x,f12.9,5x,f12.9,5x,f12.9) c c elseif(relc.eq.'sr') then -c +c rx1 = dble(rx(1,na)) rx2 = dble(rx(2,na)) expr = 0.5d0 + sqrt( dfloat(l*(l+1)) +1 - dble(fsc*z(na))**2 ) y0(l) = dcmplx(rx1**expr,0.d0) y1(l) = dcmplx(rx2**expr,0.d0) - call pgenll1m(l, e, hx(na), rx(1,na), vxr(1,ns), bx(1,ns), + call pgenll1m(l, e, hx(na), rx(1,na), vxr(1,ns), bx(1,ns), & kmx(na), kplx(na), rs(na), px0(1,np), ps0(nn), & dps0(nn), ramf00, stmat0, y0(l),y1(l)) c @@ -14014,7 +14015,7 @@ c if(l.eq.0) ilm = 1 do il = 1, ilm c - if(il.eq.1) then + if(il.eq.1) then do i = 1, kmx(na) vxp(i) = vxr(i,ns) + float(l)*vxso(i,ns) enddo @@ -14023,12 +14024,12 @@ c expr = 0.5d0 + sqrt( dfloat(l+1)**2 -dble(fsc*z(na))**2 ) y0(l) = dcmplx(rx1**expr,0.d0) y1(l) = dcmplx(rx2**expr,0.d0) - call pgenll1m(l, e, hx(na), rx(1,na), vxp, bx(1,ns), - & kmx(na), kplx(na), rs(na), ppx(1,np), - & ps1(nn), dps1(nn), ramf01, stmat1, + call pgenll1m(l, e, hx(na), rx(1,na), vxp, bx(1,ns), + & kmx(na), kplx(na), rs(na), ppx(1,np), + & ps1(nn), dps1(nn), ramf01, stmat1, & y0(l),y1(l)) if(na.eq.nas) - & write(81,1) 'rp', na, l, real(stmat1), 1.0/stmat1, + & write(81,1) 'rp', na, l, real(stmat1), 1.0/stmat1, & real(ramf01), e else do i = 1, kmx(na) @@ -14040,9 +14041,9 @@ c if(l.eq.0) expr = 0.5d0 +sqrt( 1.0d0 -dble(fsc*z(na))**2) y0(l) = dcmplx(rx1**expr,0.d0) y1(l) = dcmplx(rx2**expr,0.d0) - call pgenll1m(l, e, hx(na), rx(1,na), vxa, bx(1,ns), - & kmx(na), kplx(na), rs(na), pax(1,np), - & ps2(nn), dps2(nn), ramf02, stmat2, + call pgenll1m(l, e, hx(na), rx(1,na), vxa, bx(1,ns), + & kmx(na), kplx(na), rs(na), pax(1,np), + & ps2(nn), dps2(nn), ramf02, stmat2, & y0(l),y1(l)) c endif @@ -14071,7 +14072,7 @@ c c endif 1 format(a3,2i5,10e13.5) 30 format(5i3,8e13.5) -c +c c 45 continue 60 continue @@ -14082,8 +14083,8 @@ c c c calculate singular solution inside muffin tin sphere for the absorbing c atom, matching to shf in interstitial region -c - if(calctype.eq.'els'.and.nks.eq.3) +c + if(calctype.eq.'els'.and.nks.eq.3) & write(6,*)' store irregular solution' 90 nl=0 lmsing=5 @@ -14093,7 +14094,7 @@ c c if(nks.eq.3) write(6,*)' nst =',nst,' nlst =',nlst l=-1 ml=lmaxn(nas)+1 - if (ml.lt.3) ml = 3 + if (ml.lt.3) ml = 3 kpp = kmx(nas) -2 arg=xe*rx(kpp,nas) call cshf2(arg,xe,ml,sbfx,dsbfx) @@ -14111,10 +14112,10 @@ c pkmx1 = cmplx(shfx(l+1))*arg1/pi c call pgenll2( l, e, hx(nas), rx(1,nas), vx(1,nas), bd, - & kpp, px(1,np), pkmx, pkmx1 ) + & kpp, px(1,np), pkmx, pkmx1 ) - call pgenll2( l, e, hx(nas), rx(1,nas), vxr(1,nas), - & bx(1,nas), kpp, px0(1,np), pkmx, pkmx1 ) + call pgenll2( l, e, hx(nas), rx(1,nas), vxr(1,nas), + & bx(1,nas), kpp, px0(1,np), pkmx, pkmx1 ) ilm = 2 if(l.eq.0) ilm = 1 @@ -14125,10 +14126,10 @@ c enddo c do il = 1, ilm - if(il.eq.1) - & call pgenll2( l, e, hx(nas), rx(1,nas), vxp, + if(il.eq.1) + & call pgenll2( l, e, hx(nas), rx(1,nas), vxp, & bx(1,nas), kpp, ppx(1,np), pkmx, pkmx1 ) - if(il.eq.2) + if(il.eq.2) & call pgenll2( l, e, hx(nas), rx(1,nas), vxa, & bx(1,nas), kpp, pax(1,np), pkmx, pkmx1 ) enddo @@ -14156,7 +14157,7 @@ c c c - subroutine pgenll1m(l, en, h, rx, v, b, kmax, kplx, rs, + subroutine pgenll1m(l, en, h, rx, v, b, kmax, kplx, rs, & p, ps, dps, ramf, stmat, y0, y1 ) c c @@ -14178,11 +14179,11 @@ c c dimension rx(kmax) c - double precision dfl, a, hd, hsq12, rxi, den, arb2, + double precision dfl, a, hd, hsq12, rxi, den, arb2, & alphad, betad, rlv, amv complex*16 dvi c - complex*16 um(0:kmax), vm(0:kmax), + complex*16 um(0:kmax), vm(0:kmax), & am(0:kmax), bm(0:kmax) c c @@ -14203,7 +14204,7 @@ c betad = dble(beta) den = dble(en) dfl = dble(float(l)) - a = (dfl + 1)*dfl + a = (dfl + 1)*dfl hd = dble(h) hsq12 = hd*hd/12.d0 c @@ -14212,8 +14213,8 @@ c arb2 = (alphad*rxi + betad)**2 dvi = dcmplx(v(i)) am(i) = 1.d0 + 1.d0/arb2 * ( rxi**2 * (den-dvi) - a - - & betad*(alphad*rxi + betad/4.d0)/arb2 )*hsq12 - bm(i) = 2.d0*(6.d0 - 5.d0*am(i)) + & betad*(alphad*rxi + betad/4.d0)/arb2 )*hsq12 + bm(i) = 2.d0*(6.d0 - 5.d0*am(i)) enddo do i = 2, kmax-1 @@ -14227,17 +14228,17 @@ c pd(1) = y0 * sqrt( alphad + betad/dble(rx(1)) ) pd(2) = y1 * sqrt( alphad + betad/dble(rx(2)) ) do i = 2, kmax - 1 - pd(i+1) = (pd(i) - um(i)*pd(1))/vm(i) + pd(i+1) = (pd(i) - um(i)*pd(1))/vm(i) enddo c do i = 1, kmax - pd(i) = pd(i)*sqrt(dble(rx(i))/(alphad*dble(rx(i))+betad) ) * + pd(i) = pd(i)*sqrt(dble(rx(i))/(alphad*dble(rx(i))+betad) ) * & dble(fsc)/2.0D0 /sqrt(dcmplx(b(i)))/ dble(rx(i)) p(i) = cmplx(pd(i)) enddo c kplx3 = kplx - 3 - call interp(rx(kplx3),p(kplx3),7,rs,ps,dps,.true.) + call interp(rx(kplx3),p(kplx3),7,rs,ps,dps,.true.) c x=dps/ps ramff=cmplx(sbf(l+1))*x-cmplx(dsbf(l+1)) @@ -14260,15 +14261,15 @@ c complex v(kmax), p(kmax), b(kmax), pkmx, pkmx1 dimension rx(kmax) c - double precision dfl, a, hd, hsq12, rxi, den, arb2, + double precision dfl, a, hd, hsq12, rxi, den, arb2, & alphad, betad c complex*16 um(0:kmax), vm(0:kmax), am(0:kmax), bm(0:kmax) - complex*16 dvi, dnm + complex*16 dvi, dnm c data pi/3.14159265/, fsc/7.29735E-3/ c -c calculate coefficients um(m) and vm(m). +c calculate coefficients um(m) and vm(m). c vm(kmax) = (0.d0,0.d0) @@ -14278,7 +14279,7 @@ c betad = dble(beta) den = dble(en) dfl = dble(float(l)) - a = (dfl + 1)*dfl + a = (dfl + 1)*dfl hd = dble(h) hsq12 = hd*hd/12.d0 c @@ -14287,15 +14288,15 @@ c arb2 = (alphad*rxi + betad)**2 dvi = dcmplx(v(i)) am(i) = 1.d0 + 1.d0/arb2 * ( rxi**2 * (den-dvi) - a - - & betad*(alphad*rxi + betad/4.d0)/arb2 )*hsq12 - bm(i) = 2.d0*(6.d0 - 5.d0*am(i)) + & betad*(alphad*rxi + betad/4.d0)/arb2 )*hsq12 + bm(i) = 2.d0*(6.d0 - 5.d0*am(i)) enddo do i = kmax-1, 2, -1 dnm = ( bm(i) - am(i+1)*vm(i+1) ) vm(i) = am(i-1) / dnm um(i) = am(i+1) * um(i+1) / dnm -c write(6,*) vm(i), um(i) +c write(6,*) vm(i), um(i) enddo @@ -14337,7 +14338,7 @@ c C amass=0.0d0 beta=0.0d0 -c +c call scfdat (title, ifr, iz, ihole, xion, amass, beta, iprint, 1 vcoul, rho0, dum1, dum2, enp, eatom) c @@ -14378,16 +14379,16 @@ c c write(6,*) ' ---' do i = 1, 2 -c +c ityhole = ihole1 c if(i.eq.2) ityhole = ihole2 ----- corrected 23th June 2017 - if(i.eq.2) then + if(i.eq.2) then ityhole = ihole2 xion = 1.0d0 endif c if(i.eq.1) write(6,*) ' total energy for atom in ground state ' - if(i.eq.2) write(6,*) ' total energy for atom with a hole in ', + if(i.eq.2) write(6,*) ' total energy for atom with a hole in ', & edge, ' edge' c @@ -14400,12 +14401,12 @@ c write(6,*) ' calculated ionization energy for edge ', edge, & ' = ', cip*13.6, ' eV' c -c.....Find out energy distance between edges and construct two edge +c.....Find out energy distance between edges and construct two edge c dipole cross section c xion=1.0d0 c - if(edge.eq.'k '.or.edge.eq.'l1'.or.edge.eq.'m1'.or.edge.eq.'n1') + if(edge.eq.'k '.or.edge.eq.'l1'.or.edge.eq.'m1'.or.edge.eq.'n1') & go to 15 if(edge.eq.'l2'.or.edge.eq.'l3') then ihole1 = 3 @@ -14428,7 +14429,7 @@ c endif c do i = 1, 2 - + ityhole = ihole1 if(i.eq.2) ityhole = ihole2 c @@ -14438,7 +14439,7 @@ c c detot = (etot(1) - etot(2))*2.0d0 detot = sign(detot,1.0d0) - if(edge.eq.'l2'.or.edge.eq.'l3') then + if(edge.eq.'l2'.or.edge.eq.'l3') then write(6,*) ' energy distance between edges l2 and l3 = ', & real( etot(1) - etot(2) )* 27.2, 'eV' elseif(edge.eq.'m2'.or.edge.eq.'m3') then @@ -14464,7 +14465,7 @@ C C c.....this subroutine calculates the radial matrix elements d(i) c.....(i=1,2) for lfin=l0i-1 (i=1) and lfin=l0i+1 (i=2) both for -c.....the regular (dmxx) and irregular solution (dmxx1) using a +c.....the regular (dmxx) and irregular solution (dmxx1) using a c.....linear-log mesh c common/mtxele/ nstart,nlast @@ -14495,7 +14496,7 @@ C C COMMON /PDQX/ PX(RDX_,F_),DPX(RDX_,F_),PSX(F_),DPSX(F_),RAMFX(N_) C COMPLEX PX,DPX,PSX,DPSX,RAMFX c - COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), + COMMON /PDQX/PX(RDX_,F_), PX0(RDX_,F_), PPX(RDX_,F_), & PAX(RDX_,F_), RAMFNR(N_), RAMFSR(N_), RAMFSOP(N_), & RAMFSOA(N_) COMPLEX PX, PX0, PPX, PAX, RAMFNR, RAMFSR, RAMFSOP, RAMFSOA @@ -14537,7 +14538,7 @@ C C C c*************************************************************************** -c note that here rpix(k) = r**3*pi(k). +c note that here rpix(k) = r**3*pi(k). c wf rpix(k) is already normalized c (see subroutine corewf) c*************************************************************************** @@ -14607,14 +14608,14 @@ c c dmx(i) = (cri(kx)/ramf(nstart+np))**2*(l0i-1+i) DO 117 K=1,KX 117 RID(K)=RPIX(K)*PX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 118 K=1,KX 118 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 119 K=1,KX 119 RID(K)=RPIX(K)*PX(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) + CALL DEFINT1(RID,DH,KX,SMX1,ID) DMXX1(I) = (SMX0 + SMX1)*(L0I-1+I)/RAMFNR(NSTART+NP) c dmx1(i) = (dx+dx1)*(l0i-1+i)/ramf(nstart+np) c @@ -14626,14 +14627,14 @@ c DMXX(I) = (CRI(KX)/RAMFSR(NSTART+NP))**2*(L0I-1+I) DO 120 K=1,KX 120 RID(K)=RPIX(K)*PX0(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 121 K=1,KX 121 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 122 K=1,KX 122 RID(K)=RPIX(K)*PX0(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) + CALL DEFINT1(RID,DH,KX,SMX1,ID) DMXX1(I) = (SMX0 + SMX1)*(L0I-1+I)/RAMFSR(NSTART+NP) c else if(relc.eq.'so') then @@ -14644,14 +14645,14 @@ c DMXX(I) = (CRI(KX)/RAMFSOP(NSTART+NP))**2*(L0I-1+I) DO 123 K=1,KX 123 RID(K)=RPIX(K)*PPX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 124 K=1,KX 124 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 125 K=1,KX 125 RID(K)=RPIX(K)*PPX(K,NP)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) + CALL DEFINT1(RID,DH,KX,SMX1,ID) DMXX1(I) = (SMX0 + SMX1)*(L0I-1+I)/RAMFSOP(NSTART+NP) C DO K=1,KX @@ -14661,25 +14662,25 @@ C DMXXA(I) = (CRI(KX)/RAMFSOA(NSTART+NP))**2*(L0I-1+I) DO 126 K=1,KX 126 RID(K)=RPIX(K)*PAX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 127 K=1,KX 127 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,DX,ID) + CALL DEFINT1(RID,DH,KX,DX,ID) DO 128 K=1,KX 128 RID(K)=RPIX(K)*PAX(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,DX1,ID) + CALL DEFINT1(RID,DH,KX,DX1,ID) DMXXA1(I) = (DX + DX1)*(L0I-1+I)/RAMFSOA(NSTART+NP) c endif - + 100 continue C c write(6,*) ' radialx matrix elements from shell li = ', l0i c write(6,*) (real(dmxx(l)),aimag(dmxx(l)),l=1,2) c write(6,*) (real(dmxx1(l)),aimag(dmxx1(l)),l=1,2) C -C.....CALCULATE RADIAL QUADRUPOLAR TRANSITION MATRIX ELEMENT +C.....CALCULATE RADIAL QUADRUPOLAR TRANSITION MATRIX ELEMENT C DO K = 1, KX RPIX(K) = RPIX(K) * RX(K,NQ) @@ -14699,19 +14700,19 @@ c DO 216 K=1,KX 216 RID(K)=RPIX(K)*PX(K,NP+1)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) CALL INTEGRCM(RID,DH,KX,CRI,ID) - QMXX(M) = (CRI(KX)/RAMFNR(NSTART+NP))**2 + QMXX(M) = (CRI(KX)/RAMFNR(NSTART+NP))**2 c dmx(i) = (cri(kx)/ramf(nstart+np))**2*(l0i-1+i) DO 217 K=1,KX 217 RID(K)=RPIX(K)*PX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 218 K=1,KX 218 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 219 K=1,KX 219 RID(K)=RPIX(K)*PX(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) - QMXX1(M) = (SMX0 + SMX1)/RAMFNR(NSTART+NP) + CALL DEFINT1(RID,DH,KX,SMX1,ID) + QMXX1(M) = (SMX0 + SMX1)/RAMFNR(NSTART+NP) c dmx1(i) = (dx+dx1)*(l0i-1+i)/ramf(nstart+np) c else if(relc.eq.'sr') then @@ -14719,53 +14720,53 @@ c RID(K)=RPIX(K)*PX0(K,NP+1)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) ENDDO CALL INTEGRCM(RID,DH,KX,CRI,ID) - QMXX(M) = (CRI(KX)/RAMFSR(NSTART+NP))**2 + QMXX(M) = (CRI(KX)/RAMFSR(NSTART+NP))**2 DO 220 K=1,KX 220 RID(K)=RPIX(K)*PX0(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 221 K=1,KX 221 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 222 K=1,KX 222 RID(K)=RPIX(K)*PX0(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) - QMXX1(M) = (SMX0 + SMX1)/RAMFSR(NSTART+NP) + CALL DEFINT1(RID,DH,KX,SMX1,ID) + QMXX1(M) = (SMX0 + SMX1)/RAMFSR(NSTART+NP) c else if(relc.eq.'so') then DO K=1,KX RID(K)=RPIX(K)*PPX(K,NP+1)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) ENDDO CALL INTEGRCM(RID,DH,KX,CRI,ID) - QMXX(M) = (CRI(KX)/RAMFSOP(NSTART+NP))**2 + QMXX(M) = (CRI(KX)/RAMFSOP(NSTART+NP))**2 DO 223 K=1,KX 223 RID(K)=RPIX(K)*PPX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 224 K=1,KX 224 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,SMX0,ID) + CALL DEFINT1(RID,DH,KX,SMX0,ID) DO 225 K=1,KX 225 RID(K)=RPIX(K)*PPX(K,NP)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,SMX1,ID) - QMXX1(M) = (SMX0 + SMX1)/RAMFSOP(NSTART+NP) + CALL DEFINT1(RID,DH,KX,SMX1,ID) + QMXX1(M) = (SMX0 + SMX1)/RAMFSOP(NSTART+NP) C DO K=1,KX RID(K)=RPIX(K)*PAX(K,NP+1)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) ENDDO CALL INTEGRCM(RID,DH,KX,CRI,ID) - QMXXA(M) = (CRI(KX)/RAMFSOA(NSTART+NP))**2 + QMXXA(M) = (CRI(KX)/RAMFSOA(NSTART+NP))**2 DO 226 K=1,KX 226 RID(K)=RPIX(K)*PAX(K,NP+1+NPSS)*RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL INTEGRCM(RID,DH,KX,CRI1,ID) + CALL INTEGRCM(RID,DH,KX,CRI1,ID) DO 227 K=1,KX 227 RID(K)=RID(K)*CRI(K) - CALL DEFINT1(RID,DH,KX,DX,ID) + CALL DEFINT1(RID,DH,KX,DX,ID) DO 228 K=1,KX 228 RID(K)=RPIX(K)*PAX(K,NP+1)*(CRI1(KX) - CRI1(K))* & RX(K,NQ)/(ALPHA*RX(K,NQ) + BETA) - CALL DEFINT1(RID,DH,KX,DX1,ID) - QMXXA1(M) = (DX + DX1)/RAMFSOA(NSTART+NP) + CALL DEFINT1(RID,DH,KX,DX1,ID) + QMXXA1(M) = (DX + DX1)/RAMFSOA(NSTART+NP) c endif C @@ -14826,7 +14827,7 @@ C RAT=SINH(X)/(X*SBF(1)) DSBF1=SBF2*RAT GO TO 16 -11 CONTINUE +11 CONTINUE DO 12 J=1,KMAX SBFK=XF1*SBF1/X-SBF2 SBF2=SBF1 @@ -14835,8 +14836,8 @@ C IF (J.LT.JMIN) GO TO 12 SBF(K)=SBFK K=K-1 -12 CONTINUE - 15 RAT=SIN(X)/(X*SBF(1)) +12 CONTINUE + 15 RAT=SIN(X)/(X*SBF(1)) DSBF1=-SBF2*RAT 16 DO 17 K=1,MAX 17 SBF(K)=RAT*SBF(K) @@ -14913,5 +14914,5 @@ C STOP 2013 C END -C +C