From deb350a520d0a836a579aba6c2a59a02049893b1 Mon Sep 17 00:00:00 2001 From: kmdunseath <73846968+kmdunseath@users.noreply.github.com> Date: Thu, 29 Jun 2023 11:11:49 +0200 Subject: [PATCH] For spectrocopy EIG, added the possibility of using Arnoldi iteration techniques implemented in ARPACK for calculating a subset of eigenvalues rather than all of them or just the largest one. This is invoked by setting algorithm='arnoldi' in the constructor MSSPEC. calc = MSSPEC(spectroscopy='EIG', algorithm='arnoldi') This version finds NEV eigenvalues with the largest magnitudes. The value of NEV is set in the code to be 2% of the total number. The algorithm also sets the number of basis vectors, NCV, to be 2*NEV. Options for allowing the user to set these values as well as other convergence criteria should be added in later versions. Files src/msspec/calculator.py and src/msspec/parameters.py have been modified to accept algorithm='arnoldi'. A new directory, src/msspec/spec/fortran/eig/ar, has been created and contains various Fortran files implementing the Arnoldi iteration. Some files are written in Fortran90 free-format and use some features of the Fortran 2003 standard. Compilation requires the ARPACK library to be installed and accessible using -larpack Compilation rules and options have been introduced/modified where necessary in the files src/options.mk, src/msspec/spec/fortran/Makefile and src/msspec/spec/fortran/eig_ar.mk This version has been tested for the CaO example. --- src/msspec/calculator.py | 7 +- src/msspec/parameters.py | 4 +- src/msspec/spec/fortran/Makefile | 8 +- .../spec/fortran/eig/ar/arnoldi_mod.f90 | 249 +++ src/msspec/spec/fortran/eig/ar/do_main.f | 1558 +++++++++++++++++ src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 | 291 +++ src/msspec/spec/fortran/eig/ar/eigdif_ar.f | 103 ++ src/msspec/spec/fortran/eig/ar/main.f | 22 + src/msspec/spec/fortran/eig_ar.mk | 14 + src/options.mk | 10 +- 10 files changed, 2259 insertions(+), 7 deletions(-) create mode 100644 src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90 create mode 100644 src/msspec/spec/fortran/eig/ar/do_main.f create mode 100644 src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 create mode 100644 src/msspec/spec/fortran/eig/ar/eigdif_ar.f create mode 100644 src/msspec/spec/fortran/eig/ar/main.f create mode 100644 src/msspec/spec/fortran/eig_ar.mk diff --git a/src/msspec/calculator.py b/src/msspec/calculator.py index 17f36ba..37497ef 100644 --- a/src/msspec/calculator.py +++ b/src/msspec/calculator.py @@ -95,6 +95,7 @@ from msspec.parameters import TMatrixParameters from msspec.phagen.fortran.libphagen import main as do_phagen from msspec.spec.fortran import _eig_mi from msspec.spec.fortran import _eig_pw +from msspec.spec.fortran import _eig_ar from msspec.spec.fortran import _phd_mi_noso_nosp_nosym from msspec.spec.fortran import _phd_se_noso_nosp_nosym from msspec.spec.fortran import _phd_ce_noso_nosp_nosym @@ -418,6 +419,8 @@ class _MSCALCULATOR(Calculator): do_spec = _eig_mi.run elif self.global_parameters.algorithm == 'power': do_spec = _eig_pw.run + elif self.global_parameters.algorithm == 'arnoldi': + do_spec = _eig_ar.run else: LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not " "an allowed combination.".format(self.global_parameters.spectroscopy, @@ -986,8 +989,8 @@ class _EIG(_MSCALCULATOR): _MSCALCULATOR.__init__(self, spectroscopy='EIG', algorithm=algorithm, polarization=polarization, dichroism=dichroism, spinpol=spinpol, folder=folder, txt=txt) - if algorithm not in ('inversion', 'power'): - LOGGER.error("Only the 'inversion' or the 'power' algorithms " + if algorithm not in ('inversion', 'power', 'arnoldi'): + LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms " "are supported in EIG spectroscopy mode") exit(1) self.iodata = iodata.Data('EIG Simulation') diff --git a/src/msspec/parameters.py b/src/msspec/parameters.py index f4fa39a..7af83cf 100644 --- a/src/msspec/parameters.py +++ b/src/msspec/parameters.py @@ -770,7 +770,8 @@ class GlobalParameters(BaseParameters): Parameter('algorithm', types=str, allowed_values=('expansion', 'inversion', 'correlation', - 'power'), + 'power', + 'arnoldi'), default='expansion', doc=""" You can choose the algorithm used for the computation of the scattering path operator. @@ -778,6 +779,7 @@ class GlobalParameters(BaseParameters): - '**expansion**', for the Rehr-Albers series expansion - '**correlation**', for the correlation-expansion algorithm - '**power**', for the power method approximation scheme (only for spectroscopy='EIG') + - '**arnoldi**', for computing multiple eigenvalues using Arnoldi iteration (only for spectroscopy='EIG') The series expansion algorithm is well suited for high energy since the number of terms required decreases as the energy increases. The exact solution is obtained by the matrix inversion diff --git a/src/msspec/spec/fortran/Makefile b/src/msspec/spec/fortran/Makefile index 01228ea..da1e90e 100644 --- a/src/msspec/spec/fortran/Makefile +++ b/src/msspec/spec/fortran/Makefile @@ -1,6 +1,6 @@ -.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw comp_curve clean +.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve clean -all: phd_se phd_mi phd_ce eig_mi eig_pw comp_curve +all: phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve phd_se: @+$(MAKE) -f phd_se_noso_nosp_nosym.mk all @@ -17,6 +17,9 @@ eig_mi: eig_pw: @+$(MAKE) -f eig_pw.mk all +eig_ar: + @+$(MAKE) -f eig_ar.mk all + comp_curve: @+$(MAKE) -f comp_curve.mk all @@ -26,4 +29,5 @@ clean:: @+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@ @+$(MAKE) -f eig_mi.mk $@ @+$(MAKE) -f eig_pw.mk $@ + @+$(MAKE) -f eig_ar.mk $@ @+$(MAKE) -f comp_curve.mk $@ diff --git a/src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90 b/src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90 new file mode 100644 index 0000000..8c2e35e --- /dev/null +++ b/src/msspec/spec/fortran/eig/ar/arnoldi_mod.f90 @@ -0,0 +1,249 @@ +!==============================================================================! + module arnoldi_mod +!==============================================================================! + implicit none + private +! + integer, parameter :: sp = kind(1.0) + integer, parameter :: dp = kind(1.0d0) + integer, parameter :: zp = kind((0.0d0,0.0d0)) +! + complex(zp), parameter :: one = complex(1.0d0, 0.0d0) + complex(zp), parameter :: zero = complex(0.0d0, 0.0d0) +! +! Public procedures +! + public :: arnoldi_iteration +! +! Private data +! +! ARPACK's debug common block +! + integer :: logfil = 6, ndigit = -3, mcaupd = 1 + integer :: mgetv0 = 0, msaupd = 0, msaup2 = 0, msaitr = 0, mseigt = 0, & + msapps = 0, msgets = 0, mseupd = 0, mnaupd = 0, mnaup2 = 0, & + mnaitr = 0, mneigh = 0, mnapps = 0, mngets = 0, mneupd = 0, & + mcaup2 = 0, mcaitr = 0, mceigh = 0, mcapps = 0, mcgets = 0, & + mceupd = 0 +! + common /debug/ logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, & + msapps, msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, & + mnapps, mngets, mneupd, mcaupd, mcaup2, mcaitr, mceigh, & + mcapps, mcgets, mceupd +!==============================================================================! + contains +!==============================================================================! +! Public procedures +!==============================================================================! + subroutine arnoldi_iteration (ndim, a, nev, d) +! +! Use ARPACK routines to find a few eigenvalues (lambda) and corresponding +! eigenvectors (x) for the standard eigenvalue problem: +! +! A*x = lambda*x +! +! where A is a general NDIM by NDIM complex matrix +! +! This subroutine is based on an ARPACK test program adapted by Logan Boudet +! as part of his M1 project in Rennes in 2022. +! + integer, intent(in) :: ndim + complex(zp), intent(in) :: a(:,:) + integer, intent(inout) :: nev ! number of eigenvalues required + complex(zp), intent(out) :: d(:) ! vector of required eigenvalues +! +! Local data +! + character(1), parameter :: bmat = "I" ! standarg eigenvalue problem + character(2), parameter :: which = "LM" ! find NEV eigenvalues of +! ! largest magnitude +! + integer :: ido, ierr, info, j, lworkl, nconv, ncv + integer :: iparam(11) + integer :: ipntr(14) + real(dp) :: tol + complex(zp) :: sigma + logical :: rvec +! + real(dp), allocatable :: rd(:,:) + real(dp), allocatable :: rwork(:) + complex(zp), allocatable :: ax(:) + complex(zp), allocatable :: resid(:) + complex(zp), allocatable :: v(:,:) + complex(zp), allocatable :: workd(:) + complex(zp), allocatable :: workev(:) + complex(zp), allocatable :: workl(:) + logical, allocatable :: select(:) +! +! External BLAS/LAPACK functions used +! + real(dp), external :: dznrm2, dlapy2 +! +! + write(6,*) "----------------- BEGIN OF ARNOLDI ITERATION -----------------" +! +! NCV is is the largest number of basis vectors that will be used in the +! Implicitly Restarted Arnoldi Process. Work per major iteration is +! proportional to NDIM*NCV*NCV. +! +! Note: we must have NCV >= NEV + 2, and preferably NCV >= 2*NEV +! +! ncv = max(ceiling(1.125*nev + 15), nev+3) + ncv = 2*nev +! + iparam(11) = 0 + ipntr(14) = 0 +! +! stopping criteria; machine precision is used if tol <= 0 +! + tol = 0.0_dp +! +! Algorithm mode +! + iparam(1) = 1 ! exact shift strategy + iparam(3) = 300 ! maximum number of iterations + iparam(7) = 1 ! use mode1 of ZNAUPD +! + lworkl = ncv*(3*ncv + 5) + allocate(rwork(ncv)) + allocate(resid(ndim), v(ndim,ncv)) + allocate(workd(3*ndim), workev(2*ncv), workl(lworkl)) + allocate(select(ncv)) +! +! IDO is the reverse communication parameter used to determine action to be taken +! on return from ZNAUPD. Its initial value must be 0. +! + ido = 0 +! +! On entry, INFO == 0 instructs ZNAUPD to use a random starting vector. +! To specify a particular starting vector, set INFO to a non-zero value. +! The required startng vector should then be supplied in RESID. +! + info = 0 +! +! Main reverse communication loop +! + do + call znaupd(ido, bmat, ndim, which, nev, tol, resid, ncv, v, ndim, & + iparam, ipntr, workd, workl, lworkl, rwork, info) + + if (abs(ido) /= 1) exit +! +! Matrix-vector multiplication y = A*x +! Initial vector x is stored starting at workd(ipntr(1)) +! The result y should be stored in workd(ipntr(2)) +! + call matvec(a, ndim, ndim, workd(ipntr(1)), workd(ipntr(2))) +! + end do +! + if (info < 0) then + write(6,*) + write(6,*) "Error: znaupd returned info = ", info + write(6,*) "Check the documentation of znaupd for more information" + write(6,*) + stop + end if +! +! No fatal errors, post-process with ZNEUPD to extract computed eigenvalues. +! Eigenvectors may be also computed by setting RVEC = .TRUE.) +! + if (info == 1) then + write(6,*) + write(6,*) "Maximum number of iterations reached." + write(6,*) + else if (info == 3) then + write(6,*) + write(6,*) "No shifts could be applied during implicit Arnoldi update" + write(6,*) "Try increasing NCV." + write(6,*) + end if +! + rvec = .false. +! + call zneupd(rvec, 'A', select, d, v, ndim, sigma, workev, bmat, ndim, & + which, nev, tol, resid, ncv, v, ndim, iparam, ipntr, workd, & + workl, lworkl, rwork, ierr) +! + if (ierr /= 0) then + write(6,*) + write(6,*) "Error: zneupd returned ierr = ", ierr + write(6,*) "Check the documentation of zneupd for more information" + write(6,*) + stop + end if +! +! Eigenvalues are returned in the one dimensional array D and if RVEC == .TRUE. +! the corresponding eigenvectors are returned in the first NCONV == IPARAM(5) +! columns of the two dimensional array V +! + nconv = iparam(5) +! + if (rvec) then +! +! Compute the residual norm || A*x - lambda*x || for the NCONV accurately +! computed eigenvalues and eigenvectors +! + allocate(rd(ncv,3), ax(ndim)) +! + do j = 1, nconv + call matvec(a, ndim, ndim, v(:,j), ax) + call zaxpy(ndim, -d(j), v(:,j), 1, ax, 1) + rd(j,1) = real(d(j), dp) + rd(j,2) = aimag(d(j)) + rd(j,3) = dznrm2(ndim, ax, 1) / dlapy2(rd(j,1), rd(j,2)) + end do +! + call dmout(6, nconv, 3, rd, 2*nev, -6, & + "Ritz values (Real, Imag) and relative residuals") +! + deallocate(rd, ax) +! + end if +! + write(6,*) + write(6,*) "SUMMARY" + write(6,*) "=======" + write(6,*) + write(6,*) "Size of the matrix is ", ndim + write(6,*) "The number of Ritz values requested is ", nev + write(6,*) "The number of Arnoldi vectors generated (NCV) is ", ncv + write(6,*) "Portion of the spectrum: ", which + write(6,*) "The number of converged Ritz values is ", nconv + write(6,*) "The number of implicit Arnoldi update iterations taken is ", iparam(3) + write(6,*) "The number of OP*x is ", iparam(9) + write(6,*) "The convergence criterion is ", tol + write(6,*) +! + nev = nconv +! + write(6,*) "------------------ END OF ARNOLDI ITERATION ------------------" +! + deallocate(rwork) + deallocate(resid, v) + deallocate(workd, workev, workl) + deallocate(select) +! + return + end subroutine arnoldi_iteration +!==============================================================================! +! Private procedures +!==============================================================================! + subroutine matvec (a, n, lda, x, y) +! +! Compute the matrix-vector product a*x, storing the result in y +! + complex(zp), intent(in) :: a(:,:) + integer, intent(in) :: n + integer, intent(in) :: lda + complex(zp), intent(in) :: x(*) + complex(zp), intent(out) :: y(*) +! +! + call zgemv('n', n, n, one, a, lda, x, 1, zero, y, 1) +! + return + end subroutine matvec +!==============================================================================! + end module arnoldi_mod +!==============================================================================! diff --git a/src/msspec/spec/fortran/eig/ar/do_main.f b/src/msspec/spec/fortran/eig/ar/do_main.f new file mode 100644 index 0000000..b8ee9a8 --- /dev/null +++ b/src/msspec/spec/fortran/eig/ar/do_main.f @@ -0,0 +1,1558 @@ +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 + SUBROUTINE DO_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 DIM_MOD + 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 +C! READ(*,776) NFICHLEC +C! READ(*,776) ICOM +C! DO JF=1,NFICHLEC +C! READ(*,777) INDATA(JF) +C! ENDDO +C +C.......... Loop on the data files .......... +C + NFICHLEC=1 + ICOM = 5 + DO JFICH=1,NFICHLEC +C! OPEN(UNIT=ICOM, FILE=INDATA(JFICH), STATUS='OLD') + OPEN(UNIT=ICOM, FILE='../input/spec.dat', 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_AR + 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) +C! 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 +C! IF((ISOM.NE.0).OR.(NFICHLEC.EQ.1)) CLOSE(IUO1) + IF(ISOM.NE.0) CLOSE(IUO2) +CKMD STOP + return +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,E18.6,E18.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 diff --git a/src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 b/src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 new file mode 100644 index 0000000..b84f634 --- /dev/null +++ b/src/msspec/spec/fortran/eig/ar/eig_mat_ar.f90 @@ -0,0 +1,291 @@ +!==============================================================================! + subroutine eig_mat_ar (je, e_kin) +!==============================================================================! +! This subroutine stores the G_o T kernel matrix and computes a subset of the +! the eigenvalues with largest magnitude using Arnoldi iteration methods from +! ARPACK +! + use dim_mod, only: n_gaunt, nl_m + use coor_mod, only: natyp, ncorr, n_prot, sym_at + use outfiles_mod, only: outfile2 + use outunits_mod, only: iuo1, iuo2 + use trans_mod, only: lmax, vk, tl +! + use arnoldi_mod, only: arnoldi_iteration +! + implicit none +! + integer, parameter :: sp = kind(1.0) + integer, parameter :: dp = kind(1.0d0) + integer, parameter :: cp = kind((0.0,0.0)) + integer, parameter :: zp = kind((0.0d0,0.0d0)) +! +! Subroutine arguments +! + integer, intent(in) :: je + real(sp), intent(in) :: e_kin +! +! Local data +! + integer, parameter :: ibess = 3 + integer, parameter :: nprint = 10 + real(dp), parameter :: pi = acos(-1.0_dp) + real(dp), parameter :: small = 0.0001_dp + complex(zp), parameter :: ic = complex(0.0d0,1.0d0) + complex(zp), parameter :: zeroc = complex(0.0d0,0.0d0) + complex(zp), parameter :: four_pi_i = 4.0_dp*pi*ic +! + integer :: jatl, jlin, jtyp, jnum, lj, lmj, mj, nbtypj + integer :: katl, klin, ktyp, knum, lk, lmk, mk, nbtypk + integer :: j, jp, l, lin, l_max, l_min, m + integer :: ndim, n_dot, n_eig, nev, nfin, nltwo, npr, n_xmax + real(sp) :: eig, xj, yj, zj, xmax_l + real(dp) :: attkj, xkj, ykj, zkj, rkj, zdkj, krkj + complex(zp) :: expkj, sum_l, tlk +! + integer, save :: iout2, iout3 +! + real(sp), allocatable :: w1(:), w2(:) + real(dp), allocatable :: gnt(:) + complex(zp), allocatable :: hl1(:), sm(:,:), ylm(:,:), w(:) + character(:), allocatable :: outfile, path +! +! + if (je == 1) then +! +! Name of second output file where eigenvalues will be written +! + n_dot = index(outfile2, '.') + outfile = outfile2(1:n_dot)//'egv' + path = outfile2(1:n_dot)//'pth' + open(newunit=iout2, file=outfile, status='unknown') + open(newunit=iout3, file=path, status='unknown') +! + end if +! +! Construction of the multiple scattering kernel matrix G_o T. +! Elements are stored using a linear index LINJ representing (J,LJ) +! +! First compute Go T array dimension +! + jlin = 0 + do jtyp = 1, n_prot + nbtypj = natyp(jtyp) + lmj = lmax(jtyp,je) + do jnum = 1, nbtypj + do lj = 0, lmj + do mj = -lj, lj + jlin = jlin + 1 + end do + end do + end do + end do +! + ndim = jlin + write(6,*) "GoT matrix sm has dimension ", ndim +! + allocate(sm(ndim,ndim)) + sm = zeroc +! + nltwo = 2*nl_m + allocate(ylm(0:nltwo, -nltwo:nltwo)) + ylm = zeroc +! + allocate(hl1(0:nltwo)) + hl1 = zeroc +! + allocate(gnt(0:n_gaunt)) + gnt = 0.0_dp +! + 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) + do lj = 0, lmj + do mj = -lj, lj + jlin = jlin + 1 +! + klin=0 + do ktyp = 1, n_prot + nbtypk = natyp(ktyp) + lmk = lmax(ktyp,je) + do knum = 1, nbtypk + katl = ncorr(knum,ktyp) +! + if (katl /= jatl) then + xkj = real(sym_at(1,katl) - xj, dp) + ykj = real(sym_at(2,katl) - yj, dp) + zkj = real(sym_at(3,katl) - zj, dp) + rkj = sqrt(xkj*xkj + ykj*ykj + zkj*zkj) + krkj = real(vk(je), dp)*rkj + attkj = exp(-aimag(cmplx(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) + end if +! + do lk = 0,lmk + l_min = abs(lk-lj) + l_max = lk + lj + tlk = cmplx(tl(lk,1,ktyp,je)) + do mk = -lk, lk + klin = klin + 1 +! + sm(klin,jlin) = zeroc + sum_l = zeroc + if (katl /= jatl) then + call gaunt2(lk, mk, lj, mj, gnt) + do l = l_min, l_max, 2 + m = mj - mk + if (abs(m) <= l) then + sum_l = sum_l + (ic**l)*hl1(l)*ylm(l,m)*gnt(l) + end if + end do + sum_l = sum_l*attkj*four_pi_i + else + sum_l = zeroc + end if +! + sm(klin,jlin) = tlk*sum_l +! + end do + end do + end do + end do + end do + end do + end do + end do +! + deallocate(ylm, hl1, gnt) +! +! Compute subset of eigenvalues of SM using ARPACK +! +! NEV is the number of eigenvalues required, set to 2% of NDIM +! + nev = ceiling(0.02*ndim) + allocate(w(nev)) +! + call arnoldi_iteration(ndim, sm, nev, w) +! + deallocate(sm) +! +! Save results to filestream for easy access from python +! + call save_eigenvalues(w, nev, e_kin) +! +! Write results to OUTFILE on unit IOUT2 +! + write(iout2,75) + write(iout2,110) + write(iout2,80) e_kin + write(iout2,110) + write(iout2,75) + write(iout2,105) + write(iout2,75) +! + allocate(w1(nev), w2(nev)) +! + n_eig = 0 + xmax_l = 0.0 + n_xmax = 0 + do lin = 1, nev + eig = real(abs(w(lin))) + write(iout2,100) real(w(lin)), aimag(w(lin)), eig + if ((eig-xmax_l) > 0.0001) n_xmax = lin + xmax_l = max(xmax_l, eig) + w1(lin) = eig + if (eig > 1.000) then + n_eig = n_eig + 1 + end if + end do +! + write(iout2,75) + write(iout2,85) xmax_l + write(iout2,90) n_xmax + write(iout2,95) w(n_xmax) + write(iout2,75) +! +! Summarize results in main output file +! + call ordre(nev, w1, nfin, w2) +! + write(iuo1,5) + 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 >= 1) then + if (n_eig == 1) then + write(iuo1,25) ndim + else + write(iuo1,30) n_eig, ndim + end if + end if +! + write(iuo1,65) n_xmax + write(iuo1,70) w(n_xmax) + write(iuo1,10) + write(iout3,100) real(w(n_xmax)), aimag(w(n_xmax)) +! + npr = nprint/5 + write(iuo1,10) + write(iuo1,10) + write(iuo1,35) 5*npr + write(iuo1,10) + 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 == 0) then + write(iuo1,50) + else + write(iuo1,55) + end if + write(iuo1,10) + write(iuo1,10) + write(iuo1,60) +! + deallocate(w, w1, w2) +! + return +! + 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 OF A TOTAL OF ',I8,6X,'-') + 30 format(11X,'-',4X,I5,' EIGENVALUES ARE > 1 OF A TOTAL OF ',I8,2X,'-') + 35 format(11X,'-',11X,'THE ',I3,' LARGEST 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,'-',4X,'SPECTRAL RADIUS OF THE KERNEL MATRIX : ',F8.5,3X,'-') + 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,F9.5,2X,F9.5,2X,F9.5) +105 format(7X,'EIGENVALUES :',3X,'MODULUS :') +110 format(2X,'-------------------------------') +!==============================================================================! + end subroutine eig_mat_ar +!==============================================================================! diff --git a/src/msspec/spec/fortran/eig/ar/eigdif_ar.f b/src/msspec/spec/fortran/eig/ar/eigdif_ar.f new file mode 100644 index 0000000..d5fdb36 --- /dev/null +++ b/src/msspec/spec/fortran/eig/ar/eigdif_ar.f @@ -0,0 +1,103 @@ +C +C======================================================================= +C + SUBROUTINE EIGDIF_AR +C +C This subroutine computes some of the eigenvalues of the +C multiple scattering matrix using Arnoldi iteration as +C implemented in ARPACK. +C +C Last modified : 21 June 2023 +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 +ckmd IF(I_PWM.EQ.0) THEN +ckmd CALL EIG_MAT_MS(JE,ECIN) +ckmd ELSE +ckmd CALL SPEC_RAD_POWER(JE,ECIN) +ckmd ENDIF +C + call eig_mat_ar(je, ecin) +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 diff --git a/src/msspec/spec/fortran/eig/ar/main.f b/src/msspec/spec/fortran/eig/ar/main.f new file mode 100644 index 0000000..036f43e --- /dev/null +++ b/src/msspec/spec/fortran/eig/ar/main.f @@ -0,0 +1,22 @@ + SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, + & NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, + & NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, + & N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_) + + USE DIM_MOD + IMPLICIT INTEGER (A-Z) +CF2PY INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_ +CF2PY INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_ +CF2PY INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_ +CF2PY INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_ + + CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_, + & NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_, + & NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_, + & N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_) + + CALL DO_MAIN() + CALL CLOSE_ALL_FILES() + + END SUBROUTINE + diff --git a/src/msspec/spec/fortran/eig_ar.mk b/src/msspec/spec/fortran/eig_ar.mk new file mode 100644 index 0000000..85ea6a3 --- /dev/null +++ b/src/msspec/spec/fortran/eig_ar.mk @@ -0,0 +1,14 @@ +memalloc_src := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f +cluster_gen_src := $(wildcard cluster_gen/*.f) +common_sub_src := $(wildcard common_sub/*.f) +renormalization_src := $(wildcard renormalization/*.f) +#eig_common_src := $(wildcard eig/common/*.f) +eig_common_src := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f)) +eig_ar_src := $(wildcard eig/ar/*.f) +eig_ar_src_f90 := $(wildcard eig/ar/*.f90) + +SRCS = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(eig_common_src) $(eig_ar_src_f90) $(eig_ar_src) +MAIN_F = eig/ar/main.f +SO = _eig_ar.so + +include ../../../options.mk diff --git a/src/options.mk b/src/options.mk index b254cb1..1128a5d 100644 --- a/src/options.mk +++ b/src/options.mk @@ -31,7 +31,7 @@ IFORT_FFLAGS_DBG = ################################################################################ # F2PY CONFIGURATION # ################################################################################ -F2PYFLAGS = --opt=-O2 -llapack +F2PYFLAGS = --opt=-O2 -llapack -larpack F2PYFLAGS_DBG = --debug-capi --debug ################################################################################ @@ -103,7 +103,7 @@ endif FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX)) -OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS)))) +OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f90,%.o, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS))))) .PHONY: clean obj all info @@ -141,6 +141,12 @@ $(BUILDDIR)/%.o: %.f $(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION) +$(BUILDDIR)/%.o: %.f90 + @echo "Compiling $@..." + mkdir -p $(basename $@) + $(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION) + + $(SO): $(OBJS) $(MAIN_F) @echo "building Python binding $@..." mkdir -p $(BUILDDIR)