Table of Contents

GA.MUFSPEC.TEXT

This information on the muffin-tin potential program comes from a file dated 14/09/81 from Daresbury Laboratory.  It breaks down the program into the component subroutines and functions and gives descriptions of the input, output, and purpose of each.  WARNING!!!  Mufpot has changed since this was written!!!  I have made no attempt to update this information!!!

Nonetheless, I believe the information about the program and how it operates is useful and so am including it here.  The independent web pages of the Muffin-Tin Potential Program should be used for sources of current information.

Created: April 10, 1999 ---- Last Updated: April 13, 1999
By Mark D. Pauli

_______________________________________________________________

MUFFIN-TIN POTENTIAL PROGRAM
----------------------------------------------------
----------------------------------------------------
PURPOSE:
TO CALCULATE THE MUFFIN-TIN POTENTIAL FOR ATOMS IN A CRYSTAL LATTICE
ACCORDING TO THE PRESCRIPTION OF MATTHEISS (REF 1,2), STARTING WITH
THE ATOMIC WAVE FUNCTIONS (REF 3,4), AND INCLUDING THE OPTION OF A 
MADELUNG CORRECTION FOR IONIC MATERIALS.  THE EXCHANGE CONTRIBUTION IS
TAKEN PROPORTIONAL TO:
			ALPHA X (CHARGE DENSITY)**1/3
WHERE ALPHA IS AN INPUT PARAMETER.  A SAMPLING TECHNIQUE PROVIDES THE
MUFFIN-TIN ZERO: AN AVERAGE VALUE OF THE POTENTIAL OUTSIDE THE
MUFFIN-TIN SPHERE.
THE NEAREST-NEIGHBOR INFORMATION REQUIRED BY THE PROGRAM IS COMPUTED
FROM INPUT DATA OF ATOMIC POSITIONS WITHIN THE UNIT CELL.
THE EVENTUAL OUTPUT OF THE PROGRAM IS A TABULATION OF THE MUFFIN-TIN
POTENTIAL FOR EACH ATOM, IN RYDBERGS, ON LOUCKS' EXPONENTIAL GRID:
			RX(I) = EXP(-8.8+0.05(I-1))

IN ADDITION, SUBROUTINES MAY BE CALLED TO CALCULATE SCATTERING PHASE 
SHIFTS OF ATOMIC EXAFS MATRIX ELEMENTS.

REFERENCES:
1)   LOUCKS T L, (1967), AUGMENTED PLANE WAVE METHOD, BENJAMIN, NY.
2)   MATTHEISS L F, (1964), PHYS. REV. 133, A1399
3)   HERMAN F AND SKILLMAN S, (1965), ATOMIC STRUCTURE CALCULATIONS,
      PRENTICE HALL, ENGLEWOOD CLIFFS, NJ.
      WOOD J H AND BORING M, (1974)
       COMP. PHYS. COMMUNICATIONS, 7, 73
       ""     ""     ""           10,436
4)   CLEMENTI E AND RUETTI C, (1974), ATOMIC DATA AND NUCLEAR DATA
      TABLES, VOL14, 177-237, ACADEMIC.
      (EQUATIONS REFERRED TO IN THE SPEC. ARE TAKEN FROM REF 1)
5)   LIBERMAN D A, CROMER D T AND WABER J T. (1971)
      COMP. PHYS. COMMUNICATIONS, 2, 107
      ""     ""     ""            2, 471
      ""     ""     ""            9, 129
6)   HARTREE D R, (1957), THE CALCULATION OF ATOMIC STRUCTURES,
      CHAPMAN HALL, LONDON.
7)   SLATER J C, (1960), QUANTUM THEORY OF ATOMIC STRUCTURE,
      VOL 1, MC GRAW HILL, NEW YORK.

INPUT ON UNIT5, CRYSTALLOGRAPHIC DATA:
--------------

TITLE		FORMAT(9A8)			JOB TITLE
SPA		FORMAT(F8.4)			LATTICE CONSTANT (A.U.)
RC(3,3)		FORMAT(3F8.4)			R(I,J) HOLDS I'TH COORDINATE OF
						THE J'TH AXIS OF UNIT CELL.
IPX, IREL	FORMAT(2I4)			SWITCH FOR EXAFS ROUTINE (IPX=1)
						SWITCH FOR L.S TERM IN PHASE
NR		FORMAT(I4)			NUMBER OF INEQUIVALENT ATOMS
						(ATOMIC TYPES) IN UNIT CELL.
FOR EACH ATOMIC TYPE:				---------------------------------
NAME		FORMAT(2A8)							I
NRR,Z,ZC,RMT	FORMAT(I4,3F8.4)		NUMBER OF ATOMS IN UNIT CELL,   I
						ATOMIC NUMBER, VALENCE,         I
						MUFFIN-TIN RADIUS (A.U.)        I
RK(3,NRR)	FORMAT(3F8.4)			CARTESIAN COORDINATES OF ATOMS  I
						IN UNIT CELL, IN UNITS OF SPA,  I
						RK(I,J) HOLDS I'TH COORDINATE   I
						OF J'TH ATOM IN THE UNIT CELL   I
						---------------------------------
NHM, ALPHA	FORMAT(I4,F8.4)			MUFFIN-TIN ZERO PARAMETER,
						VALUE OF ALPHA FOR EXCHANGE.
						THE MUFFIN-TIN ZERO IS ESTIMATED
						FROM A SAMPLING GRID WHOSE MOST
						ACCURATE MESH IS GIVEN BY NHM**3
						POINTS SPANNING THE UNIT CELL.
						IF NHM=0 THE MUFFIN-TIN ZERO
						VINT IS GIVEN BY ALPHA.
FOR EACH ATOMIC TYPE:			--------------------------------------------
L,NE,E,DE	FORMAT(2I4,2F8.4)		L IS MAXIMUM ANGULAR MOMENTUM FOR  I
						PHASE SHIFT ROUTINE OR INITIAL     I
						STATE MOMENTUM FOR EXAFS ROUTINE,  I
						NUMBER OF ENERGY POINTS, INITIAL   I
						VALUES AND INCREMENT FOR:          I
						ENERGY AT WHICH PHASE SHIFT IS     I
						CALCULATED, OR ENERGY OF FINAL     I
						STATE IN EXAFS CALCULATION (BOTH   I
						IN HARTREES RELATIVE TO MUFFIN-    I
						TIN ZERO LEVEL).  IF IPX IS NOT    I
						UNITY AND L IS LESS THAN ZERO      I
						NO CALCULATION IS PERFORMED FOR    I
						THAT ATOM.                         I
