C PSTG2NX:: Parallel Version 1.2 of stg2nx 13/08/09 C C. P. Ballance and N. R. Badnell C C V1.1 Initial release C V1.2 Simplified MPI C C C N. R. BADNELL UoS v1.11 - QUB vx.x 28/11/08 C C PROGRAM NXRAD/STG2NX C C *** THIS IS THE RADIAL STAGE OF RMATRX NX C !*********************************************************************** module param implicit real*8(a-h,o-z) INCLUDE 'PARAM' end module param !*********************************************************************** module comm_interface implicit none include "mpif.h" public comm_init ! Initialize MPI public comm_barrier ! MPI barrier public comm_finalize ! Terminate MPI integer, public :: iam integer, public :: nproc SAVE private integer :: mpicom CONTAINS !----------------------------------------------------------------------------- subroutine comm_init() implicit none integer :: ier mpicom = MPI_COMM_WORLD call mpi_init(ier) call mpi_comm_rank(mpicom, iam, ier) call mpi_comm_size(mpicom, nproc, ier) return end subroutine comm_init !----------------------------------------------------------------------------- subroutine comm_barrier() implicit none integer :: ier call mpi_barrier(mpicom, ier) return end subroutine comm_barrier !----------------------------------------------------------------------------- subroutine comm_finalize() implicit none integer :: ier call mpi_finalize(ier) return end subroutine comm_finalize !----------------------------------------------------------------------------- end module comm_interface C C*********************************************************************** C PROGRAM MAIN use param use comm_interface, only : iam,nproc,comm_init,comm_barrier, A comm_finalize IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM3=IDIM9+1,KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9, * KDIM6=(IDIM4+1)/2,KDIM7=1, * KDIM8=1, * KDIM9=IDIM1+IDIM1-1,KDIM10=IDIM7*2, * KDIM11=1,KDIM12=IDIM1*IDIM12, * KDIM13=1) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) COMMON /BASDER/ FM,TLC,WR,IMEMM,ITST,JR,KM,MMM,NBTP1,NG COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON /BNDORB/ P(IDIM9,KDIM10) COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /COPY / ICOPY1,ICOPY2,ITOTAL,IPSEUD COMMON /CORE / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2), * MAXNC(IDIM2),MAXPN(IDIM2),ICHECK COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4 COMMON /FACT / GAMMA(IDIM14) COMMON /FUNVAL/ FRH(KDIM3),U(KDIM3),X COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI * ,BVALUE COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /POTVAL/ POVALU(KDIM10) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5) COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS COMMON /YKSTOR/ YK(IDIM7),RK(IDIM7),RK1(IDIM7),JN2,JL2,JN4,JL4,JKM C 1000 FORMAT(TR30,' END OF NXRAD'/TR30,' ------------'///) C C kickstart mpi C call comm_init() c write(0,*)'starting proc', iam call cpu_time(timei) C C READ IN AND WRITE OUT THE BASIC DATA. INITIALISE CERTAIN ARRAYS C AND VARIABLES NFACT=IDIM14 CALL SHRIEK(NFACT) C C REREAD STGANG DATA AND C CALCULATE RANGE OF CONTINUUM ANGULAR MOMENTUM LRANGE2 C FROM CONFIGURATION DATA C CALL RESETL(MJS,MINIM,ICONTN) C C FOR A CONTINUATION FROM RMATRX READ IN DATA, BOUND ORBITALS, C AND THE POTENTIAL FROM RMATRX FILE NX1 C IF (ICONTN .EQ. 1) THEN CALL STGRRX(MINIM,MJS) ELSE CALL STGRRD (MJS) C EVALUATE THE BOUND ORBITALS AT THE MESH POINTS. CALL EVALUE (MINIM) C GENERATE THE STATIC POTENTIAL CALL POTF END IF C C EVALUATE THE CONTINUUM ORBITALS CALL BASORB(MINIM) C C CALCULATE THE BUTTLE CORRECTION C CALL NEWBUT(MINIM,MJS) C C WRITE SERIAL FILE FOR STGHAM C CALL WRITAP C GENERATE THE RADIAL INTEGRALS CALL GENINT (MINIM) WRITE(IWRITE,1000) C call cpu_time(timef) time=timef-timei time = time/60.0 WRITE(IWRITE,999) TIME 999 FORMAT(//1X,'CPU TIME=',F9.3,' MIN') C call comm_finalize() C STOP END C C======================================================================= C SUBROUTINE BASFUN(NBT,LC,NODES,RA,BSTO,WINIT,DELTA,ETA) use param use comm_interface, only : iam,nproc C C NEW NUMERICAL R-MATRIX ORBITAL ROUTINE C*********************************************************************** C C OPERATING INSTRUCTIONS C C*********************************************************************** C C 1. THE USER MUST PROVIDE THE FOLLOWING INPUT DATA... C C NBT....... THE NUMBER OF FUNCTIONS TO WHICH THE SOLUTION IS TO BE C ORTHOGONALIZED...... IF NBT.GT.5 IS REQUIRED THE FIRST C INDEX IN THE ARRAYS US,P,U,DU,FR,FRH,FRM,ALAMDA,DELT, C SDELT,UNAME,ADEL,ADL,BDL SHOULD BE INCREASED. C C LC........ THE ANGULAR MOMENTUM VALUE C C RA........ THE BOUNDARY RADIUS C C BSTO...... THE VALUE OF THE LOGARITHMIC DERIVATIVE AT X=RA C C WINIT..... THE INITIAL ENERGY OR POTENTIAL MULTIPLE VALUE C .... IF ETA=0.0 OR IS LESS THAN 1.0E-8 THE WAVE C FUNCTION WILL BE EVALUATED AT THIS VALUE. IF ETA C IS GREATER THAN 1.0E-8 (TYPICALLY ETA=0.00001) THEN THE C PROGRAM WILL ITERATE FROM WR OR PR=WINIT TO THE VALUE C OF WR OR PR WHICH GIVES A SOLUTION SATISFYING THE C LOGARITHMIC BOUNDARY CONDITION AT X=RA C C DELTA..... THE INCREMENT IN THE ENERGY OR POTENTIAL MULTIPLE USED C FOR OBTAINING THE DERIVATIVE IN NEWTONS METHOD....THIS C SHOULD BE OF THE SAME ORDER AS ETA C C ETA....... THE PROGRAM WILL STOP ITERATING TO FIND THE EIGENVALUE C WHEN THE CORRECTION TO THE ENERGY OR POTENTIAL MULTIPLE C BECOMES LESS THAN ETA. C C IREAD..... THE INPUT PERIPHERAL NUMBER C C IWRITE.... THE OUTPUT PERIPHERAL NUMBER C C HINT...... THE BASIC INTEGRATION STEP LENGTH C C NIX....... THE NUMBER OF CHANGES OF INTEGRATION STEP OR THE NUMBER C OF SUBINTERVALS INTO WHICH THE INTERVAL X=0 TO RA IS C DIVIDED C C IHX(I),I=1,NIX.... THE MULTIPLE OF THE BASIC INTEGRATION STEP IN C EACH SUBINTERVAL C C IRX(I),I=1,NIX.... THE TOTAL NUMBER OF INTEGRATION STEPS TO THE C END OF THE ITH SUBINTERVAL C C C 2. THE USER MUST PROVIDE THE POTENTIAL FUNCTION AND ORTHOGONALISATION C FUNCTIONS AND STORE THEM IN THE ARRAYS POVALU(1200) AND P(6,1200) C C ** THESE MUST BE EVALUATED AS FOLLOWS --- MUST BEING EMPHASISED ** C C THE ODD ELEMENTS POVALU(2N-1) AND P(I,2N-1),N=1,IRX(NIX), SHOULD C CONTAIN THE FUNCTION VALUES AT THE HALF-MESH POINTS. C C THE EVEN ELEMENTS POVALU(2N) AND P(I,2N),N=1,IRX(NIX), SHOULD C CONTAIN THE VALUES AT THE MESH POINTS. C C*********************************************************************** C C DEBUGGING PARAMETERS C C*********************************************************************** C C IF THESE ARE SET EQUAL TO ZERO THERE IS NO INTERMEDIATE PRINT-OUT C C NBUG1.....IF THIS IS NON-ZERO THE INTERMEDIATE INTEGRATIONS AND C ENERGIES ARE OUTPUT C C NBUG3.....IF THIS IS NON-ZERO THE ARRAYS FOR THE DETERMINATION C OF THE MISMATCH ARE OUTPUT C C*********************************************************************** C C OUTPUT RESULTS C C*********************************************************************** C C ORB(K) , K=1,IRX(NIX)+1 CONTAINS THE FINAL SOLUTION C C NODES.... IS THE NUMBER OF NODES IN THE FUNCTION C C EIGEN IS THE EIGENVALUE C C ALAMDA(I),I=1,NBTP1 CONTAINS THE NORMALISED LAGRANGE MULTIPLIERS C C ANORMI IS THE STURMIAN NORMALISATION COEFFICIENT FOR R.GT.RA C C BVALUE IS THE LOGARITHMIC DERIVATIVE VALUE AT R=RA C C*********************************************************************** C C OTHER VARIABLE AND ARRAY DEFINITIONS C C*********************************************************************** C C SEE THE LONG WRITE-UP C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C C PARAMETER(LARGE=MZMAX,VSMALL=1.D1**(-MZMAX),QLARGE=1.D1**MZMX2) C PARAMETER(IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG) PARAMETER(KDIM3=IDIM9+1,KDIM10=IDIM7*2) C COMMON /BASDER/ FM,TLC,WR,IMEMM,ITST,JR,KM,MMM,NBTP1,NG COMMON /BNDORB/ P(IDIM9,KDIM10) COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /FUNVAL/ FRH(KDIM3),U(KDIM3),X COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI * ,BVALUE C ***** COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS COMMON /POTVAL/ POVALU(KDIM10) C DIMENSION AVAL(2),BNDRY(2),DBNDRY(2) DIMENSION ORTH(KDIM3),YIN(KDIM3),YOUT(KDIM3),DYIN(KDIM3), * DYOUT(KDIM3),US(KDIM3,IDIM7),DUS(KDIM3,IDIM7),DU(KDIM3), * FR(KDIM3),FRM(KDIM3),DELT(KDIM3,KDIM3) DIMENSION ADL(KDIM3+1),SDELT(KDIM3+1,KDIM3+1) DIMENSION ADEL(KDIM3+2,KDIM3+2),BDL(KDIM3+2) C CHARACTER*4 UNAME(4) C DATA UNAME/'U(1)','U(2)','U(3)','U(4)'/ C DATA UNAME/'U(1)','U(2)','U(3)','U(4)','U(5)','U(6)','U(7)'/ C 1001 FORMAT(8E15.6) 1002 FORMAT(/5H IEN=,I2,5X,5H JEN=,I2,5X,11H ENERGY WR=,F15.7,5H RYDS) 1003 FORMAT(54H **** DEBUGGING PRINT-OUT IN BASFUN FOR NBUG1=0 ****) 1004 FORMAT(/ 13H ERROR - IRX(,I1,27H)IS AN ODD NUMBER IN BASFUN/) 1005 FORMAT(/8X,4H R ,5(11X,A4)/) 1006 FORMAT(/20H OUTWARD INTEGRATION) 1007 FORMAT(/19H INWARD INTEGRATION) 1008 FORMAT(/19H INWARD INTEGRATION,I2) 1009 FORMAT(54H **** DEBUGGING PRINT-OUT IN BASFUN FOR NBUG3=0 ****) 1010 FORMAT(1H+,44X,2H /) 1011 FORMAT(/ 53H ERROR - NO CONVERGENCE IN BASFUN AFTER 99 ITERATIONS */) 1012 FORMAT(/ 55H ERROR - HINT,IHX,IRX AND RA ARE INCOMPATIBLE IN BASF *UN/) 1013 FORMAT(20X,7F15.7) 1014 FORMAT(/20H SOLUTION FROM MA01A/) 1015 FORMAT(/17H ARRAYS FOR MA01A/) 1016 FORMAT(/4H X1=,F15.8) C C THE LARGEST AND SMALLEST FLOATING POINT NUMBERS ON THIS C MACHINE, LARGE AND VSMALL, ARE SET IN A PARAMETER C STATEMENT. QLARGE IS THE SQUARE ROOT OF LARGE. C THEY ARE USED IN CONSTRUCTING A SCALING FACTOR TO PREVENT C UNDERFLOW AND OVERFLOW AT HIGH VALUES OF L. C C CHECK COMPATIBILITY OF HINT,IHX,IRX AND RA C S=ZERO J=0 DO 1 I=1,NIX H=DBLE(IHX(I))*HINT S=S+H*DBLE(IRX(I)-J) J=IRX(I) 1 CONTINUE IF(ABS(S-RA)-HALF*HINT) 3,3,2 2 WRITE(IWRITE,1012) RETURN C C CHECK THAT IRX(I) ARE EVEN INTEGERS C 3 DO 5 I=1,NIX J=IRX(I) M=MOD(J,2) IF(M) 4,5,4 4 WRITE(IWRITE,1004)I RETURN 5 CONTINUE C C EVALUATE AND INITIALIZE SOME COMMONLY USED PARAMETERS C ONEPT5=ONE+HALF IMATCH=IRX(NIX)-20 IF(ETA.GT.ZERO) IMATCH=IRX(NIX)-10 C C**** CHANGE TO ALLOW OUTWARD INTEGRATION TO REACH BOUNDARY C WILL HAPPEN AT HIGH L WHEN POINT OF INFLECTION IS OUTSIDE C RMATRIX BOUNDARY. VMB JAN 1986 C IF (ETA .EQ. 0.0 .AND. WINIT .GE. 0.0 .AND. NBT .EQ. 0) THEN IMATCH=IRX(NIX) END IF C IF(MOD(IMATCH,2).EQ.1) IMATCH=IMATCH+1 FM=ZERO ITST=0 NBTP1=NBT+1 NG=NBT I9=IRX(NIX)+1 DO 6 I=1,I9 ORB(I)=ZERO DORB(I)=ZERO 6 CONTINUE LANG=LC TLC=DBLE(LC*(LC+1)) LCM=LC+1 OVERLP=1.D0/DBLE(LCM) INST=1 HIMT=HINT*DBLE(IHX(INST)) C YD=HIMT**LCM C DYD=DBLE(LCM)*YD/HIMT IEN=1 JEN=1 WR=WINIT EIGEN=WR C C BEGINNING OF LOOP TO FIND THE ENERGY WR, OR POTENTIAL MULTIPLE C PR, WHICH GIVES A FUNCTION SATISFYING THE LOGARITHMIC BOUNDARY C CONDITION. THIS IS DONE BY USING NEWTONS METHOD TO FIND THE ZERO C OF THE FUNCTION BVAL-BSTO. THIS LOOP IS LEFT ONLY WHEN THE C ENERGY, OR POTENTIAL MULTIPLE, INCREMENT DVAL IS LESS THAN ETA. C C NOTE.... THIS LOOP IS ENTERED ONCE ONLY IF ETA=0.0 C DO 99 IEN=1,99 IF(IEN.LT.99) GO TO 7 WRITE(IWRITE,1011) RETURN 7 WR1=WR DO 97 JEN=1,2 8 IF(NBUG1.NE.0) THEN WRITE(IWRITE,1003) WRITE(IWRITE,1010) WRITE(IWRITE,1002)IEN,JEN,WR WRITE(IWRITE,1006) WRITE(IWRITE,1005)(UNAME(I),I=1,NBTP1) END IF C C C ****MODIFICATIONS FOR LARGE L C C CALCULATE TWOZ C TWOZ=0.0D0 IF (POVALU(2) .GT. TINY) THEN NIZ=NINT(POVALU(1)/POVALU(2)-HALF) IF (NIZ .GE. 1) THEN TWOZ=POVALU(1)*HIMT*HALF END IF END IF ZSC=TWOZ*HALF C C ESTIMATE THE POINT OF INFLECTION. IF IT IS BEYOND C THE RADIUS AT WHICH R**(L+1) OVERFLOWS CALCULATE C SCALING FACTOR C SCAST=1.0D0 IF (ETA .EQ. 0.0D0 .AND. WR .GT. 0.0D0 .AND. NBT .EQ. 0 1 .AND. LCM .GT. 3) THEN IF (ZSC .NE. 0.0D0) THEN WRSC=WR/ZSC**2 ELSE WRSC=WR END IF PTINF=(SQRT(1.D0+WRSC*TLC)-1.D0)/WRSC/ZSC RLIMIT=10**(LARGE*OVERLP) IF (PTINF .GE. RLIMIT) THEN SCAST1=RLIMIT/XR(I9) IF (SCAST1 .LT. 1.0) THEN SCAST=SCAST1**LCM C WRITE(IWRITE,'('' SCAST='',E14.5)') SCAST END IF END IF END IF C XTEST=XR(2) NSTART=1 IF (SCAST*XTEST**LCM .EQ. ZERO) THEN I2=-2 DO 205 INTX=1,NIX I1=I2+5 I2=IRX(INTX)-3 DO 200 I=I1,I2 IF (SCAST*XR(I)**LCM.GT. VSMALL)THEN NSTART=I-1 INST=INTX C WRITE(IWRITE,'('' NSTART='',I4,''INST='',I4)') NSTART,INST GO TO 201 END IF 200 CONTINUE 205 CONTINUE WRITE(IWRITE,'('' NO STARTING POSITION FOUND FOR L='',I4)')LANG STOP END IF 201 CONTINUE C C INITIALIZATION OF FUNCTION AT X=0.0 WHICH IS USED AS THE FIRST C POINT IN THE SIMPSONS RULE NORMALIZATION AND ORTHOGONALIZATION C INTEGRATIONS. C*****INITIALISATION AT X=XR(NSTART) TO OVERCOME UNDERFLOW C 9 X=XR(NSTART) DO 10 II=1,NBTP1 U(II) = ZERO DU(II) = ZERO 10 CONTINUE IF(JEN.EQ.2) GO TO 12 DO 11 I=1,NBTP1 C **** DO 202 III=1,NSTART US(I,III)=ZERO DUS(I,III)=ZERO 202 CONTINUE C 11 CONTINUE 12 IF(NBUG1.EQ.0) GO TO 13 WRITE(IWRITE,1001)X,(U(NV),NV=1,NBTP1) 13 CONTINUE C C EVALUATE THE FUNCTION AND DERIVATIVE AT HINT AND STORE THE C FUNCTION C***** EVALUATE AT X=XR(NSTART+1) C LP1=LCM LP2 = LC+2 C**** HIMT=HINT*IHX(INST) X=XR(NSTART+1) H=HIMT DO 19 K=1,NBTP1 IF(K.EQ.NBTP1) GO TO 17 C U(K) = P(K,2)*H*H/DBLE(4*LC+6) C***** CHANGE FOR START AT NSTART C***** CAUTION THIS HAS NOT BEEN TESTED BECAUSE EXPECT C***** NO UNDERFLOW PROBLEM FOR L IBUFF1+1 BUFF1(I)=ZERO 30 CONTINUE WRITE(JBUFF1,REC=NBUFF1) BUFF1 ENDIF C C WRITE NO EXCHANGE POINTERS TO RAD3 C WRITE(ITAPE3) ((NCCD(LDIF,L),LDIF=1,LLPDIM),L=1,LRANG2) RETURN END C======================================================================= SUBROUTINE GENMBB(NBBPOL,IBBPOL) use param use comm_interface, only : iam,nproc C C GENERATES ALL THE BOUND-BOUND MULTIPOLE INTEGRALS C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM4=MZLMX, * IDIM7=MZPTS,IDIM8=MZAMP, * IDIM11=MZSLT,IDIM12=MZNR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM6=(IDIM4+1)/2,KDIM9=IDIM1+IDIM1-1, * KDIM12=IDIM1*IDIM12) PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9) C COMMON /CORE / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2), * MAXNC(IDIM2),MAXPN(IDIM2),ICHECK COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1, * JBUFF1,JBUFF2 DIMENSION X(2) DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6) C 1000 FORMAT (' LAM L1 L1P N1 N1P RVAL') 1001 FORMAT (2X, 5I5, E13.5) C IF (NBUG6 .GT. 0) THEN WRITE (IWRITE,1000) END IF C C C ZEROIZE ARRAYS C DO 3,I=1,LRANG1 DO 2,J=1,LRANG1 DO 1,K=1,LAMIND IBBPOL(I,J,K)=-999 1 CONTINUE 2 CONTINUE 3 CONTINUE NBBPOL=0 LUP=LRANG1-1 DO 10,L1=0,LUP MAXC1=MAXNC(L1+1) MAXHF1=MAXNHF(L1+1) IF (MAXC1.EQ.MAXHF1) GOTO 10 DO 20,L2=L1,LUP MAXC2=MAXNC(L2+1) MAXHF2=MAXNHF(L2+1) IF (MAXC2.EQ.MAXHF2) GOTO 20 LAMLO=ABS(L1-L2) IF (LAMLO.EQ.0) THEN LAMLO=2 ENDIF LAMST=(LAMLO+1)/2-1 LAMUP=MIN(L1+L2,LAMAX) DO 30,LAM=LAMLO,LAMUP,2 LAMST=LAMST+1 IBBPOL(L1+1,L2+1,LAMST)=IBUFF1+1+(LBUFF1*NBUFF1) IF (IBUFF1+1.GT.LBUFF1) THEN PRINT*,' ERROR IN CODE STOP IN GENMBB' STOP ENDIF DO 40,N1=MAXC1+1,MAXHF1 IF (L1.EQ.L2) THEN N2UP=N1 ELSE N2UP=MAXHF2 ENDIF DO 50,N2=MAXC2+1,N2UP CALL RADINT(N1,L1+1,N2,L2+1,LAM,X(1)) IUP=1 DO 60,I=1,IUP NBBPOL=NBBPOL+1 IBUFF1=IBUFF1+1 BUFF1(IBUFF1)=X(I) IF (IBUFF1.EQ.LBUFF1) THEN WRITE(JBUFF1,REC=NBUFF1) BUFF1 NBUFF1=NBUFF1+1 IBUFF1=0 ENDIF C WRITE (IWRITE,1001) LAM,L1,L2,N1,N2,X(I) IF (NBUG6 .GT. 0) THEN WRITE (IWRITE,1001) LAM,L1,L2,N1,N2,X(I) END IF 60 CONTINUE 50 CONTINUE 40 CONTINUE 30 CONTINUE 20 CONTINUE 10 CONTINUE RETURN END C C********************************************************** SUBROUTINE LSQ(P,Q,A,IMAX) use param use comm_interface, only : iam,nproc IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C C DIMENSION P(6),Q(6),X(3,3),C(3),A(3) C C 600 FORMAT(/18H COEFFICIENTS ARE ,3F14.7) 601 FORMAT(' X(I,J) = ',F14.7) 602 FORMAT(/'C(J) = ',F14.7) C C N=IMAX C DO 100 I=1,3 DO 101 J=1,3 X(I,J)=0.0D0 L=I+J-2 IF(L.NE.0) GOTO 105 X(I,J)=N GOTO 101 105 DO 102 K=1,N X(I,J)=X(I,J)+P(K)**L 102 CONTINUE C WRITE(6,601) X(I,J) 101 CONTINUE 100 CONTINUE C DO 103 J=1,3 C(J)=0.0D0 L=J-1 IF(L.NE.0) GOTO 106 DO 107 K=1,N C(J)=C(J)+Q(K) 107 CONTINUE GOTO 103 106 DO 104 K=1,N C(J)=C(J)+(P(K)**L)*Q(K) 104 CONTINUE C WRITE(6,602) C(J) 103 CONTINUE C CALL MA01A(X,B,3,0,1,3,1) CALL PMUL(X,C,A,1,3) C C C RETURN END C*********************************************************************** SUBROUTINE MA01A(A,B,M,N,M1,IA,IB) use param use comm_interface, only : iam,nproc C C SOLUTION OF SIMULTANEOUS EQUATIONS AND OR MATRIX INVERSION C C C A THE M*M MATRIX OF LEFT HAND SIDES OR THE MATRIX BEING C INVERTED. OVERWRITTEN ON EXIT BY THE INVERSE MATRIX C C B THE M*N MATRIX OF THE RIGHT HAND SIDES. OVERWRITTEN C ON EXIT BY SOLUTIONS C C M THE ORDER OF THE A-MATRIX. THIS MUST NOT BE GREATER C THAN THE DIMENSION OF THE C AND IND ARRAYS. THE UPPER C LIMIT CAN BE EXTENDED BY RECOMPILING WITH LARGER C DIMENSIONS FOR THE PRIVATE ARRAYS C AND IND C C N THE NUMBER OF THE RIGHT HAND SIDES IN THE C SIMULTANEOUS EQUATIONS C C M1 =0 ONLY SIMULTANEOUS EQUATIONS ARE SOLVED IF N.GT.0 C IF N=0 A FURTHER ENTRY TO MA01A WITH M1.LT.0 C REQUIRED TO OBTAIN THE INVERSE OF A C .GT.0 MATRIX INVERSION. IN ADDITION SIMULTANEOUS C EQUATIONS ARE SOLVED IF N.GT.0 C .LT.0 ONLY USED IF PREVIOUS ENTRY TO MA01A C WITH M1=0. IN THIS CASE THE MATRIX INVERSION IS C COMPLETED C C IA DEFINES THE DIMENSIONS OF THE ARRAY WHERE C THE A-MATRIX IS STORED C C IB DEFINES THE SECOND DIMENSION OF THE ARRAY WHERE C B-MATRIX IS STORED C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM9=MZLAG) PARAMETER(KDIM3=IDIM9+1) C DIMENSION A(IA,IA),B(IA,IB),C(KDIM3+2),IND(KDIM3+2) C ZERO=DBLE(0) ONE=DBLE(1) MM=M-1 IF(M1) 32,1,1 C C CHECK FOR ZERO DIAGONAL ELEMENTS C 1 IF(MM) 55,2,2 2 DO 3 I=1,M IF(A(I,I).EQ.ZERO) A(I,I)=DBLE(10)**(-35) 3 CONTINUE C C DO THE TRIVIAL CASE OF A 1*1 MATRIX C IF(MM) 55,4,7 4 IF(N-1) 5,6,7 5 A(1,1)=ONE/A(1,1) GO TO 55 6 B(1,1)=B(1,1)/A(1,1) GO TO 5 C C THIS IS NOT A TRIVIAL CASE C C FIND THE FIRST PIVOTAL ELEMENT AND STORE THE CORRESPONDING ROW C NUMBER IN I4. IND DEFINES THE ORDER OF THE ROWS OF THE ORIGINAL C A-MATRIX BEFORE ROW INTERCHANGE C 7 AMAX=ZERO DO 9 I=1,M IND(I)=I IF(ABS(A(I,1))-AMAX) 9,9,8 8 AMAX=ABS(A(I,1)) I4=I 9 CONTINUE C C EACH TIME THROUGH THE FOLLOWING LOOP THE A-MATRIX IS C REDUCED BY ONE C IF(MM) 24,24,10 10 DO 23 J=1,MM C C INTERCHANGE THE I4TH AND THE JTH ROWS OF THE A-MATRIX AND STORE C ORDER IN IND IF I4 .NE.J C IF(I4-J) 15,15,11 11 ISTO=IND(J) IND(J)=IND(I4) IND(I4)=ISTO DO 12 K=1,M STO=A(I4,K) A(I4,K)=A(J,K) A(J,K)=STO 12 CONTINUE C C INTERCHANGE THE I4TH AND THE JTH ROWS OF THE B-MATRIX IF N.GT. 0 C IF(N) 15,15,13 13 DO 14 K=1,N STO=B(I4,K) B(I4,K)=B(J,K) B(J,K)=STO 14 CONTINUE C C THE JTH ROW NOW CONTAINS THE PIVOTAL ELEMENT IN THE JTH POSITION C ELIMINATE THE JTH ELEMENT FROM EACH OF THE REMAINING ROWS OF THE C A-MATRIX AND THE B-MATRIX AND STORE THE MULTIPLIERS IN THE LOWER C TRIANGLE C 15 AMAX=ZERO J1=J+1 DO 22 I=J1,M A(I,J)=A(I,J)/A(J,J) DO 18 K=J1,M A(I,K)=A(I,K)-A(I,J)*A(J,K) IF(K-J1) 16,16,18 C C FIND THE NEXT PIVOTAL ELEMENT AND STORE THE CORRESPONDING ROW C NUMBER IN I4 C 16 IF(ABS(A(I,K))-AMAX) 18,18,17 17 AMAX=ABS(A(I,K)) I4=I 18 CONTINUE 19 IF(N) 22,22,20 20 DO 21 K=1,N B(I,K)=B(I,K)-A(I,J)*B(J,K) 21 CONTINUE 22 CONTINUE 23 CONTINUE C C THE ELIMINATION IS NOW COMPLETE AND THE A-MATRIX HAS BEEN C REDUCED TO THE PRODUCT OF AN UPPER AND LOWER TRIANGLE MATRIX C 24 IF(N) 31,31,25 C C NOW CARRY OUT THE BACK SUBSTITUTION AND STORE RESULT IN THE C B-MATRIX IF THERE IS AT LEAST ONE RIGHT HAND SIDE C 25 DO 30 I1=1,M I=M+1-I1 DO 29 J=1,N IF(M-I) 28,28,26 26 I2=I+1 DO 27 K=I2,M B(I,J)=B(I,J)-A(I,K)*B(K,J) 27 CONTINUE 28 B(I,J)=B(I,J)/A(I,I) 29 CONTINUE 30 CONTINUE 31 IF(M1) 55,55,32 C C REPLACE THE A-MATRIX BY ITS INVERSE WHEN M1.NE. ZERO C C FIRST INVERT THE LOWER TRIANGLE MATRIX AND STORE ON ITSELF C 32 IF(MM) 40,40,33 33 DO 39 I1=1,MM I=M+1-I1 I2=I-1 DO 37 J1=1,I2 J=I2+1-J1 J2=J+1 W1=-A(I,J) IF(I2-J2) 36,34,34 34 DO 35 K=J2,I2 W1=W1-A(K,J)*C(K) 35 CONTINUE 36 C(J)=W1 37 CONTINUE DO 38 K=1,I2 A(I,K)=C(K) 38 CONTINUE 39 CONTINUE C C NOW INVERT THE UPPER TRIANGLE MATRIX AND FORM THE ORIGINAL C A-MATRIX APART FROM COLUMN INTERCHANGE. THIS OVERWRITES THE C ORIGINAL A-MATRIX C 40 DO 50 I1=1,M I=M+1-I1 I2=I+1 W=ONE/A(I,I) DO 48 J=1,M IF(I-J) 41,42,43 41 W1=ZERO GO TO 44 42 W1=ONE GO TO 44 43 W1=A(I,J) 44 IF(I1-1) 47,47,45 45 DO 46 K=I2,M W1=W1-A(I,K)*A(K,J) 46 CONTINUE 47 C(J)=W1 48 CONTINUE DO 49 J=1,M A(I,J)=C(J)*W 49 CONTINUE 50 CONTINUE C C RE-ORDER THE COLUMNS OF THE INVERSE A-MATRIX TO COINCIDE WITH C THE ORDER OF THE ROWS OF THE A-MATRIX ON INPUT C DO 54 I=1,M 51 IF(IND(I)-I) 52,54,52 52 J=IND(I) DO 53 K=1,M STO=A(K,I) A(K,I)=A(K,J) A(K,J)=STO 53 CONTINUE ISTO=IND(J) IND(J)=J IND(I)=ISTO GO TO 51 54 CONTINUE C 55 RETURN END C======================================================================= SUBROUTINE MESH use param use comm_interface, only : iam,nproc C C NEW VERSION WRITTEN BY VAL BURKE 7 OCTOBER 1985 C C THIS ROUTINE AUTOMATICALLY GENERATES THE INTEGRATION MESH, C BASED ON THE NUCLEAR BEHAVIOUR OF BOUND ORBITALS, THE NUMBER C OF CONTINUUM ORBITAL LOOPS AND THE CURRENT ARRAY DIMENSIONS C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C C PARAMETER(MSHDIM=16) C PARAMETER(IDIM1=MZLR1, * IDIM6=MZNIX,IDIM7=MZPTS, * IDIM11=MZSLT,IDIM12=MZNR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM12=IDIM1*IDIM12) C COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM C DATA MINFAC/1/,MAXFAC/2/ C NCORSE=NRANG2*MSHDIM IF (NCORSE .LT. 96) NCORSE=96 C C CALCULATE THE COARSEST MESH, HMAX, AND THE MESH REQUIRED C NEAR THE ORIGIN. CALCULATE HMIN = HMAX/2**M WHERE M+1 IS C THE NUMBER OF STEP SIZES C HMAX=RA/NCORSE HINNER=.025D0/NZ DELTA=HINNER/5.0D0 M=0 HMIN=HMAX 1 IF (HMIN .GT. HINNER+DELTA) THEN M=M+1 HMIN=HMIN/2 GO TO 1 END IF C HINT=HMIN NIX=M+1 IF (NIX .GT. IDIM6) CALL RECOV1(6,NIX,IDIM6) C C SET UP THE IHX ARRAY C IH=1 DO 2 I=1,NIX IHX(I)=IH IH=IH+IH 2 CONTINUE C C CONSIDER SEPARATELY M .LT. 4 AND M .GE. 4 C NA IS THE NUMBER OF STEPS AT EACH STEP SIZE C IT IS A MULTIPLE OF 16 AND CAN TAKE VALUES C FROM 16*MAXFAC DOWN TO 16*MINFAC C NAFAC=MAXFAC C IF (M .GE. 4) THEN 3 NA=16*NAFAC NTOT=NCORSE+(M-1)*NA+NA/8 IF (NTOT .GE. IDIM7) THEN NAFAC=NAFAC-1 IF (NAFAC .LT. MINFAC) CALL RECOV1(7,NTOT,IDIM7) GO TO 3 END IF C C SET UP IRX ARRAY C IRX(1)=NA IRX(2)=NA+IRX(1) IRX(3)=NA+IRX(2) IRX(4)=NA+NA/8+IRX(3) DO 4 I=5,M IRX(I)=NA+IRX(I-1) 4 CONTINUE C ELSE C MPOW2=2**M 5 NA=16*NAFAC NTOT=NCORSE+M*NA-NA*(MPOW2-1)/MPOW2 IF (NTOT .GE. IDIM7) THEN NAFAC=NAFAC-1 IF (NAFAC .LT. MINFAC) CALL RECOV1(7,NTOT,IDIM7) GO TO 5 END IF C C FILL IRX ARRAY C IA=0 DO 6 I=1,M IA=IA+NA IRX(I)=IA 6 CONTINUE C END IF C C NUMBER OF STEPS AT EACH STEP SIZE MUST BE EVEN C IF (MOD(NTOT,2) .NE. 0) THEN NTOT=NTOT+1 IRX(M)=IRX(M)+2 END IF C IRX(NIX)=NTOT C C PERFORM CHECK C RVAL=ZERO IR=0 DO 7 I=1,M+1 RVAL=RVAL+HINT*(IRX(I)-IR)*IHX(I) IR=IRX(I) 7 CONTINUE CW WRITE (IWRITE,'('' RVAL='',E14.7)')RVAL C RETURN END C********************************************************************* SUBROUTINE NEWBUT(MINIM,MJS) use param use comm_interface, only : iam,nproc IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG, * IDIM11=MZSLT,IDIM12=MZNR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM3=IDIM9+1,KDIM10=IDIM7*2, * KDIM12=IDIM1*IDIM12) C LOGICAL OK C COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON/BASIC/ BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, 1 LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /BNDORB/ P(IDIM9,KDIM10) COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN COMMON/INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /ORBOUT/ ORB(IDIM7),DORB(IDIM7),EIGEN,ALAMDA(KDIM3),ANORMI * ,BVALUE COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH C DIMENSION EN(IDIM3,IDIM2),AMP(IDIM3,IDIM2) DIMENSION EL(6),RT(6),RCN(6),A(3),EK(6) C 1000 FORMAT(' BUTTLE CORRECTION COEFFICIENTS FOR EACH CHANNEL L'/ 1 ' L',12X,'A(0)',16X,'A(1)',16X,'A(2)'/) 1001 FORMAT(' ',I4,3F20.6) C 6000 FORMAT(' TWO BUTTLE POLES ARE EQUAL-CONTINUE WITHOUT CORRECTION') 6001 FORMAT(' SUBROUTINE NEWBUT - ENERGY RANGE TOO SMALL TO AVOID BUTTL *E POLE, L = ',I3) 6003 FORMAT(' '/28X,' SUBROUTINE NEWBUT'/28X,' -----------------'/ 1 ' EK2MAX = ',F14.7,' EK2MIN = ',F14.7) 6004 FORMAT(' SUBROUTINE NEWBUT - CANNOT EVALUATE BUTTLE CORRECTION FOR * L = ',I3,' AND ENERGY ',F14.7,'RYD.') C C LOOP OVER CHANNEL ANGULAR MOMENTA TO CALCULATE BUTTLE POLES C UP TO LRANG2 C C IPTS=2*IRX(NIX) I9=IRX(NIX)+1 C WRITE(IWRITE,6003) EK2MAX,EK2MIN WRITE(IWRITE,1000) C C SET COEFF(3,1) TO INDICATE TO THE ASYMPTOTIC PROGRAM C WHICH FITTING PROCEDURE IS USED C IF (MINIM .GT. 1) THEN IF (MJS .EQ. 1) THEN COEFF(3,1)=-20000.0D0 ELSE COEFF(3,1)=0.0D0 END IF END IF C DO 50 L=MINIM,LRANG2 NPOLES=NVARY(L) C C SET UP ENERGY POINTS IN EK() C C WRITE(IWRITE,6005) IMAX=6 EKINT=IMAX-1 EKINT=(EK2MAX-EK2MIN)/EKINT E=EK2MIN-EKINT DO 99 I=1,IMAX E=E+EKINT EK(I)=E 99 CONTINUE C L1=L-1 NBT=0 IF(L.LE.LRANG1) NBT=MAXNLG(L)-L+1 IF (NBT .GT. 0) READ(60,REC=L) P ITEST=0 IF(L.GT.LRANG1) GOTO 51 ITEST=MAXNHF(L)-MAXNLG(L) 51 IF (ITEST .GT. 0) GO TO 52 C C IF MAXNHF(L).EQ.MAXNLG(L) THE BUTTLE POLES AND AMPLITUDES ARE EQUAL C TO THE CONTINUUM EIGENVALUES AND FUNCTIONS EVALUATED IN STG1 AND C STORED IN EIGENS AND ENDS. C C DO 53 N=1,NPOLES AMP(N,L)=ENDS(N,L) EN(N,L)=EIGENS(N,L) 53 CONTINUE GO TO 60 C C IF MAXNLG(L).LT.MAXNHF(L) EVALUATE BUTTLE POLES C 52 ITEST=0 OK=.TRUE. DO 54 N=1,NPOLES NC=N-ITEST IF(NC) 55,55,56 55 W=EN(N,L) EH=0.0D0 GOTO 57 56 W=EIGENS(NC,L) EH=0.0D0 57 CALL BASFUN(NBT,L1,NODE,RA,BSTO,W,DELTA,EH) EN(N,L)=EIGEN AMP(N,L)=ORB(I9) N1=N-1 IF(N1) 54,54,58 58 DO 59 NLN=1,N1 IF(ABS(EIGEN-EN(NLN,L)).LT.TINY) OK=.FALSE. 59 CONTINUE 54 CONTINUE IF (OK) GO TO 60 WRITE(IWRITE,6000) RETURN C C CHECK THAT THE CHANNEL ENERGIES DO NOT COINCIDE WITH THE C BUTTLE POLES AND THEN CALCULATE COEFFICIENTS FOR QUADRATIC FIT. C C C SET CHANNEL ENERGIES AND TEST AGAINST POLES C C C FROM OPACITY VERSION OF STGR C C ENSURE THAT EL(I) NOT TOO CLOSE TO POLES. EXCLUDE RANGE C ((1-X)*E(J)+X*E(J-1)) TO ((1-X)*E(J)+X*E(J+1)) C WHERE E(J)=EN(J,L) AND X=0.2 C 60 DO 63 I=1,IMAX EL(I)=EK(I) 63 CONTINUE X=0.2D0 NA=1 DO 35 I=1,IMAX C C FIND FIRST POLE ABOVE EL(I) C DO 20 N=NA,NPOLES IF (EN(N,L) .GT. EL(I)) THEN NN=N GO TO 30 END IF 20 CONTINUE 30 NA=NN EA=EN(NA,L) C C CASE OF ALL POLES ABOVE C IF (NA .EQ. 1) THEN EP=EN(2,L) D=X*(EP-EA) IF ((EA-EL(I)) .LT. D) EL(I)=EA-D C C SOME POLES BELOW C ELSE C C NEAREST POLE BELOW C EB=EN(NA-1,L) D=X*(EA-EB) IF ((EL(I)-EB) .LT. D) EL(I)=EB+D IF ((EA-EL(I)) .LT. D) EL(I)=EA-D END IF 35 CONTINUE C C CHECK NUMBER OF INDEPENDENT VALUES C INDEP=1 DO 36 N=2,IMAX IF (EL(N-1) .NE. EL(N)) INDEP=INDEP+1 36 CONTINUE IF (INDEP .LT. 3) THEN WRITE (IWRITE,6001) L1 STOP END IF C C C ARRAY EL() NOW CONTAINS ACCEPTABLE ENERGIES FOR ANGULAR MOMENTUM C L . NOW EVALUATE THE BUTTLE CORRECTION AT THESE ENERGIES. C 67 DO 68 II=1,IMAX IF(ABS(EL(II)).LT.TINY)EL(II)=EL(II)+TWO*TINY EH=ZERO CALL BASFUN(NBT,L1,NODE,RA,BSTO,EL(II),DELTA,EH) RT(II)=BVALUE-BSTO IF(ABS(RT(II)).GT.TINY) GO TO 692 WRITE(IWRITE,6004) L1,EL(II) C IFLAG=1 STOP 692 SUM=ZERO DO 69 N=1,NPOLES IF(AMP(N,L).EQ.0.0D0) GO TO 69 SUM=SUM+AMP(N,L)*AMP(N,L)/(EN(N,L)-EL(II)) 69 CONTINUE SUM=SUM/RA R=ONE/RT(II) RCN(II)=R-SUM 68 CONTINUE C C FIT BUTTLE CORRECTION USING MJS METHOD IF MJS IS NON ZERO C UNLESS BSTO IS ZERO C USE FIT TO QUADRATIC IF MJS =0 OR BSTO IS ZERO C IF (BSTO .NE. ZERO .OR. MJS .EQ. 0) THEN C CALL LSQ(EL,RCN,A,IMAX) DO 691 I=1,3 COEFF(I,L)=A(I) 691 CONTINUE C DELTB=0.0D0 DO 693 I=1,IMAX DD=ABS(RCN(I)-A(1)-EL(I)*(A(2)+EL(I)*A(3)))/RCN(I) 693 IF(DELTB.LT.DD)DELTB=DD C ELSE C C FOR BSTO.NE.0. USE IMPROVED FITTING PROCEDURE CALL BUTFIT(IMAX,EL,RCN,RA,EIGENS(NPOLES,L),ALPHA,BETA, + NBUT,DELTB) COEFF(1,L)=ALPHA COEFF(2,L)=BETA COEFF(3,L)=-10000*NBUT C ENDIF C WRITE(IWRITE,1001) L1,(COEFF(I,L),I=1,3) C C ACCURACY OF FIT WRITE(IWRITE,6010)DELTB 6010 FORMAT(25X,'LARGEST ERROR IN FIT =',E15.6) C C 50 CONTINUE RETURN END C======================================================================= SUBROUTINE ONEELE(N11,L1,N12,L2,ALBVAL) use param use comm_interface, only : iam,nproc C C CALCULATE THE ONE-ELECTRON MATRIX ELEMENT BETWEEN ORBITALS C DEFINED BY THE QUANTUM NUMBERS (N1,L1-1) AND (N2,L2-1), C AND STORE THE RESULT IN ALBVAL. C WHEN BOTH ORBITALS ARE CONTINUUM ORBITALS THE SCHMIDT C COEFFICIENTS ARE USED TO EXPRESS ONE CONTINUUM ORBITAL IN TERMS C OF ORBITALS SATISFYING A DIFFERENTIAL EQUATION AND BOUND ORBITALS C WHEN ONE OR BOTH ORBITALS ARE BOUND THEN C ANALYTIC DIFFERENTIATION OF THE BOUND ORBITAL IS USED. C SEE COMMENTS IN EVALUE AND CORECT FOR DETAILS OF THE C CORRECTION APPLIED TO THE BOUND ORBITALS. C NRB: C ADDED MASS-VELOCITY CONTRIBUTION IN THE HIGH-L LIMIT I.E. C NO BOUND-BOUND OR BOUND-CONTINUUM, JUST CONTINUUM-CONTINUUM C WITH NO LAGRANGE OR SCHMIDT CONTRIBUTION. VALID FOR A CONTINUATION C RUN WHERE (NORMALLY) THE MIN CONTINUUM L EXCEEDS THE MAX BOUND L. C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9,KDIM10=IDIM7*2, * KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) C COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON /COPY / ICOPY1,ICOPY2,ITOTAL,IPSEUD COMMON /CORE / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2), * MAXNC(IDIM2),MAXPN(IDIM2),ICHECK COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /POTVAL/ POVALU(KDIM10) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5) COMMON /NRBREL/IRELMV C C RETURN IF THE ANGULAR MOMENTA ARE NOT COMPATIBLE C N1=N11 N2=N12 ALBVAL = ZERO IF(L1.NE.L2) RETURN C C DEFINE THE POSITION OF THE MODEL POTENTIAL IN THE POTHAM ARRAY. C LLL=LPOSX(L1) 1 ZN =DBLE(NZ+NZ) ALSQ=DBLE(L1*(L1-1)) C C FIND IF ONE ORBITAL IS BOUND. IF SO PUT IT INTO THE SECOND C POSITION SO THAT THE ORBITAL DEFINED BY (N2,L2) IS BOUND. C MAXHF = L1-1 IF(L1.LE.LRANG1) MAXHF = MAXNHF(L1) IF(N2.LE.MAXHF) GO TO 2 IF(N1.GT.MAXHF) GO TO 11 N3=N1 N1=N2 N2=N3 2 I1=IPOS(N1,L1) I2=IPOS(N2,L2) IF (NUMERC.EQ.0) THEN M=0 DO 40,I=1,L2-1 DO 50,K=I,MAXNHF(I) M=M+1 50 CONTINUE 40 CONTINUE M2=NCO(M+N2-L2+1) CALL CORECT(N2,L2,M2,R0,SIGMA,C1) ENDIF IST = 1 R = ZERO C C CARRY OUT THE INTEGRATION USING SIMPSONS RULE C DO 10 I=1,NIX H =HINT*DBLE(IHX(I)) MI=1 SUM = ZERO IFI = IRX(I) DO 9 JP=IST,IFI R = R+H C C DETERMINE THE CONTRIBUTION INVOLVING THE ANGULAR MOMENTUM C AND THE NUCLEAR CHARGE C RINV = ONE/R J1=JP+1 IF(IPSEUD.GT.0) POVAL=POTHAM(LLL,JP) IF(IPSEUD.LT.0)POVAL=TWO*POTHAM(LLL,JP)*RINV IF(IPSEUD.EQ.0) POVAL=ZN*RINV X1=UJ(J1,I2)*(POVAL-ALSQ*RINV*RINV) C C DETERMINE THE CONTRIBUTION INVOLVING THE SECOND DERIVATIVE C ACTING ON THE BOUND-ORBITAL C IF (NUMERC.EQ.0) THEN C C ANALYTIC BOUND ORBITALS ARE BEING USED C PV=ZERO ALPHA=(R-R0)/SIGMA IF(-ALPHA.LT.TWELVE)PV=-(TWO*ALPHA*ALPHA-ONE)*(TWO/(SIGMA*SIGMA * ))*C1*EXP(-ALPHA*ALPHA) DO 3 JK=1,M2 IR = IRAD(L2,N2,JK) ZEE=-ZE(L2,N2,JK) ATERM=DBLE(IR*(IR-1))+TWO*ZEE*DBLE(IR)*R+ZEE*ZEE*R*R PV=PV+C(L2,N2,JK)*EXP(ZEE*R)*ATERM*R**(IR-2) 3 CONTINUE ELSE C C NUMERICAL BOUND ORBITALS ARE BEING USED C PV=UDP(I2,J1) ENDIF C C DETERMINE THE SIMPSONS RULE WEIGHTING FACTOR C GO TO (4,5),MI 4 MI = 2 WEIGHT = FOUR GO TO 8 5 MI=1 IF(JP.EQ.IFI) GO TO 6 WEIGHT = TWO GO TO 8 6 IF(I.EQ.NIX) GO TO 7 WEIGHT = THREE GO TO 8 7 WEIGHT = ONE C C ADD IN THE CONTRIBUTIONS C 8 SUM=SUM+UJ(J1,I1)*(X1+PV)*WEIGHT 9 CONTINUE IST = IFI+1 ALBVAL = ALBVAL-SUM*H/SIX 10 CONTINUE RETURN C C BOTH ORBITALS ARE CONTINUUM ORBITALS. THE SCHIMDT COEFFICIENTS C ARE USED TO EXPRESS ONE CONTINUUM ORBITAL IN TERMS OF ORBITALS C SATISFYING A DIFFERENTIAL EQUATION AND BOUND ORBITALS. THE C SECOND DIFFERENTIATION CAN THEN BE CARRIED OUT ANALYTICALLY. C 11 CONTINUE I1=IPOS(N1,L1) I2=IPOS(N2,L2) C C ITEST DETERMINES WHETHER SCHMIDT ORTHOGONALIZATION IS USED C MAXLG = MAXHF IF(L1.LE.LRANG1) MAXLG = MAXNLG(L1) ITEST = MAXHF-MAXLG C C FIRST EVALUATE CONTRIBUTION FROM THE CONTINUUM ORBITALS USING C SIMPSONS RULE C 12 SUM = ZERO X = ZERO IST = 1 DO 17 I=1,NIX H=HINT*DBLE(IHX(I)) MI=1 ADD = ZERO IFI = IRX(I) DO 16 J=IST,IFI JP=J+1 X = X+H C WEIGHT GO TO (14,15),MI 14 MI = 2 WEIGHT=FOUR GO TO 13 15 MI=1 IF(J.EQ.IFI)GO TO 69 WEIGHT=TWO GO TO 13 69 IF(I.EQ.NIX)GO TO 68 WEIGHT=THREE GO TO 13 68 WEIGHT=ONE C ADD 13 IF(IPSEUD.GT.0) POVAL2=POTHAM(LLL,J) IF(IPSEUD.LT.0)POVAL2=TWO*POTHAM(LLL,J)/X IF(IPSEUD.EQ.0) POVAL2=ZN/X POVAL1=POVALU(2*J) ADD=ADD+UJ(JP,I1)*UJ(JP,I2)*WEIGHT*(POVAL2-POVAL1) 16 CONTINUE IST = IFI+1 ADD = ADD*H/THREE SUM=SUM+ADD 17 CONTINUE ALBVAL = -SUM/TWO C C NOW ADD IN ENERGY TERM C N1P=N1-MAXHF+ITEST N2P=N2-MAXHF+ITEST C N3P=N1-MAXHF N4P=N2-MAXHF IF(ITEST.EQ.0) GO TO 21 DO 18 I=1,N4P I3 = I+ITEST ALBVAL=ALBVAL+B(N1P,I3)*B(N2P,I3)*EIGENS(I,L1)/TWO 18 CONTINUE DO 20 I=1,ITEST DO 19 J=1,N4P J3=J+ITEST ALBVAL=ALBVAL+B(N1P,I)*B(N2P,J3)*OVRLAP(J,I)*EIGENS(J,L1)/TWO 19 CONTINUE 20 CONTINUE GO TO 22 21 IF(N1.NE.N2) GO TO 22 ALBVAL = ALBVAL+EIGENS(N1P,L1)/TWO C C NRB: C ADD-IN THE MASS-VELOCITY CONTRIBUTION (CONTINUUM-CONTINUUM ONLY) C 22 IF(IRELMV.EQ.0)GO TO 79 C TEST ZERO-OUT UNPHYSICALLY LARGE DIAGONAL ELEMENT (AS PER STG1 GEN1CC) IF(N1.EQ.N2)GO TO 79 SUM = ZERO IST = 1 POVAL2 = EIGENS(N1P,L1) + EIGENS(N2P,L2) DO 77 I=1,NIX H=HINT*DBLE(IHX(I)) MI=1 ADD = ZERO IFI = IRX(I) DO 76 J=IST,IFI JP=J+1 C WEIGHT GO TO (74,75),MI 74 MI = 2 WEIGHT=FOUR GO TO 78 75 MI=1 IF(J.EQ.IFI)GO TO 72 WEIGHT=TWO GO TO 78 72 IF(I.EQ.NIX)GO TO 73 WEIGHT=THREE GO TO 78 73 WEIGHT=ONE C ADD 78 POVAL1=POVALU(2*J) ADD=ADD+UJ(JP,I1)*UJ(JP,I2)*WEIGHT*(POVAL2+POVAL1)*POVAL1 76 CONTINUE IST = IFI+1 ADD = ADD*H/THREE SUM=SUM+ADD 77 CONTINUE C C NOW ADD IN ENERGY TERM C IF(N1.EQ.N2)SUM = SUM + EIGENS(N1P,L1) * EIGENS(N2P,L1) C SUM = -SUM*FSC*FSC/EIGHT ALBVAL = ALBVAL + SUM C C END OF CONTINUUM-CONTINUUM MASS-VELOCITY CONTRIBUTION C 79 IF(ITEST.EQ.0) RETURN C C NOW ADD IN CONTRIBUTION FROM THE BOUND ORBITALS USING C SIMPSONS RULE C DO 31 I=1,ITEST N2=MAXLG+I I2=IPOS(N2,L2) M=0 DO 80,K1=1,L2-1 DO 90,K2=K1,MAXNHF(K1) M=M+1 90 CONTINUE 80 CONTINUE M2=NCO(M+N2-L2+1) CALL CORECT(N2,L2,M2,R0,SIGMA,C1) IST = 1 R = ZERO C DO 30 J=1,NIX H=HINT*DBLE(IHX(J)) MI = 1 SUM = ZERO IFI=IRX(J) DO 29 K=IST,IFI R=R+H C C DETERMINE THE CONTRIBUTION INVOLVING THE ANGULAR MOMENTUM AND C THE POTENTIAL C RINV=ONE/R X1=UJ(K+1,I2)*(POVALU(2*K)-ALSQ*RINV*RINV) C C DETERMINE THE CONTRIBUTION INVOLVING THE SECOND DERIVATIVE C ACTING ON THE BOUND ORBITAL C PV=ZERO ALPHA=(R-R0)/SIGMA IF(-ALPHA.LT.TWELVE)PV=-(TWO*ALPHA*ALPHA-ONE)*(TWO/(SIGMA*SIGMA)) * *C1*EXP(-ALPHA*ALPHA) DO 23 JK=1,M2 IR=IRAD(L2,N2,JK) ZEE=-ZE(L2,N2,JK) ATERM=DBLE(IR*(IR-1))+TWO*ZEE*DBLE(IR)*R+ZEE*ZEE*R*R PV=PV+C(L2,N2,JK)*EXP(ZEE*R)*ATERM*R**(IR-2) 23 CONTINUE C GO TO (24,25),MI 24 MI=2 WEIGHT = FOUR GO TO 28 25 MI=1 IF(K.EQ.IFI) GO TO 26 WEIGHT = TWO GO TO 28 26 IF(J.EQ.NIX) GO TO 27 WEIGHT = THREE GO TO 28 27 WEIGHT = ONE C 28 SUM = SUM+UJ(K+1,I1)*(X1+PV)*WEIGHT 29 CONTINUE IST = IFI+1 ALBVAL=ALBVAL-SUM*H*B(N2P,I)/SIX 30 CONTINUE 31 CONTINUE C RETURN END C C********************************************************************* C SUBROUTINE PMUL(A,B,C,K,N) use param use comm_interface, only : iam,nproc IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C C C MULTIPLICATION OF VECTOR BY MATRIX OR MATRIX BY MATRIX C K=1 FOR B,C VECTORS, K=N FOR B,C MATRICES C OUTPUTS RESULTS FOR C(N,N) OR C(N,1) C C DIMENSION A(N,N),B(N,K),C(N,K),AB(N) DIMENSION A(3,3),B(3,K),C(3,K),AB(3) DO 20 J=1,K DO 20 I=1,N DO 18 KL=1,N 18 AB(KL)=A(I,KL)*B(KL,J) SUM=0.0D0 DO 19 KL=1,N 19 SUM=SUM+AB(KL) 20 C(I,J)=SUM RETURN END C===================================================================== SUBROUTINE POTF use param use comm_interface, only : iam,nproc C C THIS ROUTINE AUTOMATICALLY CALCULATES THE STATIC POTENTIAL OF THE C LOWEST POSSIBLE TARGET CONFIGURATION WITH THE GIVEN RADIAL C ORBITALS. C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM10=IDIM7*2, * KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3), PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(P75=0.75D0,P125=0.125D0,SIXT=16.0D0,TNINE=9.0D0) COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /POTVAL/ POVALU(KDIM10) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS 1001 FORMAT (///' POVALU'/) 1002 FORMAT (' ', 4E16.7) C DO 60,I=1,IRX(NIX)*2 POVALU(I)=ZERO 60 CONTINUE DO 10,I=1,LRANG1 DO 20,J=I,MAXNHF(I) NQ=IPOS(J,I) NUMORB=ISTAT(I,J) IF (NUMORB.EQ.0) THEN GOTO 20 ELSE K=IRX(NIX)+1 SUM1=ZERO SUM2=ZERO POVALU(IRX(NIX)*2)=POVALU(IRX(NIX)*2)+ONE*NUMORB DO 30,L=2*IRX(NIX)-2,2,-2 K=K-1 SUM1=SUM1+(STEP(K+1)/TWO)*(UJ(K,NQ)**2+UJ(K+1,NQ)**2) SUM2=SUM2+(STEP(K+1)/TWO)*(UJ(K,NQ)**2/XR(K)+ * UJ(K+1,NQ)**2/XR(K+1)) POVALU(L)=POVALU(L)+(ONE-SUM1+XR(K)*SUM2)*NUMORB 30 CONTINUE ENDIF 20 CONTINUE 10 CONTINUE DO 40,I=2,2*IRX(NIX),2 POVALU(I)=TWO*(DBLE(NZ)-POVALU(I)) 40 CONTINUE C NOW INTERPOLATE THE MISSING POINTS USING A 4-POINT INTERPOLATION C FORMULA.... A 3-POINT FORMULA IS USED AT THE ORIGIN POVALU(1)=P75*(DBLE(NZ)+POVALU(2))-P125*POVALU(4) IX=5 A=TNINE/SIXT B=ONE /SIXT POVALU(3)=A*(POVALU(4)+POVALU(2))-B*(POVALU(6)+TWO*DBLE(NZ)) DO 50,I=1,NIX DO 70,J=IX,IRX(I)*2-3,2 POVALU(J)=A*(POVALU(J+1)+POVALU(J-1))-B*( * POVALU(J+3)+POVALU(J-3)) 70 CONTINUE IF (I.NE.NIX) THEN C CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED K1=2*IRX(I)+1 K2=K1-2 POVALU(K1)=A*(POVALU(K1+1)+POVALU(K1-1))-B*( * POVALU(K1+3)+POVALU(K1-5)) POVALU(K2)=A*(POVALU(K2+1)+POVALU(K2-1))-B*( * POVALU(K2+2)+POVALU(K2-3)) IX=K1+2 GOTO 50 ELSE C FINAL POINT TO BE INTERPOLATED K1=IRX(NIX)*2-1 POVALU(K1)=A*(POVALU(K1+1)+POVALU(K1-1))-B*( * POVALU(K1-3)+DBLE(2*(NZ-NELC))) ENDIF 50 CONTINUE K=1 DO 80,I=2,2*IRX(NIX),2 K=K+1 POVALU(I)=POVALU(I)/XR(K) POVALU(I-1)=POVALU(I-1)/(XR(K)-STEP(K)/TWO) 80 CONTINUE C IF (NBUG5 .EQ. 1) THEN WRITE(IWRITE,1001) WRITE (IWRITE,1002) (POVALU(I), I=1,IRX(NIX)*2) END IF C RETURN END C======================================================================= SUBROUTINE RADINT(N1,L1,N2,L2,K,X) use param use comm_interface, only : iam,nproc C C EVALUATES X, THE RADIAL MULTIPOLE INTEGRAL OF ORDER K BETWEEN C TWO NUMERICAL ORBITALS SPECIFIED BY THE QUANTUM NUMBERS (N1,L1-1) C AND (N2,L2-1). C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS, * IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3), PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS X=ZERO J1=IPOS(N1,L1) J2=IPOS(N2,L2) DO 10,I=1,IRX(NIX)+1 X=X+WT(I)*UJ(I,J1)*UJ(I,J2)*XR(I)**K 10 CONTINUE RETURN END C===================================================================== SUBROUTINE RECOV1(I,IDAT,ICURR) use param use comm_interface, only : iam,nproc C C THIS ROUTINE IS CALLED ONLY IN THE CASE OF ARRAY OVERFLOW. C C I IS THE POSITION IN THE IDRAY AND IPRAY LISTS HOLDING C THE NAME OF THE RELEVANT DIMENSION PARAMETER AND C PREPROCESSOR PARAMETER. C IDAT IS THE ARRAY SIZE REQUIRED BY THE DATA C ICURR IS THE CURRENT ARRAY SIZE C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C CHARACTER*6 IDRAY CHARACTER*3 IPRAY C PARAMETER(IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO, C * IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM, C * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC) C PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4,ISDIM1=MZSET, C 1 ISDIM2=MZNCS,ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF, C 2 ITDIM6=MZCHS,IPDIM1=MZSPN,ICDIM1=MZNCF,ICDIM2=MZORB) C PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) COMMON /INFORM/ IREAD,IWRITE,IPUNCH DIMENSION IDRAY(27),IPRAY(27) DATA IDRAY/'IDIM1','IDIM2','IDIM3','IDIM4','IDIM5','IDIM6', 1 'IDIM7','IDIM8','IDIM9','IDIM10','IDIM11','IDIM12', 2 'IDIM13','IDIM14','ILDIM1','ILDIM3','ISDIM1', 3 'ISDIM2','ITDIM1','ITDIM2','ITDIM3','ITDIM6', 4 'IPDIM1','ICDIM1','ICDIM2','IRDIM1','LBUFF1'/ DATA IPRAY/'LR1','LR2','NR2','LMX','NPO','NIX', 1 'PTS','AMP','LAG','SHM','SLT','NR1', 2 'ORB','FAC','LR3','LR4','SET', 3 'NCS','TAR','NSS','CHF','CHS', 4 'SPN','NCF','ORB','RKC','BUF'/ 1000 FORMAT(/' * ARRAY OVERFLOW *'/ 1 /1X,A6,' (MZ',A3,') SHOULD BE INCREASED FROM', 2 I7,' TO AT LEAST ',I7) C 1001 FORMAT(/' PROGRAM TERMINATES IN RECOV1'/) 1002 FORMAT(/' CHECK TO SEE IF OTHER ARRAYS ARE GOING TO BE EXCEEDED') C WRITE(IWRITE,1000)IDRAY(I),IPRAY(I),ICURR,IDAT 1 WRITE(IWRITE,1001) CALL EXIT (0) C RETURN END C======================================================================= SUBROUTINE RESETL(MJS0,MINIM,ICONTN) use comm_interface, only : iam,nproc, A comm_barrier,comm_finalize IMPLICIT REAL*8 (A-H,O-Z) include 'PARAM' include 'mpif.h' C C READS STGANG DATA. CALCULATES LRANG2, THE RANGE OF C VALUES OF CONTINUUM ANGULAR MOMENTUM +1 C ARRANGES CONFIGURATIONS INTO SETS HAVING THE SAME L,S,PI C THE CONFIGURATIONS ARE ASSUMED TO BE READ IN ALREADY C GROUPED ACCORDING TO L,S,PI C FOLLOWS THE CODING OF CUSETL IN STGANG C C c INCLUDE 'PARAM' C CHARACTER*1 NUM(0:9) CHARACTER*10 FILEA CHARACTER*4 TITLE,CONT CHARACTER*3 RELOP C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ICDIM3=ICDIM2+ICDIM2-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) C COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM C COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN COMMON /NRB/MAXC,LFIXN0 COMMON /NRBREL/IRELMV COMMON /buff/status(MPI_STATUS_SIZE) C C DATA CONT/'CONT'/ DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DIMENSION LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1) DIMENSION NJCOMP(ICDIM2),LJCOMP(ICDIM2) DIMENSION NOCCSH(ICDIM1),NOCORB(ICDIM2,ICDIM1), 1 NELCSH(ICDIM2,ICDIM1),J1QNRD(ICDIM3,3,ICDIM1) DIMENSION LVAL(ILDIM2),ILVAL(ILDIM2),NSCOL(ILDIM2) DIMENSION TITLE(18) DIMENSION IIIL(128),IIIU(128),IWAVE(128) C NAMELIST/STGNX/MINLT,MAXLT,LRGLLO,LRGLUP,MAXC,LFIXN,MJS,NPTYIN X ,RELOP,NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 X ,IWAVE C 1000 FORMAT(//' ',10X,'STGANG DATA') 1001 FORMAT(' LSET ISPSET IPISET NCSET') 1002 FORMAT(10(2X,I6)) 1005 FORMAT (' HAMILTONIAN MATRICES ARE TO BE CALCULATED OVER VALUES' 1 /' OF THE TOTAL ANGULAR MOMENTUM FROM ',I2,' TO ',I2/ 2 ' AND TOTAL PARITY EVEN AND ODD IN EACH CASE'/) C 1007 FORMAT(///48X,'TARGET INPUT DATA') 1008 FORMAT(///' NUMBER OF CONFIGURATIONS =',I3) 1009 FORMAT(' NUMBER OF OCCUPIED SHELLS IN THESE CONFIGURATIONS'/30I3) 1010 FORMAT(' CONFIGURATION',I3/6X,' OCCUPIED ORBITALS ARE',32X,10I3) 1011 FORMAT(5X,' NUMBER OF ELECTRONS IN RESPECTIVE OCCUPIED SHELLS', 15X,10I3) 1012 FORMAT(5X,' COUPLING SCHEME') 1013 FORMAT(5X,3I3,6(I10,2I3)) 1014 FORMAT(/' THERE ARE',I3,' ORBITALS WHOSE N,L VALUES ARE'/ 1 6X,10(I8,I3)) 1019 FORMAT(//' ****** LRANG2 ******'/ 1 ' L+1 TAKES VALUES BETWEEN ',I3,' AND ',I3) C C THE FOLLOWING FORMAT STATEMENTS ARE TO READ THE CARD INPUT DATA. C 2000 FORMAT(12I5) 2001 FORMAT(18A4) C if(iam.eq.0)then open(73,file='sizeNX.dat',form='formatted',status='unknown') endif C C NRB NPTYMN=0 NPTYMX=1 NPTYIN=-1 IRELMV=0 IWAVE=0 C C IREAD=7 IWRITE=6 C i1=iam/100 i2=(iam-100*(iam/100))/10 i3=iam-(100*(iam/100))-i2*10 FILEA='rout2nx'//NUM(i1)//NUM(i2)//NUM(i3) C OPEN(UNIT=IREAD,FILE='dstgnx',STATUS='UNKNOWN') OPEN(UNIT=IWRITE,FILE=FILEA,STATUS='UNKNOWN') REWIND (IREAD) READ(IREAD,2001) (TITLE(I),I=1,18) IF (TITLE(1) .EQ. CONT) THEN INX2=36+(iam+500) C OPEN (INX2,FILE='tapenx2',ACCESS='SEQUENTIAL', OPEN (INX2,FILE='NX2.DAT',ACCESS='SEQUENTIAL', 1 STATUS='OLD',FORM='UNFORMATTED') REWIND (INX2) ICONTN=1 READ(INX2,ERR=301)MAXN,MAXORB,LRANG3,NSETS1,MAXNCF,NAST,MAXNST, 1 NCHSUM,MAXNCH,ISRAN3,NCFG,IRELMV 301 IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1) ELSE ICONTN=0 READ(IREAD, *) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, 1 NBUG8,NBUG9 END IF C CW WRITE(IWRITE,1000) C MAXC=1 LFIXN0=-1 IF(ICONTN.EQ.1)THEN MINLT=0 MAXLT=0 LRGLLO=0 LRGLUP=0 MJS=-1 LFIXN=-1 RELOP=' ' NBUG1=0 NBUG2=0 NBUG3=0 NBUG4=0 NBUG5=0 NBUG7=0 NBUG8=0 NBUG9=0 C READ(IREAD,STGNX) C C Reminder to self ( CPB ! insert distribution here ) C C 1/ straight split of L ( many proc option ) C C 2/ iwave (list number of L per proc - fewer proc option) C C Discussion : If the straight split of all L is carried out C no further input is required. If many processors C are available, does it matter to you whether only C a few work hard right up till the end? C C The continuum basis size for increasingly L C decreases in size. However, if maxc is large C to begin with, ie 65 , then it takes a long C time to fall off (ie maxc=42 at L=50 ), then C stick with option 1. C C BUT if you have a maxc=35, by L=50 you can have C a maxc=15 ( ie 15/35 compared with 42/65 ), C which leads to the higher partial waves finishing C a long time before the low, and ultimately poor 'load balancing' C C then option 2/ is for you C C ie in namelist IWAVE=2,3,4,5,6,6,7,7, for 8 procs C C output : file sizeNX.dat is required for pstg3nx C C IRANGE=MAXLT-MINLT+1 if((iwave(1).eq.0).and.(iwave(2).eq.0))then ISPLIT=INT(IRANGE/nproc) if(isplit.eq.0)then write(0,*)'Too many processors specified for MINLT,' x ,' MAXLT. Reduce nproc to',irange stop 'Too many processors for L-range...' endif ISPLIT2=MOD(IRANGE,nproc) KK=0 DO II=MINLT,MAXLT-ISPLIT2,ISPLIT KK=KK+1 IIIL(KK)=II IIIU(KK)=II+ISPLIT-1 ENDDO IIIL(nproc)=MAXLT-ISPLIT2-ISPLIT+1 IIIU(nproc)=MAXLT else DO II=1,nproc IIIL(II)=MINLT IIIU(II)=MINLT+iwave(ii)-1 C write(0,*)IIIL(II),IIIU(II) MINLT=IIIU(II)+1 ENDDO if(iam.eq.(nproc-1))then !Check MAXLT if(IIIU(nproc).ne.MAXLT)then write(0,*)'STG1NX maxlt=', MAXLT,' is not achieved:',IIIU(nproc) stop ' Inconsistent partial wave distribution over procs' endif endif endif C C set values for proc iam C MINLT=IIIL(iam+1) MAXLT=IIIU(iam+1) c write(0,*)'iam=',iam,' MINLT=',MINLT,' MAXLT=',MAXLT C C write PW info to sizeNX.dat C IF(iam.eq.0)then DO II=1,nproc write(73,2000)IIIL(II),IIIU(II) ENDDO ENDIF C C IF(RELOP.EQ.'MVD'.OR.RELOP.EQ.'YES')IRELMV=1 IF(RELOP.EQ.'NO')IRELMV=0 IF(NPTYIN.EQ.0)NPTYMX=0 IF(NPTYIN.EQ.1)NPTYMN=1 LFIXN0=LFIXN IF(MJS.LT.0)MJS=1 MJS0=MJS IF(MINLT.GT.0)LRGLLO=MINLT IF(MAXLT.GT.0)LRGLUP=MAXLT ELSE READ (IREAD,*) LRGLLO,LRGLUP ENDIF C CW WRITE (IWRITE,1005) LRGLLO,LRGLUP IF (ILDIM3 .LE. LRGLUP) CALL RECOV1(16,LRGLUP+1,ILDIM3) C CW WRITE(IWRITE,1007) C C READ IN THE QUANTUM NUMBERS DEFINING THE CONFIGURATIONS C IF (ICONTN .EQ. 1) THEN READ(INX2) (NOCCSH(I),I=1,NCFG) CW WRITE(IWRITE,1008)NCFG DO 101 I1=1,NCFG N1=NOCCSH(I1) READ(INX2) (NOCORB(J,I1),J=1,N1),(NELCSH(J,I1),J=1,N1) M1=N1+N1-1 READ(INX2)((J1QNRD(J,K,I1),K=1,3),J=1,M1) 101 CONTINUE READ(INX2)(NJCOMP(I),LJCOMP(I),I=1,MAXORB) REWIND (INX2) ELSE READ(IREAD,*) MAXORB,NELC,NSET,NKEY IF (ICDIM2 .LT. MAXORB) CALL RECOV1(25,MAXORB,ICDIM2) READ(IREAD,*)(NJCOMP(I),LJCOMP(I),I=1,MAXORB) IREAD1=IREAD IF (NKEY .NE. 2) THEN IF (NKEY .EQ. -2) THEN READ(IREAD,*) NOCC IF (NOCC .GT. ICDIM1) CALL RECOV1(24,NOCC,ICDIM1) READ(IREAD,*) (NOCCSH(I),I=1,NOCC) DO 91 I=1,NOCC JACT=NOCCSH(I) READ(IREAD,*) READ(IREAD,*) 91 CONTINUE ELSE READ(IREAD,*) NOPTN READ(IREAD,*) READ(IREAD,*) DO 15 M=1,NOPTN READ(IREAD,*) 15 CONTINUE END IF DO 16 N=1,NSET READ(IREAD,*) LL,LSPN,LPTY 16 CONTINUE IREAD1=36 OPEN(UNIT=36,FILE='CONF.DAT',ACCESS='SEQUENTIAL', 1 STATUS='OLD',FORM='UNFORMATTED') REWIND (36) END IF C C READ IN THE QUANTUM NUMBERS DEFINING THE CONFIGURATIONS C READ(IREAD1,*)NCFG CW WRITE(IWRITE,1008)NCFG IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1) READ(IREAD1,*) (NOCCSH(I),I=1,NCFG) CW WRITE(IWRITE,1009) (NOCCSH(I),I=1,NCFG) DO 192 I=1,NCFG N=NOCCSH(I) READ(IREAD1,*) (NOCORB(J,I),J=1,N) CW WRITE(IWRITE,1010) I,(NOCORB(J,I),J=1,N) READ(IREAD1,*) (NELCSH(J,I),J=1,N) CW WRITE(IWRITE,1011) (NELCSH(J,I),J=1,N) M=N+N-1 READ(IREAD1,*) ((J1QNRD(J,K,I),K=1,3),J=1,M) CW WRITE(IWRITE,1012) CW WRITE(IWRITE,1013)((J1QNRD(J,K,I),K=1,3),J=1,M) 192 CONTINUE IF (NKEY .NE. 2) REWIND(36) C END IF C CW WRITE(IWRITE,1001) C C IS WILL COUNT THE SETS C LRANG3 WILL CONTAIN MAX. CONFIGURATION ANGULAR MOMENTUM + 1 C IS=0 LRANG3=0 C C LOOP OVER ALL CONFIGURATIONS READ C DO 2 IC=1,NCFG ISH=NOCCSH(IC) ISH2=ISH+ISH-1 IPARTY=0 C C LOOP OVER THE SHELLS OF THIS CONFIG. TO FIND THE PARITY C DO 3 ISHEL=1,ISH IL=NOCORB(ISHEL,IC) IPARTY=IPARTY+LJCOMP(IL)*NELCSH(ISHEL,IC) 3 CONTINUE C IPARTY=MOD(IPARTY,2) LCFG=(J1QNRD(ISH2,2,IC)-1)/2 ISPIN= J1QNRD(ISH2,3,IC) C IF(IS .NE. 0) THEN IF (LCFG .EQ. LSET(IS) .AND. 1 ISPIN .EQ. ISPSET(IS) .AND. 2 IPARTY .EQ. IPISET(IS)) THEN C C FILL CURRENT SET C ICS=ICS+1 IF (ISDIM2 .LT. ICS) CALL RECOV1(18,ICS,ISDIM2) NCSET(IS)=ICS GO TO 2 END IF END IF C C START NEW SET C IS=IS+1 IF (ISDIM1 .LT. IS) CALL RECOV1(17,IS,ISDIM1) ICS=1 LSET(IS)=LCFG LRANG3=MAX(LRANG3,LCFG) ISPSET(IS)=ISPIN IPISET(IS)=IPARTY NCSET(IS)=1 C 2 CONTINUE C NSET=IS C DO 4 I=1,NSET CW WRITE(IWRITE,1002)LSET(I),ISPSET(I),IPISET(I),NCSET(I) 4 CONTINUE LRANG3=LRANG3+1 IF (ILDIM1 .LT. LRANG3) CALL RECOV1(15,LRANG3,ILDIM1) LLPDIM=LRANG3+LRANG3-1 C C TO DETERMINE WHICH SETS COUPLE WITH EACH VALUE OF THE C CONTINUUM ANGULAR MOMENTUM WITHIN THE RANGE ALLOWED C THUS DETERMINE THE RANGE REQUIRED FOR THE CONTINUUM C ANGULAR MOMENTUM C MINIM=999 LRANG2=1 DO 5 LRGL=LRGLLO,LRGLUP DO 6 NPTY0=NPTYMN,NPTYMX IF(NPTYIN.LT.0)THEN NPTY=NPTY0 ELSE NPTY=MOD(LRGL,2)+NPTY0 NPTY=MOD(NPTY,2) ENDIF C C FIND RANGE OF CONTINUUM ANGULAR MOMENTUM AND C SET UP ILVAL(KL) - THE POSSIBLE VALUES OF CONTINUUM C ANGULAR MOMENTUM C LMIN=999 LMAX=0 C DO 7 IS=1,NSET LCFG=LSET(IS) LMIN=MIN(LMIN,ABS(LCFG-LRGL)) LMAX=MAX(LMAX, LCFG) 7 CONTINUE LMAX=LMAX+LRGL ILTOT=LMAX-LMIN+1 C DO 8 KL=1,ILTOT ILVAL(KL)=LMIN+KL-1 8 CONTINUE C C LOOP OVER THIS RANGE OF CONTINUUM L C DO 9 KL=1,ILTOT L=ILVAL(KL) NL=0 C DO 10 IS=1,NSET LCFG=LSET(IS) IPI=IPISET(IS) C C PARITY TEST C IF(MOD(IPI+L+NPTY,2) .NE. 0) GO TO 10 C C TRIANGLE RULE C IF(LRGL .GT. LCFG+L .OR. 1 LRGL .LT. ABS(LCFG-L)) GO TO 10 C NL=NL+1 NSCOL(KL)=1 GO TO 9 10 CONTINUE C 9 CONTINUE C C NOW DELETE VALUES OF THE CONTINUUM ANGULAR MOMENTUM WHICH C CANNOT COUPLE WITH ANY SET AND CONSTRUCT /SETL/ C KL=0 DO 11 IKL=1,ILTOT IF (NSCOL(IKL) .NE. 0) THEN KL=KL+1 LVAL(KL)=ILVAL(IKL) END IF 11 CONTINUE LTOT=KL C IF (LTOT .EQ. 0) GO TO 6 C LMIN=LVAL(1) LMAX=LVAL(LTOT) LRANG2=MAX(LMAX+1, LRANG2) MINIM=MIN (LMIN+1, MINIM) 6 CONTINUE 5 CONTINUE CW WRITE (IWRITE,1019) MINIM, LRANG2 C RETURN END C*********************************************************************** SUBROUTINE ROOT(T,FT,B,C,RELERR,ABSERR,IFLAG) use param use comm_interface, only : iam,nproc C C THIS SUBROUTINE IS TAKEN FROM THE BOOK OF: SHAMPINE AND GORDON, C 'COMPUTATIONAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS' C C ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0, WHERE F(X) C IS A CONTINUOUS REAL FUNCTION OF A SINGLE REAL VARIABLE X. THE C METHOF OF SOLUTION IS A COMBINATION OF BISECTION AND THE SECANT RULE. C C NORMAL INPUT CONSISTS OF A CONTINUOUS FUNCTION F AND IN INTERVAL C (B,C) SUCH THAT F(B)*F(C).LE.0.0. EACH ITERATION FINDS NEW VALUES OF C B AND C SUCH THAT THE INTERVALL (B,C) IS SHRUNK AND F(B)*F(C).LE.0.0. C C THE STOPPING CRITERION IS C C ABS(B-C).LE.2.0*(RELERR*ABS(B)+ABSERR) C C WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE INPUT C QUANTITIES. SET THE FLAG, IFLAG, POSITIVE TO INITIALISE THE C COMPUTATION. AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND OUTPUT, C THEY MUST BE VARIABLES IN THE CALLING PROGRAM. C IF 0.0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0. C C C THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT AS C B AND C ARE ALWAYS REDEFINED SO THAT ABS(F(B)).LE.ABS(F(C)). C C TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS IS C DONE IN THE CALLING PROGRAM. WHEN AN EVALUATION OF F IS NEEDED AT T, C ROOT RETURNS WITH IFLAG NEGATIVE. EVALUATE FT=F(T) AND CALL ROOT C AGAIN. DO NOT ALTER IFLAG. C C WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING C PROGRAM WITH IFLAG POSITIVE: C C IFLAG=1 IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET. C C =2 IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE F(B) C IS EXACTLY ZERO. THE INTERVAL (B,C) MAY NOT SATISFY C THE STOPPING CRITERION. C C =3 IF ABS(F(B)) EXCEEDS THE INPUT VALUES ABS(F(B)), C ABS(F(C)). IN THIS CASE IT IS LIKELY THAT B IS CLOSE C TO A POLE OF F. C C =4 IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVALL. C A LOCAL MININMUM MAY HAVE BEEN OBTAINED. C C =5 IF TOO MANY FUNCTION EVALUATIONS WERE MADE. C (AS PROGRAMMED, 500 ARE ALLOWED.) C C THIS CODE IS A MODIFICATION OF THE CODE Z E R O I N WHICH IS C COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, NUMERICAL C COMPUTING: AN INTRODUCTION BY L.F. SHAMPINE AND R.C. ALLEN C C*********************************************************************** C THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT * C ROUNDOFF ERROR U. IT IS CALCULATED IN THE SUBROUTINE M A C H I N * C WHICH MUST HAVE BEEN CALLED BEFORE THE FIRST CALL OF R O O T. * C NOTE....IN THIS CODE MACHIN IS NOT CALLED BUT U IS SET TO 1.0E-12 * C IN THE PARAMETER STATEMENT BELOW. * C*********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER (ZERO=0.0D0,ONE=1.0D0,EIGHT=8.0D0,HALF=0.5D0) PARAMETER (U=1.0D-12) C SAVE RE,AE,IC,ACBS,A,FA,FC,FX,KOUNT !NRB C IF (IFLAG.GE.0) GOTO 100 IFLAG=ABS(IFLAG) GOTO (200,300,400), IFLAG 100 RE=MAX(RELERR,U) AE=MAX(ABSERR,ZERO) IC=0 ACBS=ABS(B-C) A=C T=A IFLAG=-1 RETURN 200 FA=FT T=B IFLAG=-2 RETURN 300 FB=FT FC=FA KOUNT=2 FX=MAX(ABS(FB),ABS(FC)) 1 IF(ABS(FC).GE.ABS(FB)) GOTO 2 C C INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C)) C A=B FA=FB B=C FB=FC C=A FC=FA 2 CMB=HALF*(C-B) ACMB=ABS(CMB) TOL=RE*ABS(B)+AE C C TEST STOPPING CRITERION AND FUNCTION COUNT C IF(ACMB.LE.TOL) GOTO 8 IF(KOUNT.GE.500) GOTO 12 C C CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q WHERE WE ARRANGE P.GE.0. C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW. C P=(B-A)*FB Q=FA-FB IF(P.GE.ZERO) GOTO 3 P=-P Q=-Q C C UPDATE A, CHECK REDUCTION OF THE SIZE OF BRACKETING INTERVAL IS C SATISFACTORY. IF NOT, BISECT UNTIL IT IS. C 3 A=B FA=FB IC=IC+1 IF(IC.LT.4) GOTO 4 IF(EIGHT*ACMB.GE.ACBS) GOTO 6 IC=0 ACBS=ACMB C C TEST FOR TOO SMALL A CHANGE C 4 IF(P.GT.ABS(Q)*TOL) GOTO 5 C C INCREMENT BY TOLERANCE C B=B+SIGN(TOL,CMB) GOTO 7 C C ROOT OUGHT TO BE BETWEEN B AND (C+B)/2 C 5 IF(P.GE.CMB*Q) GOTO 6 C C USE SECANT RULE C B=B+P/Q GOTO 7 C C USE BISECTION C 6 B=HALF*(C+B) C C HAVE COMPLETED COMPUTATION OF NEW ITERATE B C 7 T=B IFLAG=-3 RETURN 400 FB=FT IF(FB.EQ.ZERO) GOTO 9 KOUNT=KOUNT+1 IF(SIGN(ONE,FB).NE.SIGN(ONE,FC)) GOTO 1 C=A FC=FA GOTO 1 C C FINISHED. SET IFLAG C 8 IF(SIGN(ONE,FB).EQ.SIGN(ONE,FC)) GOTO 11 IF(ABS(FB).GT.FX) GOTO 10 IFLAG=1 RETURN 9 IFLAG=2 RETURN 10 IFLAG=3 RETURN 11 IFLAG=4 RETURN 12 IFLAG=5 RETURN END C======================================================================= SUBROUTINE RS(N2,L2,N1,L1,N4,L4,N3,L3,K,RKVAL) use param use comm_interface, only : iam,nproc C C THIS SUBROUTINE EVALUATES THE FUNCTION Y(P2,P4,K/R) AND THE C SLATER INTEGRAL R(P1,P2,P3,P4,K), WHERE P1,P2,P3,P4 ARE NUMERICAL C ORBITAL FUNCTIONS DEFINED OVER A PRESCRIBED SET OF INTEGRATION C MESH POINTS. C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM3=MZNR2, * IDIM7=MZPTS, * IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS COMMON /YKSTOR/ YK(IDIM7),RK(IDIM7),RK1(IDIM7),JN2,JL2,JN4,JL4,JKM DIMENSION SUM1(IDIM7),SUM2(IDIM7),RINT(IDIM7),AR(IDIM7) C C ******** NOTES FOR THE USER ******** C C THE FUNCTIONS P1,P2,P3,P4 ARE TO BE STORED IN THE ARRAY UJ(J,I) C WITH I=1,4 RESPECTIVELY.....J IS THE INTEGRATION POINT NUMBER C WITH J=1 CORRESPONDING TO R=0.0 C C ******** ******** C C INITIALIZATION OF CONSTANTS AND ARRAY INDICES C C LOCATE THE POSITIONS OF THE P1,P2,P3,P4 FUNCTIONS IN THE UJ ARRAY C FROM THEIR N,L VALUES SPECIFIED IN THE ARGUMENT LIST C TWELTH=ONE/TWELVE J1=IPOS(N1,L1) J2=IPOS(N2,L2) J3=IPOS(N3,L3) J4=IPOS(N4,L4) C K1=K+1 C C DO NOT RECALCULATE THE YK INTEGRAL IF IT ALREADY EXISTS C IF(K-JKM) 2,1,2 1 IF(N2.NE.JN2.OR.L2.NE.JL2) GO TO 2 IF(N4.NE.JN4.OR.L4.NE.JL4) GO TO 2 GO TO 7 2 JN2=N2 JL2=L2 JN4=N4 JL4=L4 IF(K.EQ.JKM) GO TO 22 DO 20 J=2,NPTS 20 RK(J)=XR(J)**K DO 21 J=2,NPTS 21 RK1(J)=ONE/(RK(J)*XR(J)) JKM=K C C INITIALIZATION C 22 ARA=ZERO DO 3 J=2,NPTS SUM1(J)=UJ(J,J2)*UJ(J,J4) 3 CONTINUE DO 4 J=2,NPTS RINT(J)=SUM1(J)*RK(J) AR(J)=SUM1(J)*RK1(J) 4 CONTINUE C C SIMPSONS RULE INTEGRATION FROM 0 TO RA C DO 5 J=2,NPTS ARA=ARA+(WT(J)*AR(J)) 5 CONTINUE C C SIMPSONS RULE INTEGRATIONS FROM 0 TO R C YK(1)=ZERO RINT(1)=ZERO AR(1)=ZERO RIN=ZERO DO 8 J=3,NPTS,2 SUM1(J) =STEP(J) *THIRD *(RINT(J-2)+FOUR*RINT(J-1)+RINT(J)) SUM2(J) =STEP(J) *THIRD *(AR(J-2)+FOUR*AR(J-1)+AR(J)) 8 CONTINUE DO 11 J=3,NPTS,2 SUM1(J-1)=STEP(J-1)*TWELTH*(FIVE*RINT(J-2)+EIGHT*RINT(J-1)-RINT(J) 1) SUM2(J-1)=STEP(J-1)*TWELTH*(FIVE*AR(J-2)+EIGHT*AR(J-1)-AR(J)) 11 CONTINUE DO 9 J=3,NPTS,2 RINT(J-1)=RIN+SUM1(J-1) RIN=RIN+SUM1(J) RINT(J)=RIN AR(J-1)=ARA-SUM2(J-1) ARA=ARA-SUM2(J) AR(J)=ARA 9 CONTINUE DO 10 J=2,NPTS YK(J)=(RINT(J)*RK1(J)+AR(J)*RK(J))*WT(J) 10 CONTINUE C C EVALUATE THE SLATER INTEGRAL C 7 CONTINUE RKVAL=ZERO DO 12 J=2,NPTS RKVAL=RKVAL+UJ(J,J1)*UJ(J,J3)*YK(J) 12 CONTINUE RETURN END C======================================================================= SUBROUTINE SCHMDT(LP) use param use comm_interface, only : iam,nproc C C SCHMIDT ORTHOGONALIZE THE CONTINUUM ORBITALS TO THE BOUND C ORBITALS WHICH ARE NOT INCLUDED IN THE ORTHOGONALIZATION IN C BASFUN FOR ANGULAR MOMENTUM LP-1. C THE SCHMIDT COEFFICIENTS ARE STORED IN THE B-ARRAY C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9, * KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5) C C N1 IS THE NUMBER OF SCHMIDT ORTHOGONALIZED BOUND ORBITALS AND C NPTS IS THE NUMBER OF INTEGRATION POINTS C MAXHF=MAXNHF(LP) MAXLG=MAXNLG(LP) N1=MAXHF-MAXLG NPTS=IRX(NIX)+1 C C ZEROIZE THE B-ARRAY AND SET THE REQUIRED DIAGONAL ELEMENTS TO C UNITY C N=NRANG2+N1 DO 2 I=1,N DO 1 J=1,N B(I,J)=ZERO 1 CONTINUE 2 CONTINUE DO 3 I=1,N1 B(I,I)=ONE 3 CONTINUE C C CALCULATE THE SCHMIDT COEFFICIENTS C DO 8 I=1,NRANG2 N2=N1+I-1 ANORM=ONE DO 5 J=1,N2 TEMP(J)=ZERO DO 4 K=1,N1 TEMP(J)=TEMP(J)-B(J,K)*OVRLAP(I,K) 4 CONTINUE ANORM=ANORM-TEMP(J)**2 5 CONTINUE C C RETURN IF NRANG2 IS TOO LARGE C IF (ANORM .LE. TINY) THEN LP=-I GO TO 15 END IF ANORM=ONE/SQRT(ANORM) DO 7 J=1,N2 B(N2+1,J)=ZERO DO 6 K=1,N2 B(N2+1,J)=B(N2+1,J)+TEMP(K)*B(K,J) 6 CONTINUE B(N2+1,J)=B(N2+1,J)*ANORM 7 CONTINUE B(N2+1,N2+1)=ANORM 8 CONTINUE C C SCHMIDT ORTHOGONALIZE THE CONTINUUM ORBITALS C I1=IPOS(MAXHF+1,LP)-1 DO 14 I=1,NPTS DO 9 J=1,NRANG2 I2=I1+J TEMP(J)=UJ(I,I2) 9 CONTINUE DO 13 J=1,NRANG2 I2=I1+J SUM=ZERO N2=N1+J DO 12 K=1,N2 IF(K.LE.N1)GO TO 10 I3=K-N1 X1=TEMP(I3) GO TO 11 10 K1=K+MAXLG I3=IPOS(K1,LP) X1=UJ(I,I3) 11 SUM=SUM+X1*B(N2,K) 12 CONTINUE UJ(I,I2)=SUM 13 CONTINUE 14 CONTINUE C 15 RETURN END C======================================================================= SUBROUTINE SHRIEK(NFACT) use param use comm_interface, only : iam,nproc IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM14=MZFAC) COMMON /FACT / GAMMA(IDIM14) GAMMA(1)=DBLE(1) IF (NFACT.LT.2) RETURN DO 1,I=2,NFACT II=I-1 GAMMA(I)=DBLE(II)*GAMMA(II) 1 CONTINUE RETURN END C======================================================================= SUBROUTINE SKIPER(L1,L2) use param use comm_interface, only : iam,nproc C C ROUTINE FOR READING THE CONTINUUM ORBITALS AND THE SCHMIDT C COEFFICIENTS CORRESPONDING TO THE ANGULAR MOMENTA (L1-1) AND C (L2-1) FROM THE SCRATCH DISC IDISC1 INTO THE UJ ARRAY. IF ONLY C ONE ANGULAR MOMENTUM IS REQUIRED THEN L2 MUST BE EQUAL TO ZERO C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM4=IDIM3+IDIM10,KDIM5=KDIM4+IDIM9, * KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SCOEFF/ B(KDIM4,KDIM4),OVRLAP(IDIM3,IDIM10),TEMP(KDIM5) IF (L1.LE.0) THEN PRINT*,' ===== ERROR IN SKIPER ====' STOP ENDIF NPTS=IRX(NIX)+1 MAXHF=MAXNHF(L1) MAXLG=MAXNLG(L1) ITEST=MAXHF-MAXLG N1=IDIM13+1 N2=N1+NVARY(L1)-1 READ(61,REC=L1) ((UJ(I,N),I=1,NPTS),N=N1,N2) DO 10,N=1,NVARY(L1) IPOS(MAXHF+N,L1)=N1+N-1 10 CONTINUE IF (ITEST.GT.0) THEN READ(62,REC=L1) B,OVRLAP ENDIF IF (L2.EQ.0) THEN RETURN ENDIF MAXHF=MAXNHF(L2) MAXLG=MAXNLG(L2) ITEST=MAXHF-MAXLG N1=N2+1 N2=N1+NVARY(L2)-1 READ(61,REC=L2) ((UJ(I,N),I=1,NPTS),N=N1,N2) DO 12,N=1,NVARY(L2) IPOS(MAXHF+N,L2)=N1+N-1 12 CONTINUE IF (ITEST.GT.0) THEN READ(62,REC=L2) B,OVRLAP ENDIF RETURN END C********************************************************************* SUBROUTINE STGRRD (MJS) use param use comm_interface, only : iam,nproc C THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES C CERTAIN VARIABLES IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C C CHARACTER*8 FILE1,FILE2,FILE3 CHARACTER*11 FF1,FF2,FF3 CHARACTER*72 LVALUE(4)*1 character*1 NUM(0:9) DATA LVALUE/'S','P','D','F'/ DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB,IDIM14=MZFAC) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM9=IDIM1+IDIM1-1, * KDIM12=IDIM1*IDIM12) PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9) C COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1, * JBUFF1,JBUFF2 COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /COPY / ICOPY1,ICOPY2,ITOTAL,IPSEUD COMMON /CORE / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2), * MAXNC(IDIM2),MAXPN(IDIM2),ICHECK COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4 COMMON /FACT / GAMMA(IDIM14) COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS C C SET DEFAULT VALUE OF LFIXN, THE L VALUE AT WHICH THE C NUMBER OF CONTINUUM TERMS IS ALLOWED TO DECREASE BELOW C NRANG2. DATA LNDEF/10/ C C FORMAT STATEMENTS 2000 FORMAT(A) 1000 FORMAT(TR31,'***********'/TR31,'** **'/TR31,'** NXRAD **'/ * TR31,'** **'/TR31,'***********'//) 1002 FORMAT(////TR28,' SUBROUTINE STGRRD'/TR28,' -----------------') 1003 FORMAT(//' OUTPUT CHANNEL UNIT NUMBER =',I4/ * ' INPUT CHANNEL UNIT NUMBER =',I4/ * ' INPUT SEQUENTIAL UNIT NUMBER =',I4/ * ' (FOR MODEL POTENTIAL)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I4, * ' NAME= ',A8/ * ' (FOR RADIAL INTEGRALS)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I4, * ' NAME= ',A8/ * ' (FOR RK POINTER ARRAYS)'/ * ' OUTPUT SEQUENTIAL UNIT NUMBER =',I4, * ' NAME= ',A8/) 1004 FORMAT(' NBUG PARAMETERS'/12I5) 1005 FORMAT(' === NON-MODEL POTENTIAL CALCULATION ===') 1006 FORMAT(' === MODEL POTENTIAL CALCULATION ===') 1007 FORMAT(' === ANALYTIC BOUND ORBITALS BEING USED ===') 1008 FORMAT(' === NUMERICAL BOUND ORBITALS BEING USED ===') 1009 FORMAT(/' BASIC DATA'/) 1010 FORMAT(' NELC =',I3,' NZ =',I3,' LRANG1 =',I3,' NRANG2 =',I3, * ' LFIXN =',I3, ' LAMAX =',I2,' IBC =',I2) 1011 FORMAT(' MAXNHF=',35I3/) 1012 FORMAT(' MAXNLG=',35I3/) 1013 FORMAT(/' DATA DEFINING THE STATIC POTENTIAL -', 1 ' SHELL OCCUPANCY OF THE GROUND STATE') 6014 FORMAT(1X,35(1X,I1,A1)) 1014 FORMAT(1X,35I3) 1015 FORMAT(//' THE RADIAL FUNCTIONS'/TR37,' SLATER-TYPE',TR5, *' POWER OF R',TR5,' EXPONENT') 1016 FORMAT(/TR13,I2,A,' ORBITAL') 1017 FORMAT(TR37,E12.5,TR9,I3,TR9,E12.5) 1018 FORMAT(//' R-MATRIX BOUNDARY CONDITIONS, RA =',F10.5,' BSTO = *',F10.5//' AMPLITUDE OF THE FUNCTIONS AT RA'/) 1019 FORMAT(//' INTEGRATION MESH'//' NIX=',I3) 1020 FORMAT(' HINT =',E14.7) 1021 FORMAT(' IHX=',20I5) 1022 FORMAT(' IRX=',20I5//) 1023 FORMAT(' MAXNC=',20I5) 1024 FORMAT(' LPOT=',I3) 1025 FORMAT(' LPOS=',20I5) 1026 FORMAT(I5,A,TR4,E14.7) 2005 FORMAT(' *** ANGULAR AND RADIAL DATA INCONSISTENT'/ 1 ' NELC =',I3,' IN CONFIGURATIONS AND =',I3, 2 ' IN RADIAL DATA') C ETA=0.0D0 C C SET THE CHANNEL UNIT NUMBERS C IREAD=7 IWRITE=6 JBUFF1=23 JBUFF2=24 ITAPE3=25 ITAPE1=35 C i1=iam/100 i2=(iam-100*(iam/100))/10 i3=iam-(100*(iam/100))-i2*10 C FF1='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FF2='RAD2.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FF3='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) C C OPEN THE DIRECT ACCESS FILES OPEN(UNIT=JBUFF1,FILE=FF1,ACCESS='DIRECT', * STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1) OPEN(UNIT=JBUFF2,FILE=FF2,ACCESS='DIRECT', * STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2) OPEN(ITAPE3,FILE=FF3,ACCESS='SEQUENTIAL', * STATUS='UNKNOWN',FORM='UNFORMATTED') C READ IN THE BASIC DATA WRITE(IWRITE,1000) WRITE(IWRITE,1002) WRITE(IWRITE,1003)IWRITE,IREAD,ITAPE1,JBUFF1,FILE1, * JBUFF2,FILE2,ITAPE3,FILE3 READ (IREAD, *) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, * NBUG8,NBUG9 WRITE(IWRITE,1004) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, * NBUG8,NBUG9 READ (IREAD, *) IPSEUD,NUMERC IF (IPSEUD.EQ.0) THEN WRITE(IWRITE,1005) ELSE WRITE(IWRITE,1006) C EFFECTIVE POTENTIAL (2*ZA/R) IPSEUD=IABS(IPSEUD) ENDIF IF (NUMERC.EQ.0) THEN WRITE(IWRITE,1007) ELSE WRITE(IWRITE,1008) ENDIF WRITE(IWRITE,1009) READ (IREAD, *) NELC1,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC IF (NELC1 .NE. NELC) THEN WRITE (IWRITE,2005) NELC, NELC1 STOP END IF IF (LRANG1 .GT. IDIM1) CALL RECOV1(1,LRANG1,IDIM1) IF (NRANG2 .GT. IDIM3) CALL RECOV1(3,NRANG2,IDIM3) IF (LAMAX .GT. IDIM4) CALL RECOV1(4,LAMAX,IDIM4) C C CONSTRUCT NVARY(L) C IF (LFIXN .LT. LNDEF) THEN LFIXN=LNDEF WRITE(IWRITE,'('' LFIXN HAS BEEN RESET TO '',I3)')LFIXN END IF C DO 1 I=1,LFIXN NVARY(I)=NRANG2 1 CONTINUE C DO 2 I=LFIXN+1,LRANG2 NVARY(I)=0 2 CONTINUE C WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,IBC READ (IREAD, *) (MAXNHF(I),I=1,LRANG1) WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1) READ (IREAD, *) (MAXNLG(I),I=1,LRANG1) WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1) C READ IN DATA DEFINING THE STATIC POTENTIAL OF THE GROUND STATE WRITE(IWRITE,1013) NBOUND=0 DO 10,I=1,LRANG1 NHF=MAXNHF(I) NLG=MAXNLG(I) IF (IDIM9 .LE. NLG-I) CALL RECOV1(9,NLG-I+1,IDIM9) IF (IDIM10 .LT. NHF-NLG) CALL RECOV1(10,NHF-NLG,IDIM10) IF (IDIM12 .LT. NHF) CALL RECOV1(12,NHF,IDIM12) NBOUND=NBOUND+NHF-I+1 READ (IREAD, *) (ISTAT(I,J),J=I,MAXNHF(I)) WRITE(IWRITE,6014) (J,LVALUE(I),J=I,MAXNHF(I)) WRITE(IWRITE,1014) (ISTAT(I,J),J=I,MAXNHF(I)) 10 CONTINUE IF (IDIM13 .LT. NBOUND) CALL RECOV1(13,NBOUND,IDIM13) IF (NUMERC.EQ.0) THEN C ANALYTIC ORBITALS ARE TO BE USED M1=0 ZMIN=ONE/TINY WRITE(IWRITE,1015) C READ IN COEFFICIENTS DEFINING THE ANALYTIC ORBITALS DO 20,L=1,LRANG1 DO 30,N=L,MAXNHF(L) WRITE(IWRITE,1016) N,LVALUE(L) M1=M1+1 READ (IREAD, *) M IF (IDIM11 .LT. M) CALL RECOV1(11,M,IDIM11) READ (IREAD, *) (IRAD(L,N,K),K=1,M) READ (IREAD, *) (ZE (L,N,K),K=1,M) READ (IREAD, *) (C (L,N,K),K=1,M) DO 31, K=1,M IR=IRAD(L,N,K) IR2=IR+IR+1 IF (IDIM14 .LT. IR2)CALL RECOV1(14,IR2,IDIM14) 31 CONTINUE NCO(M1)=M X=ZERO C CHECK THE NORMALISATION OF EACH ANALYTIC ORBITAL DO 40,K=1,M IF (ZE(L,N,K).LT.ZMIN) THEN ZMIN=ZE(L,N,K) IMIN=IRAD(L,N,K) ENDIF X1=ZERO K2=IRAD(L,N,K)+1 Z1=ZE(L,N,K) DO 180,K1=1,M X1=X1+C(L,N,K1)*GAMMA(IRAD(L,N,K1)+K2)/ * (ZE(L,N,K1)+Z1)**(K2+IRAD(L,N,K1)) 180 CONTINUE X=X+C(L,N,K)*X1 40 CONTINUE IF (ABS(X-ONE).GT.PT01) THEN C THE COEFFICIENTS MAY BE CLEMENTI TYPE C RENORMALISE THEM DO 50,K=1,M IR=IRAD(L,N,K) ZR=TWO*ZE(L,N,K) C(L,N,K)=C(L,N,K)*SQRT(ZR/GAMMA(IR+IR+1))*ZR**IR 50 CONTINUE C CHECK NORMALISATION X=ZERO DO 190,K=1,M X1=ZERO K2=IRAD(L,N,K)+1 Z1=ZE(L,N,K) DO 182,K1=1,M X1=X1+C(L,N,K1)*GAMMA(IRAD(L,N,K1)+K2)/ * (ZE(L,N,K1)+Z1)**(K2+IRAD(L,N,K1)) 182 CONTINUE X=X+C(L,N,K)*X1 190 CONTINUE IF (ABS(X-ONE).GT.PT01) THEN C ERROR IN INPUT DATA. ORBITAL NOT NORMALISED PRINT*,' ORBITAL',N,LVALUE(L),'HAS NORMALISATION',X STOP ENDIF ENDIF Y=SQRT(X) DO 60,K=1,M C(L,N,K)=C(L,N,K)/Y WRITE(IWRITE,1017) C(L,N,K),IRAD(L,N,K),ZE(L,N,K) 60 CONTINUE 30 CONTINUE 20 CONTINUE IF (IBC.EQ.0) THEN C ESTIMATE THE BOUNDARY RADIUS RA ITEST=0 TINORB=PT001 RA=-LOG(TINORB)/ZMIN RA=RA+IMIN*LOG(RA)/ZMIN IRA=INT(RA+HALF) 220 RA=DBLE(IRA) BSTO=ZERO M1=0 SMALL=FIFTH C CHECK THAT THE BOUNDARY RADIUS IS CORRECT FOR ALL C ORBITALS. DO 70,L=1,LRANG1 DO 80,N=L,MAXNHF(L) M1=M1+1 M=NCO(M1) 100 X=ZERO DO 90,K=1,M X=X+C(L,N,K)*RA**IRAD(L,N,K)*EXP(-ZE(L,N,K)*RA) 90 CONTINUE IF (X.GT.TINORB) THEN RA=RA+SMALL ITEST=1 GOTO 100 ENDIF 80 CONTINUE 70 CONTINUE IF (ITEST.EQ.0) THEN IF (IRA.GE.1) THEN IRA=IRA-1 GO TO 220 ENDIF ENDIF ELSE READ (IREAD, *) RA,BSTO ENDIF WRITE(IWRITE,1018) RA,BSTO M1=0 DO 110,L=1,LRANG1 DO 120,N=L,MAXNHF(L) M1=M1+1 X=ZERO DO 130,K=1,NCO(M1) X=X+C(L,N,K)*RA**IRAD(L,N,K)*EXP(-ZE(L,N,K)*RA) 130 CONTINUE WRITE(IWRITE,1026) N,LVALUE(L),X 120 CONTINUE 110 CONTINUE ENDIF IF (IPSEUD.EQ.0) THEN C THIS IS A NON-MODEL POTENTIAL CALCULATION. IF (IBC.GT.1) THEN C READ IN INTEGRATION MESH READ (IREAD, *) NIX READ (IREAD, *) HINT READ (IREAD, *) (IHX(I),I=1,NIX) READ (IREAD, *) (IRX(I),I=1,NIX) ELSE C CALCULATE THE INTEGRATION MESH CALL MESH ENDIF C INITIALISE LPOT,LPOSX AND MAXNC IN NON-MODEL POTENTIAL C CALCULATION. LPOT=1 DO 150,I=1,LRANG2 MAXNC(I)=I-1 LPOSX(I)=1 150 CONTINUE ELSE C READ IN THE MESH AND MODEL POTENTIAL FROM THE C SEQUENTIAL FILE ITAPE1. OPEN(UNIT=ITAPE1,ACCESS='SEQUENTIAL',FORM='UNFORMATTED', * FILE='STG1.POT',STATUS='OLD') READ(ITAPE1) LRANG1,LRANG2,(MAXNC(I),I=1,LRANG1) READ(ITAPE1) HINT,NIX,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX) READ(ITAPE1) LPOT,(LPOSX(I),I=1,LPOT) NPTS=IRX(NIX) DO 140,I=1,LPOT READ(ITAPE1) (POTHAM(I,J),J=1,NPTS) 140 CONTINUE ENDIF WRITE(IWRITE,1019) NIX IF (IDIM6 .LT. NIX) CALL RECOV1(6,NIX,IDIM6) WRITE(IWRITE,1020) HINT WRITE(IWRITE,1021) (IHX(I),I=1,NIX) WRITE(IWRITE,1022) (IRX(I),I=1,NIX) IF (IDIM7 .LE. IRX(NIX))CALL RECOV1(7,IRX(NIX)+1,IDIM7) IF (IPSEUD .NE. 0) THEN WRITE(IWRITE,1023) (MAXNC(I),I=1,LRANG1) WRITE(IWRITE,1024) LPOT IF (IDIM8 .LT. LPOT) CALL RECOV1(8,LPOT,IDIM8) WRITE(IWRITE,1025) (LPOSX(I),I=1,LPOT) END IF C READ THE MAXIMUM ENERGY FOR THE BUTTLE CORRECTION C AND THE PARAMETER FOR BUTTLE FITTING C MJS=0 FOR QUADRATIC MJS NOT 0 FOR MJS METHOD C NB QUADRATIC IS ALWAYS USED IF BSTO IS NON ZERO C READ(IREAD, *)EK2MAX,MJS WRITE(IWRITE,'(/'' ENERGY LIMITS FOR THE BUTTLE CORRECTION''/ 1 '' MAXIMUM ENERGY = '',F14.7)')EK2MAX C IF (MJS .NE. 0) THEN WRITE(IWRITE,'('' MJS FITTING PROCEDURE REQUESTED'')') IF (BSTO .NE. 0) THEN WRITE(IWRITE,'('' BUT QUADRATIC FIT WILL BE USED SINCE'', 1 '' BSTO IS NON-ZERO'')') MJS=0 END IF ELSE WRITE(IWRITE,'('' QUADRATIC FITTING REQUESTED'')') END IF C C DEFINE THE MAXPN ARRAY . NOT USED? DO 160,I=1,LRANG1 MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1 160 CONTINUE C EXTEND THE MAXNHF AND MAXNLG ARRAYS IF (LRANG2.GT.LRANG1) THEN DO 170,I=LRANG1+1,LRANG2 MAXNHF(I)=I-1 MAXNLG(I)=I-1 170 CONTINUE ENDIF C INITIALISE THE VARIABLES LAMBC AND LAMCC IF (LAM.GE.2) THEN LAMBC=LAMAX ELSE LAMBC=0 ENDIF IF (LAM.GE.3) THEN LAMCC=LAMAX ELSE LAMCC=0 ENDIF C SET UP ARRAYS FOR USE IN SIMPSONS RULE INTEGRATION ONE3 =ONE /THREE TWO3 =TWO *ONE3 FOUR3=FOUR*ONE3 IFI=1 XR(1) =ZERO STEP(1)=ZERO DO 200,I=1,NIX H =HINT*DBLE(IHX(I)) WT(IFI) =ONE3*(H+STEP(IFI)) RSTART=XR(IFI) IST =IFI+2 IFI =IRX(I)+1 IFX =0 IF (I.GT.1) IFX=IRX(I-1) DO 210,J=IST,IFI,2 XR(J-1) =RSTART+(J-2-IFX)*H XR(J ) =RSTART+(J-1-IFX)*H STEP(J-1)=H STEP(J )=H WT(J-1 )=FOUR3*H WT(J )=TWO3*H 210 CONTINUE 200 CONTINUE WT(IFI)=ONE3*H NPTS=IFI RETURN END C======================================================================= SUBROUTINE STGRRX (MINIM,MJS) use param use comm_interface, only : iam,nproc C C THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES C CERTAIN VARIABLES USING RMATRX OUTPUT FILE IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C CHARACTER*11 FF1,FF2,FF3 CHARACTER*1 NUM(0:9) C CHARACTER*8 FILE1,FILE2,FILE3 C CHARACTER*72 LVALUE(4)*1 C DATA LVALUE/'S','P','D','F'/ C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX, * IDIM6=MZNIX,IDIM7=MZPTS,IDIM8=MZAMP,IDIM9=MZLAG,IDIM10=MZSHM, * IDIM11=MZSLT,IDIM12=MZNR1,IDIM13=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM1=IDIM13+2*IDIM3, * KDIM9=IDIM1+IDIM1-1,KDIM10=IDIM7*2, * KDIM12=IDIM1*IDIM12) PARAMETER(IADIM=IDIM2-1+IDIM3,IBDIM=IDIM12+IDIM3) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) C KDIM2=MAX(IDIM2-1+IDIM3,IDIM12+IDIM3) PARAMETER(KDIM2=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(LBUFF1=MZBUF,LBUFF2=IDIM1*IDIM1*KDIM9) PARAMETER(P75=0.75D0,P125=0.125D0,SIXT=16.0D0,TNINE=9.0D0) COMMON /BNDORB/ P(IDIM9,KDIM10) COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON /BUFFER/ BUFF1(LBUFF1),BUFF2(LBUFF2),NBUFF1,NBUFF2,IBUFF1, * JBUFF1,JBUFF2 COMMON/BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN COMMON /CONSTS/ ZERO,ONE,PT01,PT001,PT0001,TINY,PI ,FSC, * TWO ,THREE,FOUR,FIVE ,SIX ,EIGHT ,TEN,TWELVE, * HALF,THIRD,FOURTH,FIFTH,SIXTH ,EIGHTH,TENTH COMMON /COPY / ICOPY1,ICOPY2,ITOTAL,IPSEUD COMMON /CORE / POTHAM(IDIM8,IDIM7),LPOT,LPOSX(IDIM2), * MAXNC(IDIM2),MAXPN(IDIM2),ICHECK COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4 COMMON /INFORM/ IREAD,IWRITE,IPUNCH COMMON /INIT / HINT,IHX(IDIM6),IRX(IDIM6),NIX,IMATCH COMMON /MORDER/ LAMAX,LAMBC,LAMCC,LAMIND,LAM COMMON /NBUG / NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, * NBUG9 COMMON /NUMORB/ UDP(KDIM1,IDIM7),PORI(IDIM2),DPORI(IDIM2),NUMERC COMMON /ORBTLS/ UJ(IDIM7,KDIM1),IPOS(KDIM2,IDIM2) COMMON /POTVAL/ POVALU(KDIM10) COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM COMMON /SIMP / XR(IDIM7),STEP(IDIM7),WT(IDIM7),NPTS COMMON /NRB/MAXC,LFIXN0 DIMENSION PTEMP(IDIM7),UTEMP(IDIM7),IRXP(IDIM6) C C SET DEFAULT VALUE OF LFIXN, THE L VALUE ABOVE WHICH THE C NUMBER OF CONTINUUM TERMS IS ALLOWED TO DECREASE BELOW C NRANG2. DATA LNDEF/10/ C C FORMAT STATEMENTS 5002 FORMAT (' ', 4E16.7) 2000 FORMAT(A) 1000 FORMAT(TR31,'***********'/TR31,'** **'/TR31,'** NXRAD **'/ * TR31,'** **'/TR31,'***********'//) 1002 FORMAT(////TR28,' SUBROUTINE STGRRX'/TR28,' -----------------') 1003 FORMAT(//' OUTPUT CHANNEL UNIT NUMBER =',I4/ * ' INPUT CHANNEL UNIT NUMBER =',I4/ * ' INPUT SEQUENTIAL UNIT NUMBER =',I4/ * ' (FOR MODEL POTENTIAL)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I4, * ' NAME= ',A8/ * ' (FOR RADIAL INTEGRALS)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I4, * ' NAME= ',A8/ * ' (FOR RK POINTER ARRAYS)'/ * ' OUTPUT SEQUENTIAL UNIT NUMBER =',I4, * ' NAME= ',A8/) 1004 FORMAT(' NBUG PARAMETERS'/12I5) 1005 FORMAT(' === NON-MODEL POTENTIAL CALCULATION ===') 1006 FORMAT(' === MODEL POTENTIAL CALCULATION ===') 1007 FORMAT(' === ANALYTIC BOUND ORBITALS BEING USED ===') 1008 FORMAT(' === NUMERICAL BOUND ORBITALS BEING USED ===') 1009 FORMAT(/' BASIC DATA'/) 1010 FORMAT(' NELC =',I3,' NZ =',I3,' LRANG1 =',I3,' NRANG2 =',I3, * ' LFIXN =',I3, ' LAMAX =',I2) 1011 FORMAT(' MAXNHF=',35I3/) 1012 FORMAT(' MAXNLG=',35I3/) 1013 FORMAT(/' DATA DEFINING THE STATIC POTENTIAL') 1014 FORMAT(' L =',I3/(' ',35I3)) 1015 FORMAT(//' THE RADIAL FUNCTIONS'/TR37,' SLATER-TYPE',TR5, *' POWER OF R',TR5,' EXPONENT') 1016 FORMAT(/TR13,I2,A,' ORBITAL') 1017 FORMAT(TR37,E12.5,TR9,I3,TR9,E12.5) 1018 FORMAT(//' R-MATRIX BOUNDARY CONDITIONS, RA =',F10.5,' BSTO = *',F10.5//' AMPLITUDE OF THE FUNCTIONS AT RA'/) 1019 FORMAT(//' INTEGRATION MESH'//' NIX=',I3) 1020 FORMAT(' HINT =',E14.7) 1021 FORMAT(' IHX=',20I5) 1022 FORMAT(' IRX=',20I5//) 1023 FORMAT(' MAXNC=',20I5) 1024 FORMAT(' LPOT=',I3) 1025 FORMAT(' LPOS=',20I5) 1026 FORMAT(I5,A,TR4,E14.7) C ETA=0.0D0 C C SET THE CHANNEL UNIT NUMBERS IREAD=7 IWRITE=6 JBUFF1=23 JBUFF2=24 ITAPE3=25 INX1=35+(iam+4000) C i1=iam/100 i2=(iam-100*(iam/100))/10 i3=iam-(100*(iam/100))-i2*10 C FF1='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FF2='RAD2.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FF3='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) C C OPEN THE DIRECT ACCESS FILES OPEN(UNIT=JBUFF1,FILE=FF1,ACCESS='DIRECT', 1 STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1) OPEN(UNIT=JBUFF2,FILE=FF2,ACCESS='DIRECT', 1 STATUS='UNKNOWN',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2) OPEN(ITAPE3,FILE=FF3,ACCESS='SEQUENTIAL', 1 STATUS='UNKNOWN',FORM='UNFORMATTED') C OPEN (INX1,FILE='NX1.DAT',ACCESS='SEQUENTIAL',STATUS='OLD', 1 FORM='UNFORMATTED') C OPEN (INX1,FILE='tapenx1',ACCESS='SEQUENTIAL',STATUS='OLD', C WRITE(IWRITE,1000) WRITE(IWRITE,1002) WRITE(IWRITE,1003)IWRITE,IREAD,ITAPE1,JBUFF1,FILE1, * JBUFF2,FILE2,ITAPE3,FILE3 WRITE (IWRITE,1004)NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, 1 NBUG8,NBUG9 C WRITE(IWRITE,1009) READ (INX1) LRANG1,NRANG2,LAMAX,NIX,NPTS,NRKPTS READ (INX1) NELC,NZ,LFIXN,IPSEUD C EFFECTIVE CHARGE IPSEUD=-IABS(IPSEUD) IF(MAXC.GT.1)NRANG2=MAXC IF (LRANG1 .GT. IDIM1) CALL RECOV1(1,LRANG1,IDIM1) IF (NRANG2 .GT. IDIM3) CALL RECOV1(3,NRANG2,IDIM3) IF (LAMAX .GT. IDIM4) CALL RECOV1(4,LAMAX,IDIM4) IF (MINIM .LE. LRANG1) THEN WRITE (IWRITE,'('' *** NOEXCHANGE PROGRAM NOT VALID'', / 1 '' *** BOUND TERMS ARE BEING IGNORED ***'')') END IF C C OPEN FILE TO STORE BOUND ORBITALS REQUIRED IN LAGRANGE C ORTHOGONALISATION C IF (LRANG1 .GE. MINIM) THEN OPEN(60,ACCESS='DIRECT',STATUS='SCRATCH', 1 FORM='UNFORMATTED',RECL=MACDIM*IDIM9*KDIM10) END IF C C CONSTRUCT NVARY(L) C IF (LFIXN .LT. LNDEF) THEN LFIXN=LNDEF IF(LFIXN0.GT.0)LFIXN=LFIXN0 WRITE(IWRITE,'('' LFIXN HAS BEEN RESET TO '',I3)')LFIXN END IF C IF(MAXC.GT.1.AND.LFIXN.LT.MINIM)LFIXN=MINIM IF(LFIXN0.GT.0)LFIXN=LFIXN0 C DO 1 I=1,LFIXN NVARY(I)=NRANG2 1 CONTINUE C DO 2 I=LFIXN+1,LRANG2 NVARY(I)=0 2 CONTINUE C WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX NUMERC=1 READ (INX1) (MAXNHF(I),I=1,LRANG1), (MAXNLG(I),I=1,LRANG1) WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1) WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1) NBOUND=0 DO 10,I=1,LRANG1 NHF=MAXNHF(I) NLG=MAXNLG(I) IF (IDIM9 .LE. NLG-I) CALL RECOV1(9,NLG-I+1,IDIM9) IF (IDIM10 .LT. NHF-NLG) CALL RECOV1(10,NHF-NLG,IDIM10) IF (IDIM12 .LT. NHF) CALL RECOV1(12,NHF,IDIM12) NBOUND=NBOUND+NHF-I+1 10 CONTINUE IF (IDIM13 .LT. NBOUND) CALL RECOV1(13,NBOUND,IDIM13) C READ (INX1) HINT,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX) WRITE(IWRITE,1019) NIX IF (IDIM6 .LT. NIX) CALL RECOV1(6,NIX,IDIM6) WRITE(IWRITE,1020) HINT WRITE(IWRITE,1021) (IHX(I),I=1,NIX) WRITE(IWRITE,1022) (IRX(I),I=1,NIX) IF (IDIM7 .LE. IRX(NIX))CALL RECOV1(7,IRX(NIX)+1,IDIM7) LPOT=1 DO 150 I=1,LRANG2 MAXNC(I)=I-1 LPOSX(I)=1 150 CONTINUE IF(IPSEUD.NE.0)THEN READ(INX1) (MAXNC(I),I=1,LRANG1) READ(INX1) LPOT,(LPOSX(I),I=1,LPOT) DO 155 I=1,LPOT READ(INX1)(POTHAM(I,J),J=1,IRX(NIX)) 155 CONTINUE ENDIF C READ (INX1) RA,BSTO,EK2MAX WRITE (IWRITE,1018) RA,BSTO IF(MAXC.GT.1)EK2MAX=ZERO WRITE(IWRITE,'(/'' ENERGY LIMITS FOR THE BUTTLE CORRECTION''/ 1 '' MAXIMUM ENERGY = '',F14.7)')EK2MAX C IF (MJS .NE. 0) THEN WRITE(IWRITE,'('' MJS FITTING PROCEDURE REQUESTED'')') IF (BSTO .NE. 0) THEN WRITE(IWRITE,'('' BUT QUADRATIC FIT WILL BE USED SINCE'', 1 '' BSTO IS NON-ZERO'')') MJS=0 END IF ELSE WRITE(IWRITE,'('' QUADRATIC FITTING REQUESTED'')') END IF C C DEFINE THE MAXPN ARRAY . NOT USED? DO 160,I=1,LRANG1 MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1 160 CONTINUE C EXTEND THE MAXNHF AND MAXNLG ARRAYS IF (LRANG2.GT.LRANG1) THEN DO 170,I=LRANG1+1,LRANG2 MAXNHF(I)=I-1 MAXNLG(I)=I-1 170 CONTINUE ENDIF C C READ IN BOUND ORBITALS C NQ=0 NPTS=IRX(NIX) DO 301 L=1,LRANG1 MAXHF=MAXNHF(L) DO 302 N=L,MAXHF NQ=NQ+1 IPOS(N,L)=NQ READ (INX1) (UJ(I,NQ),I=1,NPTS+1) 302 CONTINUE 301 CONTINUE C C READ IN POTENTIAL C READ (INX1) (POVALU(I),I=1,NPTS*2) C C SET UP ARRAYS FOR USE IN SIMPSONS RULE INTEGRATION ONE3 =ONE /THREE TWO3 =TWO *ONE3 FOUR3=FOUR*ONE3 IFI=1 XR(1) =ZERO STEP(1)=ZERO DO 200,I=1,NIX H =HINT*DBLE(IHX(I)) WT(IFI) =ONE3*(H+STEP(IFI)) RSTART=XR(IFI) IST =IFI+2 IFI =IRX(I)+1 IFX =0 IF (I.GT.1) IFX=IRX(I-1) DO 210,J=IST,IFI,2 XR(J-1) =RSTART+(J-2-IFX)*H XR(J ) =RSTART+(J-1-IFX)*H STEP(J-1)=H STEP(J )=H WT(J-1 )=FOUR3*H WT(J )=TWO3*H 210 CONTINUE 200 CONTINUE WT(IFI)=ONE3*H NPTS=IFI C C INTERPOLATE THE BOUND ORBITALS AND STORE ON SCRATCH FILE C READY FOR LAGRANGE ORTHOGONALISATION C A=TNINE/SIXT B=ONE /SIXT DO 401 L=1,LRANG1 MAXLG=MAXNLG(L) IF (L .GE. MINIM) THEN C C TAKE THE FACTOR R**L OUTSIDE FOR THE FIRST IFIN POINTS C IF (L .EQ. 1) THEN NIX1=1 IFIN=4 IRXP(1)=3 GO TO 403 END IF XRFIN=L*(L-1) XRFIN=XRFIN/(2.D0*NZ) IR1=IRX(1)+1 IF (XRFIN .LT. XR(IR1)) THEN NIX1=1 IFIN=NINT(XRFIN/HINT) IF (IFIN .LT. 4) IFIN=4 IRXP(1)=IFIN-1 GO TO 403 END IF IRXP(1)=IRX(1) DO 402 NI=2,NIX IR2=IRX(NI)+1 IF (XRFIN .LE. XR(IR2)) THEN NIX1=NI IFIN=NINT((XRFIN-XR(IR1))/(HINT*IHX(NI))) IF (IFIN .LT. 3) IFIN=3 IFIN=IFIN+IR1 IRXP(NI)=IFIN-1 GO TO 403 END IF IRXP(NI)=IRX(NI) IR1=IR2 402 CONTINUE 403 KT=0 DO 404 N=L,MAXLG NQ=IPOS(N,L) KT=KT+1 DO 405 I=2,IFIN UTEMP(I)=UJ(I,NQ)/(XR(I)**L) 405 CONTINUE H=HALF*HINT PTEMP(1)=HALF*(3*UTEMP(2)-UTEMP(3))*H**L PTEMP(2)=(P75*UTEMP(3)+P125*(3*UTEMP(2)-UTEMP(4))) 1 *(H+HINT)**L IX=3 DO 406 I=1,NIX1 H=HALF*HINT*IHX(I) DO 407 J=IX,IRXP(I)-1 PTEMP(J)=(A*(UTEMP(J)+UTEMP(J+1))-B*( 1 UTEMP(J-1)+UTEMP(J+2))) 2 *(XR(J)+H)**L 407 CONTINUE IF (I.NE.NIX1) THEN C CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED K2=IRX(I) K1=K2+1 PTEMP(K1)=(A*(UTEMP(K1)+UTEMP(K1+1))-B*( 1 UTEMP(K1-2)+UTEMP(K1+2))) PTEMP(K2)=(A*(UTEMP(K2)+UTEMP(K2+1))-B*( 1 PTEMP(K1)+UTEMP(K2-1))) 2 *(XR(K2)+H)**L PTEMP(K1)=PTEMP(K1)*(XR(K1)+2*H)**L IX=K1+1 ELSE IX=IRXP(NIX1) END IF 406 CONTINUE DO 408 I=IX-1,IRX(NIX)+1 UTEMP(I)=UJ(I,NQ) 408 CONTINUE DO 409 I=NIX1,NIX H=HALF*HINT*IHX(I) DO 410 J=IX,IRX(I)-1 PTEMP(J)=A*(UTEMP(J)+UTEMP(J+1))-B*( 1 UTEMP(J-1)+UTEMP(J+2)) 410 CONTINUE IF (I.NE.NIX) THEN C CHANGE OF INTERVAL...DOUBLING OF STEP SIZE IS ASSUMED K2=IRX(I) K1=K2+1 PTEMP(K1)=A*(UTEMP(K1)+UTEMP(K1+1))-B*( 1 UTEMP(K1-2)+UTEMP(K1+2)) PTEMP(K2)=A*(UTEMP(K2)+UTEMP(K2+1))-B*( 1 PTEMP(K1)+UTEMP(K2-1)) IX=K1+1 ELSE K1=IRX(NIX) PTEMP(K1)=P75*UTEMP(K1) + 1 P125*(3*UTEMP(K1+1)-UTEMP(K1-1)) END IF 409 CONTINUE DO 411 I=1,NPTS-1 P(KT,I+I-1)=PTEMP(I) P(KT,I+I)=UJ(I+1,NQ) 411 CONTINUE 404 CONTINUE C IF (KT .GT. 0) WRITE(60,REC=L)P C END IF 401 CONTINUE REWIND (INX1) RETURN END C*********************************************************************** SUBROUTINE WRITAP use param use comm_interface, only : iam,nproc C C WRITES BASIC INFORMATION ONTO THE STG1R PERMANENT OUTPUT FILE C IMPLICIT REAL*8 (A-H,O-Z) C c INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2, * IDIM11=MZSLT,IDIM12=MZNR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM12=IDIM1*IDIM12) C COMMON /BASIC / BSTO,ETOT,RA,W1,WINT,WMAX,WMIN, * LRGL,NAST,NELC,NPTY,NRANG2,NVARY(IDIM2),NSPN COMMON /BASIN / EIGENS(IDIM3,IDIM2),ENDS(IDIM3,IDIM2),DELTA,ETA COMMON /BUTT/COEFF(3,IDIM2),EK2MAX,EK2MIN COMMON /DISTAP/ IDISC1,IDISC2,IDISC3,ITAPE1,ITAPE2,ITAPE3,ITAPE4 COMMON /RADIAL/ C(IDIM1,IDIM12,IDIM11),ZE(IDIM1,IDIM12,IDIM11), * IRAD(IDIM1,IDIM12,IDIM11),NCO(KDIM12), * LRANG1,LRANG2,MAXNHF(IDIM2),MAXNLG(IDIM2), * NCOEFF,NLIMIT,NZ,ISTAT(IDIM1,IDIM12),LLPDIM C C IF(ITOTAL .LE. 0) RETURN ITAPE=ITAPE3 REWIND ITAPE WRITE(ITAPE)RA,BSTO WRITE(ITAPE)(NVARY(L),L=1,LRANG2) WRITE(ITAPE)((ENDS(N,L),N=1,NVARY(L)),L=1,LRANG2) WRITE(ITAPE)((COEFF(I,L),I=1,3),L=1,LRANG2) C RETURN END