FOR CALLS TO EXAFS ROUTINE RONUC:                                                  I
NCORE, ECORE	FORMAT(I4,F11.5)		INDEX FOR INITIAL STATE WAVE       I
						FUNCTION, INITIAL STATE ENERGY     I
						IN HARTREES RELATIVE TO VACUUM,    I
						IF NCORE IS NEGATIVE OR ZERO NO    I
						CALCULATION IS PERFORMED FOR THE   I
						ATOM.                              I
						------------------------------------




INPUT ON UNIT4, ATOMIC STRUCTURE DATA
--------------

FOR EACH ATOMIC TYPE:

EITHER
1) WAVE FUNCTIONS ON HERMAN-SKILLMAN TYPE RADIAL MESH, VIA ROUTINE HSIN
WFN1		FORMAT(A8)			'HERMAN-S'
NAME,Z,NC	FORMAT(2A8/F9.4/I4)		NAME,ATOMIC NUMBER,
						NUMBER OF ATOMIC STATES
FOR EACH ATOMIC STATE:
LC,N,FRAC	FORMAT(I4/I4/F9.4)		ANGULAR MOMENTUM, NUMBER OF
						TABULATION POINTS, FRACTIONAL
						OCCUPATION
RS		FORMAT(1P5E14.7)		ATOMIC WAVE FUNCTIONS ON H-S MESH

OR
2) WAVE FUNCTIONS IN THE CLEMENTI PARAMETERISED FORM, VIA ROUTINE CLEMENTI
WFN2		FORMAT(A8)			'CLEMENTI'
NAME		FORMAT(2A8)			
IPRINT		FORMAT(I4)			.NE.0 FOR PRINTOUT OF WAVE
						FUNCTIONS
NC		FORMAT(I4)			NUMBER OF ATOMIC STATES
NS		FORMAT(I4)			NUMBER OF CLEMENTI COEFFICIENTS
NT(NS),EX(NS)	FORMAT(I4,F11.5)		CLEMENTI PARAMETERS: ORBITAL
						EXPONENTS
FOR EACH ATOMIC STATE:
LC		FORMAT(I4)			ANGULAR MOMENTUM
FAC(NS)	FORMAT(5F11.5)		CLEMENTI PARAMETERS: EXPANSION
						COEFFICIENTS
FRAC		FORMAT(F11.5)			FRACTIONAL OCCUPATION

OR
3) ATOMIC CHARGE DENSITY ON RADIAL GRID, CALCULATED IN A RELATIVISTIC
    SELF-CONSISTENT FIELD APPROACH AND INPUT VIA ROUTINE HEXIN
WFN3		FORMAT(A8)			'HEXINPUT'
NAME,N,J,RN,H	FORMAT(2A8/2I5,		NAME, NUMBER OF CHARGE DENSITY
			  F6.3,F4.0)		GRID POINTS, NUMBER OF ATOMIC
						STATES, MAXIMUM GRID RADIUS, GRID
						INTERVAL (RECIPROCAL)
RS		FORMAT(1P5E14.7)		CHARGE DENSITY ON RADIAL GRID

OR
4) MUFFIN-TIN POTENTIAL ON RADIAL GRID, CORRECTED TO MUFFIN-TIN LEVEL
WFN4		FORMAT(A8)			'POTENTIA'

FOR EACH ATOMIC TYPE:
NTAB		FORMAT(I4)			NUMBER OF TABULATED POINTS
WK1,WK2	FORMAT(5E14.6)		RADIAL GRID AND TABULATED 
						POTENTIAL, CORRECT TO M.T. ZERO


OUTPUT ON UNIT 6
----------------

TITLE		FORMAT(9A8)			JOB TITLE
RC(3,3)		FORMAT(3F8.4)			AXES OF UNIT CELL IN A.U.
						(MULTIPLIED BY SPA)
AV,OMA,RWS	FORMAT(F15.4/			UNIT CELL VOLUME, ATOMIC
			F18.4/F12.4)		VOLUME, WIGNER-SEITZ RADIUS
FOR EACH ATOMIC TYPE:				------------------------------------
IR,NAME,NRR	FORMAT(I2,2AI/IR)		ATOM TYPE NUMBER, NAME, NUMBER
						OF ATOMS OF TYPE IR PER UNIT CB
Z,ZC,RMT	FORMAT(F15.1/F21.1,		ATOMIC NUMBER, VALENCE CHARGE,
			/F14.4)			AND MUFFIN-TIN RADIUS
						------------------------------------
RX		FORMAT(5(10F11.5/))		LOUCK'S RADIAL MESH
						------------------------------------
FOR EACH ATOMIC TYPE:

IF HSIN OR CLEMIN CALLED:

NAME		FORMAT(2A8)

FOR EACH ATOMIC STATE:
LC,WFC		FORMAT(I3,5(10F11.5/))		ANGULAR MOMENTUM, WAVEFUNCTION
						X RADIUS ON RADIAL GRID
						ON RADIAL GRID
RX(NX),NX,RHO	FORMAT(F12.5,I4//,		MAXIMUM RADIAL COORDINATE,
			5(10E12.4/))		NUMBER OF GRID POINTS,
						CHARGE DENSITY ON RADIAL GRID

IF HEXIN CALLED:

NAME,RX(NX),	FORMAT(2A8,F12.4,I4/		NAME, MAXIMUM RADIAL COORDINATE
NX,RHO			5(10E12.4/))		NUMBER OF GRID POINTS, CHARGE
						DENSITY ON RADIAL GRID
						-------------------------------------

FOR ALL INPUT OPTIONS EXCEPT WFN4:		*************************************

FOR EACH ATOMIC TYPE:				-------------------------------------
NAME,ZE,ZNT,	FORMAT(2A8,F12.5,F7.1,		NAME, INTEGRATED ELECTRONIC
RX(NIX),NIX		,F12.5,I4)		CHARGE, NORMALISED CHARGE,
						MAXIMUM RADIAL COORDINATE AND
						INDEX
SIG		FORMAT(5(10E12.4/))		HARTREE POTENTIAL FOR ISOLATED
						ATOM OBTAINED BY SOLVING
						POISSON'S EQUATION (RYDBERGS)

AD(J,IR)	FORMAT(15F8.4)			DISTANCE TO J'TH RADIAL SHELL
						FOR ATOM TYPE IR
NA(J,IR)	FORMAT(15I5)			NUMBER OF ATOMS IN J'TH SHELL
IA(J,IR)	FORMAT(15I5)			ATOM TYPE IN J'TH SHELL
						--------------------------------------
ALPHA,NHM	FORMAT(F10.4//I4)		SLATER EXCHANGE PARAMETER, MUFFIN
						TIN ZERO PARAMETER
VOLC,VINT	FORMAT(F11.4,F12.4)		INTERMEDIATE MONTE CARLO ESTIMATE
						OF INTERSTITAL VOLUME AND MUFFIN TIN
						ZERO IN RYDBERGS.
NPOINT,		FORMAT(I10/F11.4,		TOTAL NUMBER OF MONTE CARLO MESH
VOLT,VOLC		F11.4)			POINTS, TRUE INTERSTITIAL VOLUME,
						MONTE CARLO VOLUME
VHAR,VEX,VINT	FORMAT(F14.5,F12.5/		AVERAGE HARTREE POTENTIAL, AVERAGE
			F12.5)			EXCHANGE POTENTIAL, MUFFIN-TIN
						ZERO LEVEL (ALL IN RYDBERGS).

FOR IONIC MATERIALS WHERE MADELUNG CORRECTION IS MADE:--------------------------------

G(I,J)		FORMAT(3F8.4)			I'TH COORDINATE OF THE J'TH
						RECIPROCAL LATTICE VECTOR
RCMIN,GMIN,	FORMAT(F10.4,F10.4,		SMALLEST LATTICE VECTOR, SMALLEST
RKMAX,TEST,AL		F10.4,E12.4/E12.4)	RECIPROCAL LATTICE VECTOR, LARGEST
						ATOMIC COORDINATE, ACCURACY,
						SEPARATION CONSTANT
ITR		FORMAT(I3)			REAL SPACE SUMMATION PARAMETER
VMM(IR)		FORMAT(5E12.4)			INTEGRAL OF POTENTIAL ABOUT ATOM
						TYPE IR OUT TO MUFFIN-TIN RADIUS
ITG		FORMAT(I3)			RECIPROCAL SPACE SUMMATION
						PARAMETER
VMM(IR)		FORMAT(5E12.4)			INTEGRAL OF POTENTIAL ABOUT ATOM
						TYPE IR OUT TO MUFFIN-TIN RADIUS
AMT		FORMAT(5E12.4)			MADELUNG CORRECTION TO MUFFIN-TIN
						ZERO LEVEL IN RYDBERGS
						--------------------------------------

FOR EACH ATOMIC TYPE:				--------------------------------------
NAME,RMT	FORMAT(2A8,F8.4)		NAME, MUFFIN-TIN RADIUS
RX,VH,VS,	FORMAT(F12.5,4E20.6)		RADIAL MESH, HARTREE POTENTIAL,
VMAD,SIG					EXCHANGE POTENTIAL, MADELUNG
						CORRECTION, TOTAL POTENTIAL
						ALL IN RYDBERGS
						--------------------------------------

FOR CALL TO EXAFS ROUTINE (IPX=1), WITH INPUT OPTIONS WFN3 AND WFN4:
FOR EACH ATOMIC TYPE:
S,FNNS		FORMAT(E10.4,E14.5)		SCATTERING VECTOR, SCATTERING
						FACTOR
						**************************************




MAIN PROGRAM
------------

	DIMENSION SIG(260,2),RHO(260,2),VH(260,2),VS(260,2),
	+ VMAD(260,2),RX(260)
	DIMENSION RC(3,3),RK(3,12),ZM(12),Z(2),ZC(2),RMT(2),
	+ NRR(2),NCON(2),NX(2)
	DIMENSION IA(15,2),NA(15,2),AD(15,2)
	DIMENSION ATL(260),XM(260)
	DIMENSION CX(260),CRHO(260),CRI(260)
	DIMENSION RRHO(260),FNNS(260,2)
	REAL*8 TITLE(9),NAME(2,2),WFN,WFN1,WFN2,WFN3,WFN4/
	DATA NGRID,MC,PI/250,15,3.1415926536/,WFN1,WFN2,WFN3,WFN4
	+ 'HERMAN-S','CLEMENTI','HEXINPUT','POTENTIA'/
	COMMON /WK/ WK1(260),WK2(260)
	COMMON /WF/ WF2(260,15,2),WC(15),LC(15)
	COMMON /REL/ IREL,ALDS,VBCON

AFTER INPUT OF THE CRYSTALLOGRAPHIC DATA ON UNIT 5, ATOMIC CHARGE
DENSITIES, RHO, ARE RETURNED FROM HSIN, CLEMIN OR HEXIN:  THESE ARE
NORMALISED TO THE TOTAL ELECTRONIC CHARGE AROUND EACH ATOM.  ROUTINE
POISSON SOLVES POISSON'S EQUATION (EQ 3.24 IN LOUCKS (1967)) TO GIVE
THE HARTREE POTENTIAL IN RYDBERGS FOR AN ISOLATED NEUTRAL ATOM IN 
ARRAY SIG.

SUBROUTINE NBR COMPUTES THE NEAREST-NEIGHBOUR INFORMATION FOR EACH
ATOMIC TYPE, WHICH IS NEEDED FOR THE SUMMATION PROCEDURE CARRIED OUT
IN SUMAX (EQ 3.26,3.28).  ELECTROSTATIC AND EXCHANGE CONTRIBUTIONS TO
THE POTENTIAL ARE LEFT IN ARRAYS VH AND VS.

SUBROUTINE MTZ SAMPLES POINTS IN THE INTERSTITIAL REGION OF THE UNIT
CELL, BETWEEN MUFFIN-TIN SPHERES, AND PROVIDES AN AVERAGE VALUE OF THE 
POTENTIAL IN THIS REGION, WHICH BECOMES THE MUFFIN-TIN ZERO LEVEL IN
RYDBERGS.
DEPENDING ON THE FINENESS OF THE SAMPLING GRID, CONTROLLED BY INTEGER
NHM, THIS CAN BE THE MOST TIME-CONSUMING PART OF THE PROGRAM.

IF THE VALENCE OF ALL ATOMS IS NOT ZERO, THEN ROUTINE MAD IS CALLED.
THE VALENCE CORRECTION, ZC, IS INCLUDED BY CALCULATING A SPHERICALLY
SYMMETRIC POTENTIAL TO REPRESENT A POINT CHARGE AT EACH ATOMIC POSITION
AND SUMMING THE RESULTING LATTICE OF POINT CHARGES BY A MADELUNG TYPE
SUMATION (EWALD SPLIT TECHNIQUE).  NOTE THAT TC IS THE VALENCE CHARGE 
OR IONICITY: >0 FOR A POSITIVE ION, <0 FOR A NEGATIVE ION.  THE MADELUNG
CORRECTION IS ITSELF CORRECTED BY A MUFFIN-TIN ZERO ESTIMATE, AND IS THEN
COMBINED WITH THE POTENTIAL ALREADY OBTAINED FOR THE NEUTRAL ATOM.

HARTREE + EXCHANGE + MADELUNG CORRECTION ARE FINALLY COMBINED TO 
PRODUCE THE MUFFIN-TIN POTENTIAL IN RYDBERGS, IN ARRAY SIG.  THESE THREE
TOGETHER WITH THE TOTAL POTENTIAL, ARE PRINTED OUT CORRECT TO THE MUFFIN-
TIN ZERO LEVEL.

IF THE EXAFS SUBROITINE RONUC IS CALLED, THE X-RAY SCATTERING FACTOR
IS TABULATED AS A FUNCTION OF THE SCATTERING VECTOR.

DEPENDING ON THE CHOICE OF OPTION THE PROGRAM CALLS ROUTINES PHCALC
OR RONUC TO CALCULATE RESPECTIVELY THE SCATTERING PHASE SHIFTS
OR ATOMIC EXAFS MATRIX ELEMENTS.  DETAILS OF THESE ROUTINES ARE
GIVEN BELOW.

ROUTINES CALLED:  HSIN, CLEMIN, HEXIN, CHGRID, POISON, NBR, SUMAX,
		  MTZ, MAD, FTSIN, PHSCALC, RONUC




SUBROUTINE HSIN
---------------

	SUBROUTINE HSIN(RHO,RX,NX,NGRID,IR)
	DIMENSION RHO(NGRID),RX(NGRID)
	REAL*8 NAME(2)
	COMMON /WK/ RR(441),RS(441),PS(260)
	COMMON /WF/ WFC(260,15,2),WC(15),LC(15)

DEALS WITH INPUT OF ATOMIC WAVEFUNCTIONS ON A HERMAN-SKILLMAN TYPE
RADIAL MESH FROM PROGRAM HERSKL, AND PRODUCES CHARGE DENSITY RHO
TABULATED ON MESH RX.  (SEE HERMAN-SKILLMAN (1965), WOOD AND BORING
(1974), AND NOTE THAT THE TABULATION IS WAVE FUNCTION X RADIUS)
THE HERMAN-SKILLMAN MESH USED HAS AN INITIAL INCREMENT:
	DR=0.0025*0.88534138/EXP(ALOG(Z)/3.0)
AND A DOUBLING PERIOD OF 40 MESH POINTS.

INPUT:	RX	RADIAL MESH (I.E.CLEMENTI)
	NGRID	DIMENSION OF RX
	IR	ATOM TYPE

OUTPUT: RHO	4 X PI X CHARGE DENSITY X (RADIUS SQUARED)
	NX	LIMIT OF TABULATION FOR RHO

NOTE THAT THE CORE WAVEFUNCTIONS STORED IN WFC IN COMMON /WF/ ARE
USED BY ROUTINE RONUC IN CALCULATION OF EXAFS MATRIX ELEMENTS.

ROUTINE CALLED:  CHGRID



SUBROUTINE CLEMIN
-----------------

	SUBROUTINE CLEMIN(RHO,RX,NX,NGRID,IR)
	DIMENSION RHO(NGRID),RX(NGRID)
	REAL*8 NAME(2)
	COMMON /WK/ EX(20),FAC(20),FNT(20),NT(20)
	COMMON /WF/ WFC(260,15,2),WC(15),LC(15)

DEALS WITH INPUT OF ATOMIC WAVEFUNCTIONS IN THE CLEMENTI PARAMTERISED
FORM, AND PRODUCES CHARGE DENSITY RHO ON RADIAL MESH RX.  (SEE CLEMENTI
AND ROETTI (1974))
AS WITH HSIN TABULATION IS WAVE FUNCTION X RADIUS.

INPUT:	RX	RADIAL MESH (I.E. CLEMENTI)
	NGRID	DIMENSION OF RX
	IR	ATOM TYPE

OUTPUT:	RHO	4 X PI X CHARGE DENSITY X (RADIUS SQUARED)
	NX	LIMIT OF TABULATION FOR RHO

NOTE THAT THE CORE WAVEFUNCTIONS STORED IN WFC IN COMMON /WF/ ARE
USED BY ROUTINE RONUC IN CALCULATION OF EXAFS MATRIX ELEMENTS.



SUBROUTINE HEXIN
----------------

	SUBROUTINE HEXIN(RHO,RX,NX,NGRID)
	DIMENSION RHO(NGRID),RX(NGRID),RR(421),RS(421)
	REAL*8 NAME(2)

DEALS WITH INPUT OF ATOMIC CHARGE DENSITY FROM RELATIVISTIC SELF
CONSISTENT FIELD CALCULATION CARRIED OUT IN PROGRAM HEX (SEE LIBERMAN,
CROMER AND WABER (1971)).  PRODUCES CHARGE DENSITY RHO ON RADIAL MESH

INPUT:	RX	RADIAL MESH (I.E. CLEMENTI)
	NGRID	DIMENSION OF RX
	IR	ATOM TYPE

OUTPUT:	RHO	4 X PI X CHARGE DENSITY X (RADIUS SQUARED)
	NX	LIMIT OF TABULATION FOR RHO

NOTE THAT THE CORE WAVEFUNCTIONS STORED IN WFC IN COMMON /WF/ ARE
USED BY ROUTINE RONUC IN CALCULATION OF EXAFS MATRIX ELEMENTS.

ROUTINE CALLED:  CHGRID



SUBROUTINE CHGRID
-----------------

	SUBROUTINE CHGRID(FX,X,NX,FY,Y,NY)
	DIMENSION FX(NX),X(NX),FY(NY),Y(NY)

PERFORMS INTERPOLATION FROM ONE TABULATION GRID TO ANOTHER, VIA A
PIECEWISE QUADRATI PROCEDURE USING AITKEN'S DIVIDED DIFFERENCE
SCHEME.

INPUT:	FX	FUNCTION TABULATED ON GRID X
	X	INPUT TABULATION GRID 
	NX 	NUMBER OF TABULATION POINTS
	Y	OUTPUT TABULATION GRID
	NY DIMENSION OF NY

OUTPUT:	FY	FUNCTION TABULATED ON GRID Y
	NY	REVISED NUMBER OF TABULATION POINTS



SUBROUTINE POISON
-----------------

	SUBROUTINE POISON(PSQ,Z,J,W)
	DIMENSION PSQ(J),W(J)
	COMMON /WK/ E(260),F(260)

SOLVES POISSON'S EQUATION TO PRODUCE COULOMB POTENTIAL IN RYDBERGS FROM
CHARGE DENSITY - SEE EQ 3.24 AND APPENDIX 1 OF LOUCKS, REF 1.  NOTE THAT
THIS ROUTINE IS SPECIALISED FOR LOUCKS' GRID RX.

INPUT:	PSQ	4 X PI X (RADIUS SQUARED) X (CHARGE DENSITY)
	Z	ATOMIC NUMBER (NUCLEAR CHARGE)
	J	LIMIT OF TABULATION

OUTPUT:	W	SQRT(RADIUS) X ELECTROSTATIC POTENTIAL (RYDBERGS)



SUBROUTINE NBR
--------------

	SUBROUTINE NBR(IA,NA,AC,NCON,NRR,NR,RC,RK,N,RMAX,MC)
	DIMENSION IC(MC,NR),NA(MC,NR),AD(MC,NR),NCON(NR),NRR(NR),
	+ RC(3,3),RK(3,N)
	COMMON /WK/ RJ(3)

PRODUCES NEAREST-NEIGHBOR INFORMATION FOR USE IN SUMMING CONTRIBUTIONS
FROM NEIGHBORING ATOMS IN THE CRYSTAL LATTICE.  UP TO MC SHELLS OF 
NEIGHBORS ARE INCLUDED AROUND EACH ATOMIC TYPE: EACH SHELL CONTAINS
ATOMS OF ONLY ONE ATOMIC TYPE.

INPUT:	NRR(IR)	NUMBER OF TYPE IR ATOMS IN UNIT CELL
	NR	NUMBER OF ATOMIC TYPES IN UNIT CELL
	RC(I,J)	I'TH COORDINATE OF J'TH AXIS OF THE UNIT CELL
	RK(I,J)	I'TH COORDINATE OF J'TH ATOM IN THE UNIT CELL
	N	TOTAL NUMBER OF ATOMS IN UNIT CELL
	RMAX	NEAREST NEIGBOR SHELLS ARE INCLUDED OUT TO A
		RADIUS RMAX
	MC	MAXIMUM NUMBER OF SHELLS ALLOWED (SET TO 15 IN
		MAIN PROGRAM DATA STATEMENT)

OUTPUT:	FOR A TYPE IR ATOM
	IA(I,IR)	ATOMIC TYPE IN I'TH NEIGHBORING SHELL
	NA(I,IR)	NUMBER OF ATOMS IN I'TH NEIGHBORING SHELL
	AD(I,IR)	RADIUS OF I'TH SHELL



SUBROUTINE SUMAX
----------------

	SUBROUTINE SUMAX(CHIAT, CHICRY,NX,NEIGHB,IA,NA,AD,JTOP,
	+ NGRID,NR)
	DIMENSION CHIAT(NGRID,NR),CHICRY(NGRID),NX(NR),IA(1),NA(1),
	+ AD(1)

PERFORMS THE SUMMATION OF CONTRIBUTIONS TO CHARGE DENSITY AND COULOMB
POTENTIAL ABOUT A PARTICULAR ATOM, FROM NEIGHBORING ATOMS IN THE 
CRYSTAL LATTICE (EQS 3.22,3.26,3.28 OF LOUCKS (1967)).  THIS ROUTINE IS
MODIFIED FROM APPENDIX 2 OF LOUCKS, TO INCLUDE MORE THAN ONE ATOMIC
TYPE IN THE UNIT CELL.  SPECIALISED FOR LOUCKS' TABULATION GRID.

INPUT:	CHIAT		POTENTIAL OR CHARGE DENSITY TABULATED ABOUT
			EACH ATOMIC TYPE
	NX 		LIMIT OF TABULATION FOR EACH TYPE
	NEIGHB		NUMBER OF NEAREST-NEIGHBOR SHELLS INCLUDED
	IA,NA,AD	NEARREST-NEIGHBOR DATA, FROM ROUTINE NBR
	JTOP		TABULATION GRID POINT FOR MUFFIN-TIN RADIUS
	NGRID		DIMENSION OF CHIAT
	NR		NUMBER OF ATOMIC TYPES

OUTPUT:	CHICRY		TABULATION OF TOTAL POTENTIAL OR CHARGE DENSITY
			ABOUT THE ATOM UNDER CONSIDERATION



SUBROUTINE MTZ
--------------

	SUBROUTINE MTZ(NHM,SIG,RHO,RX,NGRID,RMT,NRR,NX,NR,
	+ RC,RK,N,VINT,ALPHA,AV)
	DIMENSION SIG(NGRID,NR),RHO(NGRID,NR),RX(NGRID),RMT(NR)
	+ NRR(NR),NX(NR),RC(3,3),RK(3,N)
	COMMON /WK/ X(3),RB(3)
	DATA PI /3.14159265358/

PROVIDES AN ESTIMATE OF THE AVERAGE VALUE OF THE POTENTIAL IN THE
REGION OUTSIDE THE MUFFIN-TIN SPHERES, AS A MEASURE OF THE MUFFIN-TIN
ZERO LEVEL.  THE UNIT CELL IS SPANNED BY A 3-D MESH OF SAMPLE POINTS:
EACH POINT LYING WITHIN THE INTERSTITIAL REGION IS INCLUDED IN THE
AVERAGING PROCEDURE.
THE FINENESS OF THE SAMPLING GRID IS CONTROLLED BY A PARAMETER WH
WHICH IS INCREASED PROGRESSIVELY FROM 1 TO A MAXIMUM VALUE NHM, AND
THE UNIT CELL IS SPANED BY NH**3 SAMPLING POINTS.
FOR EACH VALUE OF NH, THE CURRENT MONTE-CARLO VOLUME AND MUFFIN-TIN
ZERO ESTIMATES ARE PRINTED OUT TO INDICATE THE DEGREE OF CONVERGENCE.
IN ADDITION THE TRUE INTERSTITIAL VOLUME IS PRINTED.
THE CALCULATION IS IN RYDBERGS.

INPUT:	NHM		MAXIMUM VLUE OF MONTE CARLO GRID PARAMETER NH
	SIG(I,IR)	I'TH TABULATED VALUE OF COULOMB POTENTIAL
			ABOUT A TYPE IR ATOM (RYDBERGS)
	RHO(I,IR)	I'TH TABULATED VALUE O CHARGE DENSITY ABOUT A
			TYPE IR ATOM
	RX		TABULATION GRID
	NGRID		DIMENSION OF SIG, RHO
	RMT(IR)		MUFFIN-TIN RADIUS FOR TYPE IR ATOM
	NRR(IR)		NUMBER OF TYPE IR ATOMS IN UNIT CELL
	NX(IR)		LIMIT OF TABULATION ABOUT TYPE IR ATOM
	NR		NUMBER OF ATOMIC TYPES
	RC(I,J)		I'TH COORDINATE OF J'TH AXIS OF UNIT CELL
	RK(I,J)		I'TH COORDINATE OF J'TH ATOM IN UNIT CELL
	N		TOTAL NUMBER OF ATOMS IN UNIT CELL
	ALPHA		EXCHANGE PARAMETER ALPHA
	AV		ATOMIC VOLUME (TOTAL VOLUME OF UNIT CELL)

OUTPUT:	VINT		MUFFIN-TIN ZERO LEVEL (RYDBERGS)



SUBROUTINE MAD
--------------

	SUBROUTINE MAD(VMAD,RX,NGRID,RMT,NRR,NX,NR,RC,RK,ZM,N,AV)
	DIMENSION VMAD(NGRID,NR),RX(NGRID),RC(3,3),RK(3,N),ZM(N),
	+ RMT(NR),NRR(NR),NX(NR)
	COMMON /WK/ G(3,3),VMM(5),FR(5),RA(3),GA(3)
	DATA PI,TEST/3.1415926526,1.0E-4/

CALCULATES THE SPHERICALLY AND SPATIALLY AVERAGED FIELDS FROM A LATTICE
OF POINT CHARGES AT THE POSITIONS OF THE ATOMS IN THE UNIT CELL.
FOR EACH ATOMIC TYPE THE VALUE OF THE TOTAL ELECTROSTATIC POTENTIAL:
			SIGMA J ( ZJ / MOD(R-RJ) )
FOR ALL THE ATOMS J OF THE LATTICE, IS TABULATED OUT TO THE MUFFIN-TIN
RADIUS IN VMAD, AND CORRECTED BY AN ESTIMATE OF THE MUFFIN-TIN ZERO
LEVEL.  THE SUMMATION IS ACHIEVED BY SPLITTING INTO REAL AND RECIPROCAL
SPACE CONTRIBUTIONS (EWALD SPLIT), THE SEPARATION CONSTANT BEING
CHOSEN TO MAINTAIN APPROXIMATELY EQUAL NUMBERS OF TERMS IN REAL AND
RECIPROCAL SUMATIONS.  THIS IS DONE ON THE BASIS OF A CONVERGENCE
CRITERION TEST.  SET TO 1.0E-4.

INPUT:	RX		RADIAL TABULATION GRID
	NGRID		NUMBER OF TABULATION POINTS
	RMT(IR)		MUFFIN-TIN RADIUS FOR TYPE IR ATOM
	NRR(IR)		NUMBER OF TYPE IR ATOMS IN UNIT CELL
	NR		NUMBER OF ATOMIC TYPES
	RC(I,J)		I'TH COORDINATE OF J'TH AXIS OF UNIT CELL
	RK(I,J)		I'TH COORDINATE OF J'TH ATOM IN UNIT CELL
	ZM(J)		CHARGE ON THE J'TH ATOM IN THE UNIT CELL
	N		TOTAL NUMBER OF ATOMS IN UNIT CELL
	AV		UNIT CELL VOLUME

OUTPUT:	VMAD(I,IR)	I'TH TABULATED VALUE OF THE TOTAL POTENTIAL
			ABOUT TYPE IR ATOM (RYDBERGS)



SUBROUTINE FTSIN
----------------

	SUBROUTINE FTSIN(FUN,R,NR,A,FTI,NS)
	DIMENSION F(4),P(4),D(4),FL(4),FR(4),CL(4,4),CR(4,4)
	DIMENSION FUN(260),R(260),SC(2)

FOURIER TRANSFORMS INPUT FUNCTION FUN, SUPPORTED BY ROUTINE CP.
FUN IS READ IN AT GRID POINTS GIVEN BY ARRAY R, AND IS TRANSFORMED
WITH RESPECT TO SIN(A*R) WHERE A IS A SIMPLE INPUT VARIABLE.

INPUT:	FUN		INPUT FUNCTION TO BE TRANSFORMED
	R		INPUT GRID
	NR		NUMBER OF INPUT GRID POINTS
	A		TRANSFORM VARIABLE
	NS		PARAMETER FIXING THE DEGREE OF POLYNOMIAL USED TO
			APPROXIMATE FUN, AND HENCE THE NUMBER OF SUB-
			INTERVALS WITHIN THE INTERVAL OF THE R AXIS

OUTPUT:	FTI		FOURIER TRANSFORM OF FUN W.R.T. SIN(A*R)

ROUTINE CALLED: CP



SUBROUTINE CP
-------------

	SUBROUTINE CP(P,D,CL,CR,NS,M)
	DIMENSION P(NS),D(NS),CL(M,NS),CR(M,NS)

SUPPORTS ROUTINE FTSIN IN CALCULAION OF POLYNOMIAL COEFFICIENTS







PHASE SHIFT SUBROUTINE: PHCALC
------------------------------

	SUBROUTINE PHCALC(Z,NM,VM,XM,WM,RM,E1,DE,NE,LMAX)
	DIMENSION XM(260),VM(260),Y(200),VL(200)
	DIMENSION STORE(200),PHS(50)
	COMMON /WS/ KCF,JRCF,Y,CL,E,VL
	COMMON /REL/ IREL,ALDS,VBCON
	DATA ET/3.0/

PURPOSE:

TO COMPUTE PHASE SHIFTS BY INTEGRATION OF THE SCHRODINGER EQUATION
FOR THE MUFFIN-TIN POTENTIALS CALCULATED IN THE MAIN PROGRAM.
(SEE LOUCKS (1967) P56 AND FF)

ROUTINE INTERPOLATES INPUT (MUFFIN-TIN POTENTIAL IN HARTREES X RADIAL
DISTANCE ) ON TO A HERMAN-SKILLMAN GRID AND THEN CALLS ROUTINE CACLFI TO
CALCULATE PHASE SHIFTS USING THE RUNGE-KUTTA-MERSON SUBROUTINE DA01A, THE
BESSEL FUNCTION SUBROUTINE CALCBF AND THE DIFFERENCTIAL FUNCTION ROUTINE
DYBDX.

THE MUFFIN-TIN POTENTIAL ON THE H-S GRID IS OUTPUT ON UNIT 7 AND THE 
PHASE SHIFTS ARE OUTPUT ON UNIT 8.

IF THE OPTION IS CALLED FOR (IREL=1), AN L.S TERM IS INCLUDED AND
PHASE SHIFTS ARE PRINTED FOR SPIN UP AND SPIN DOWN CASES.  THE L.S
TERM IS CUT OFF AT THE SECOND MESH POINT AND SET TO THE VALUE AT
THAT POINT TO AVOID DIVERGENCES.

INPUT:	Z		ATOMIC NUMBER
	NM		NUMBER OF DATA POINTS ON INPUT GRID
	VM		INPUT FUNCTION: MUFFIN-TIN POTENTIAL (HARTREES)
					X RADIAL DISTANCE
	XM		INPUT GRID: LOUCKS' RADIAL GRID
	WM		MUFFIN-TIN ZERO LEVEL IN HARTREES
	RM		MUFFIN-TIN RADIUS
	E1		INITIAL ENERGY FOR PHASE SHIFT CALCULATION
	DE		ENERGY INCREMENT
	NE		NUMBER OF ENERGY VALUES
	LMAX		MAXIMUM ANGULAR MOMENTUM FOR WHICH PHASE SHIFT
			IS REQUIRED

OUTPUT ON UNIT 6:
----------------

Z,NM,WM,RM	FORMAT(F5.1,I4,		ATOMIC NUMBER, NUMBER OF INPUT
			F7.4,F7.4)	GRID POINTS, MUFFIN-TIN ZERO LEVEL
					IN HARTREES, MUFFIN-TIN RADIUS
XM,VM		FORMAT(200(/F12.8,	LOUCKS' RADIAL GRID, MUFFIN-TIN
			F14.5))		POTENTIAL IN HARTREES X RADIAL
					DISTANCE ON INPUT GRID XM
ET		FORMAT(F8.4)		HERMAN-SKILLMAN GRID PARAMETER
Y,VL		FORMAT(6E12.4)		HERMAN-SKILLMAN GRID, MUFFIN-TIN
					POTENTIAL IN HARTREES X RADIAL
					DISTANCE ON GRID Y
E,PHS		FORMAT(F12.5,		ENERGY IN HARTREES, PHASE SHIFTS
			10(/5F12.5))	FOR L=0 TO L=LMAX

OUTPUT ON UNIT 7:
----------------

NMESH		FORMAT(I4)		NUMBER OF MESH POINTS IN H-S GRID
Y,VL		FORMAT(5E14.6)		HERMAN-SKILLMAN GRID AND MUFFIN-
					TIN POTENTIAL IN HARTREES ON GRID
					Y

OUTPUT ON UNIT 8:
----------------

E,PHS		FORMAT(F12.5,		ENERGY IN HARTREES, PHASE SHIFTS
			10(/5E12.5))	FOR L=0 TO L=LMAX

ROUTINES CALLED:  CHGRID,CALCFI



SUBROUTINE CALCFI
-----------------

	SUBROUTINE CALCFI(L,Z,STORE,DEL,NF,RHO)
	DIMENSION ATP(200),X(200),STORE(200)
	DIMENSION A(2),B(2),C(2),D(2),F(2),RAD(3)
	COMMON /WS/ K,JR,X,CL,E,ATP
	COMMON /REL/ IREL,ALDS,VBCON
	DATA PI,CC/3.14159265,18779.688/

CALCULATES PHASE SHIFTS BY SEPARATING SCHRODINGER EQUATION INTO TWO
FIRST ORDER EQUATIONS AND INTEGRATING USING RUNGE-KUTTA-MERSON ROUTINE
DA01A.
FOLLOWING THE METHOD OF HARTREE (1957) P81, THE TWO EQUATIONS ARE OF THE
FORM:

	D/D = RAD(2) + (L+1)RAD(1)/RAD(3)

	D/D = (2V-2E)RAD(1) - (L+1)RAD(2)/RAD(3)

WHERE

	RAD(1) = R X WAVE FUNCTION

	RAD(2) = DRAD(1)/DR - (L+1)RAD(1)/R

	RAD(3) = R

THE OPTIONAL L.S TERM IS ADDED TO THE SECOND EQUATION IN THE FORM:
	(L.S)/(C**2) X DV/DR X RAD(1)/RAD(3)	WHERE C IS 137.039 IN A.U.

THE ABOVE DIFFERENTIAL FUNCTIONS ARE CALCULATED IN ROUTINE DYBDX.

DEPENDING ON WHETHER THE L.S TERM IS INCLUDED THE STARTING CONDITIONS
ARE SET UP AS IN HARTREE (1957) P81, OR FROM THE CUTOFF L.S TERM AT
THE SECOND MESH POINTS.
INTEGRATION PROCEEDS OVER TEN VALUES OF THE BASIC GRID.  AFTER EACH RANGE
OF TEN VALUES, AND BEFORE PROCEEDING TO THE NEXT, THE STEP SIZE IS HALVED
AND THE ACCURACY OF THE CURRENT RANGE IS CHECKED.  IF THE REQUIRED
ACCURACY IS NOT REACHED AFTER HALVING THE STEP LENGTH FOUR TIMES, THE
ACHIEVED ACCURACY IS PRINTED OUT WITH A WARNING FLAG AND THE PROGRAM
PROCEEDS.  THIS OFTEN OCURS IN THE FIRST RANGE OF TEN AND CAN BE IGNORED
IF TEST IS SMALL.  SUBROUTINE DA01A ADVANCES THE INTEGRATION ACROSS ONE
GRID INTERVAL.
THE PHASE SHIFT IS THEN CALCULATED FROM THE LOGARITHMIC DERIVATIVE AT
THE MESH POINT IMMEDIATELY OUTSIDE THE MUFFIN-TIN RADIUS AND THE
LOGARITHMIC DERIVATIVE FOR THE UNSCATTERED WAVE OBTAINED USING THE BESSEL
FUNCTION ROUTINE CALCBF. THIS RESULT IS OBTAINED MODULO PI AND THE TRUE
PHASE SHIFT MAGNITUDE IS THEN CALCULATED BY DETERMINING THE NUMBER OF
ZEROES IN THE WAVEFUNCTION.

INPUT:	L		ANGULAR MOMENTUM
	Z		ATOMIC NUMBER
	NF		NUMBER OF HERMAN-SKILLMAN GRID POINTS
	RHO		MUFFIN-TIN RADIUS

OUTPUT:	STORE		(WAVEFUNCTION X RADIUS) OBTAINED FROM INTEGRATION
	DEL 		PHASE SHIFT

IF REQUIRED ACURACY NOT OBTAINED:
OUTPUT ON UNIT 6:
----------------

TEST		FORMAT(E13.5)		ACCURACY OF NONCONVERGED WAVEFTN

ROUTINES CALLED:  DA01A, CALCBF



SUBROUTINE CALCBF
-----------------

	SUBROUTINE CALCBF(Z,CL,BJ,BN)

INPUT:	Z		ARGUMENT
	CL		ORDER

OUTPUT:	BJ		SPHERICAL BESSEL 'J'
	BN		SPHERICAL BESSEL 'N'



SUBROUTINE DA01A
----------------

	SUBROUTINE DA01A(Y,E,A,B,C,D,NO,XO,DX)
	DIMENSION Y(1),E(1),A(1),B(1),C(1),D(1)

RUNGE-KUTTA-MERSON ROUTINE ADVANCES THE INTEGRATION OF A SET OF
DIFFERENTIAL EQUATIONS DY(I)/DX=F(I),I=1,N  BY ONE STEP FROM X0
TO X0+DX.

INPUT:	Y		VALUES OF FUNCTIONS AT POINT X0
	N		NUMBER OF FUNCTIONS TO BE CALCULATED
	X0		STARTING POINT FOR INTEGRATION
	DX		INTERVAL OVER WHICH INTEGRATION IS TO PROCEED

OUTPUT:	Y		VALUES OF FUNCTIONS AT POINT X0+DX
	E		TRUNCATION ERRORS
	A,B,C,D		CORRECTION FACTORS

ROUTINE CALLED:  DYBDX



SUBROUTINE DYBDX
----------------

	SUBROUTINE DYBDX(YT,F,N,XT)
	DIMENSION YT(2),F(2),ATP(200),X(200)
	COMMON /WS/ K,JR,X,CL,E,ATP
	COMMON /REL/ IREL,ALDS,VBCON

CALCULATES DIFFERENTIAL FUNCTIONS AS DESCRIBED ABOVE.
L.S TERM IS INCLUDED IF IREL=1, WITH CUTOFF BELOW SECOND GRID POINT.

INPUT:	YT		VALUE OF FUNCTIONS AT POINT XT
	N		NUMBER OF FUNCTIONS BEING CALCULATED
	XT		POINT AT WHICH DIFFERENTIAL FUNCTIONS CALCULATED

OUTPUT:	F		DIFFERENTIAL FUNCTIONS







EXAFS MATRIX ELEMENT SUBROUTINE: RONUC
--------------------------------------

	SUBROUTINE RONUC(Z,NM,VM,XM,WM,RM,E1,DE,NJE,LI,IR)
	DIMENSION XM(260),VM(260),Y(200),VL(200)
	DIMENSION STORE(200),PHS(15)
	DIMENSION PRE(100),PRM(100),CWM(260),CWL(200),DING(200)
	COMMON /WS/ KCF,JRCF,Y,CL,E,VL
	COMMON /WF/ WFC(260,15,2),WC(15),LC(15)

PURPOSE:

TO COMPUTE MATRIX ELEMENTS AND CROSS SECTOINS FOR ATOMIC EXAFS FOR AN
ARBITRARY EDGE USING THE CORE STATE WAVEFUNCTIONS INPUT IN ROUTINES
HSIN OR CLEMIN AND THE WAVEFUNCTIONS CALCULATED IN THE PHASE SHIFT
ROUTINE CALCFI.
NOTE THAT THE TWO OTHER INPUT OPTIONS ARE NOT CONSISTENT WITH CALLS TO
THIS ROUTINE AND THE SWITCH IREL FOR THE L.S TERM SHOULD BE ZERO. ALSO
THE INITIAL STATE ENERGY INPUT ON UNIT 5 IS GIVEN WITH RESPECT TO THE
VACUUM LEVEL AND IMMEDIATELY CORRECTED TO THE MUFFIN-TIN ZERO LEVEL.
AFTER PLOTTING OUT THE MUFFIN-TIN POTENTIAL AND CORE WAVE FUNCTION ON
THE HERMAN-SKILLMAN GRID (NOTE WAVEFUNCTION X RADIUS IS TABULATED), THE
ANGULAR PART OF THE MATRIX ELEMENT APPROPRIATE TO THE INITIAL CORE
STATE IS CALCULATED.  THEN FOR EACH FINAL STATE THE PHASE SHIFT
ROUTINE CALCFI IS CALLED TO PROVIDE THE FINAL STATE WAVEFUNCTION AND THE
MATRIX ELEMENTS ARE CALCULATED.  THE CROSS SECTIONS ARE THEN COMPUTED AND
PRINTED OUT.
NOTE: HOLE CORRECTIONS TO THE CHARGE DENSITY ARE HANDLED IN THE INPUT
STREAM BY USING WAVEFUNCTIONS FOR AN ATOM WITH ATOMIC NUMBER Z+1 BUT
REMOVING A CORE ELECTRON BY ADJUSTING THE OCCUPATION NUMBER OF THE 
APPROPRIATE CORE STATE.

INPUT:	Z		ATOMIC NUMBER
	NM		NUMBER OF POINTS ON INPUT GRID
	VM		INPUT FUNCTION: MUFFIN-TIN POTENTIAL IN HARTREES
					X RADIAL DISTANCE
	XM		INPUT GRID: LOUCKS' RADIAL GRID
	WM		MUFFIN-TIN ZERO LEVEL IN HARTREES
	RM		MUFFIN-TIN RADIUS
	E1		ENERGY (HARTREES) OF FINAL STATE FOR WHICH MATRIX
			ELEMENTS ARE TO BE CALCULATED
	DE		ENERGY INCREMENT
	NJE		NUMBER OF ENERGY VALUES
	LI		INITIAL STATE ANGULAR MOMENTUM
	IR		ATOM TYPE NUMBER

OUTPUT ON UNIT 6:
----------------

NCORE,ECORE,LI	FORMAT(I4,F11.5,I4)	INDEX FOR INITIAL STATE WAVE-
					FUNCTION, INITIAL STAE ENERGY
					IN HARTREES RELATIVE TO MUFFIN
					TIN ZERO LEVEL AND INITIAL STATE
					MOMENTUM
Y,VL,CWL	FORMAT(6E12.4)		HERMAN-SKILLMAN GRID, MUFFIN-TIN
					POTENTIAL X RADIAL DISTANCE AND
					INITIAL STATE WAVEFUNCTIONS
L,MF,LI,MI,	FORMAT(4I4,2F11.5)	ORDERS OF SPHERICAL HARMONICS OVER
XA,ZA					WHICH THE PRODUCT INTEGRALS XA
					AND ZA ARE SUMMED TO GIVE ANG
ANG		FORMAT(F11.5)		ANGULAR PART OF MATRIX ELEMENT
E,DEL,EMMM	FORMAT(2F9.4,E14.5)	FINAL STATE ENERGY (HARTREES),
					PHASE SHIFT FOR FINAL STATE,
					RADIAL PART OF MATRIX ELEMENT
PRE,PRM		FORMAT(F10.4,E14.5)	ENERGY (HARTREES) AT WHICH CROSS
					SECTION IS CALCULATED, EXAFS
					CROSS SECTION

ROUTINES CALLED:  CHGRID,CALCFI
FUNCTION CALLED:  BLM



FUNCTION BLM
------------

	FUNCTION BLM(IL1,IM1,IL2,IM2,IL3,IM3)
	DATA PI/3.14159265358979/

RETURNS THE INTEGRAL OF THE PRODUCT OF THREE SPHERICAL HARMONICS, EACH
OF WHICH CAN BE EXPRESSED AS A PRE-FACTOR MULTIPLIED BY A LEGENDRE
FUNCTION.  THE INTEGRAL OF THE THREE LEGENDRE FUNCTIONS FOLLOWS GAUNT'S
SUMMATION SCHEME SET OUT BY SLATER (1960) P309.

INPUT:	ILN,IMN,N=1,3	ORDERS OF THREE SPHERICAL HARMONICS

OUTPUT:	BLM		INTEGRAL OF PRODUCTS OF SPHERICAL HARMONICS





---------------------------------------------------------------------------
---------------------------------------------------------------------------

	COPIES OF MUFPOT OR OF THE INPUT PROGRAMS HERSKL AND HEX OR
	FURTHER INFORMATION ON THESE MAY BE OBTAINED BY WRITING TO:

		DR G. C. AERS
		ROOM B59
		SCIENCE RESEARCH COUNCIL
		DARESBURY LABORATORY
		DARESBURY
		WARRINGTON
		WA44AD
---------------------------------------------------------------------------
---------------------------------------------------------------------------