c PSTG3NX:: Parallel Version 1.5 of stg3nx 19/08/09 c c C. P. Ballance, D. M. Mitnik, N. R. Badnell, and D. C. Griffin c c V1.1 Initial release c V1.2 Use of pdsyevd c V1.3 Re-work gathering of surface amplitudes c V1.4 Multiple reads instead of read/broadcast c V1.5 Correct read of Buttle from pstg2nx c c C N. R. BADNELL UoS v1.11 - QUB vx.x 29/04/05 C C PROGRAM NXHAM/STG3NX C C *** THIS IS THE HAMILTONIAN STAGE OF RMATRX NX C C VERSION WHICH ALLOWS NRANG2 TO DIMINISH WITH INCREASE C IN LARGE L USING THE ARRAY NVARY C PROGRAM MAIN IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C CHARACTER*1 NUM(0:9) CHARACTER*11 ffff LOGICAL EX,exs C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO) PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM6=(IDIM4+1)/2) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ICDIM3=ICDIM2+ICDIM2-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGBUF/BUFV1(LBUFFV),BUFV2(LBUFFV), 1 IBUFV1(LBUFFV),IBUFV2(LBUFFV) COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4) COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2), 1 ICONCH(ITDIM3),KCONCH(ITDIM3), 2 L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1), 3 NSCHAN(IPDIM1),NSPREC(IPDIM1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/HPOINT/KRECH(0:ILDIM4) COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RADBUF/BUFF1(LBUFF1),BUFF2(LBUFF1) COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON/SETL /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) COMMON/STG1 /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2) C COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DATA IFLAGMAIN/0/ SAVE IFLAGMAIN C 1000 FORMAT(TR31,' END OF NXHAM'/TR31,' ------------'///) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc common/procdat/ictxt,myrow,mycol common/crashblock/istop real*8 HMAT allocatable :: HMAT(:,:) call mpi_init(ierr) call mpi_comm_rank(mpi_comm_world,iam,ierr) call mpi_comm_size(mpi_comm_world,nproc,ierr) istop=0 call cpu_time(timei) c **** parallel **** c write(0,*) ' begin proc=:',iam c 8448 CONTINUE C if(iam.eq.0)then CALL HEADER endif C C CALL RECSET(LRGLLO,LRGLUP,ICONTN) C if(iam.eq.0)then CALL BUTTLECOEFF endif C IF (ICONTN .EQ. 1) THEN CALL STGRRX ELSE CALL STGRRD END IF C c **** parallel **** inquire(file='RAD1.DAT',exist=exs) c........read input file for parallel diagonalization if(iflagmain.eq.0)then call READ3P if (istop.ne.0) go to 9999 c.......initialize the BLACS context call INITCTXT c.......bails out if this process is not a part of this context if (myrow.eq.-1) go to 9999 c **** parallel **** endif c CALL STGHRD(ICONTN) C C READ STG1 TAPE C CALL RETAP1 C CALL CRLMAT C C WRITE ASYMPTOTIC FILE C CALL WRITE1 C DO 10 LRGL=LRGLLO,LRGLUP DO 11 NPTY0=NPTYMN,NPTYMX IF(NPTYIN.LT.0)THEN NPTY=NPTY0 ELSE NPTY=MOD(LRGL,2)+NPTY0 NPTY=MOD(NPTY,2) ENDIF CALL CUSETL CALL SETCHN(LRGLLO) CALL STHMAT(LRGLLO) DO 6 ISP=1,NSP C IF (LRGL .EQ. LRGLUP .AND.NPTY0.EQ.NPTYMX.AND. ISP .EQ. NSP) 1 THEN if(.not.exs)then IFLAGMAIN=IFLAGMAIN+1 i1=iflagmain/100 i2=(iflagmain-100*(iflagmain/100))/10 i3=iflagmain-(100*(iflagmain/100))-i2*10 ffff='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) IFLAGMAIN=IFLAGMAIN-1 inquire (file=ffff,exist=ex) if(.not.ex) then MORE=0 else MORE=2 endif else MORE=0 endif ELSE MORE= 2 END IF C c **** parallel **** c.......calculate the size of the local Hamiltonian igbsize = lenham(isp) call HSIZE(igbsize,irowsize,icolsize) c write(0,*) ' NN=:',igbsize,' local=:',irowsize,'x',icolsize, c + ' proc=:',iam c.......allocate the local Hamiltonian allocate(HMAT(irowsize,icolsize),stat=ierr) if (ierr.ne.0) then print*,'program could not allocate space for HMAT' print*,' irowsize=:',irowsize,' icolsize=:',icolsize stop endif c.......initializate HMAT do icol=1,icolsize do irow=1,irowsize HMAT(irow,icol)=0.0D0 enddo enddo c **** parallel **** C C READ THE ASYMPTOTIC COEFFICIENTS AND THE HAMILTONIAN MATRIX C FROM DIRECT ACCESS FILE JBUFD. C CALCULATE ENERGY LEVELS AND SURFACE AMPLITUDES. C WRITE THE ASYMPTOTIC FILE BLOCKS FOR THIS LRGL,NPTY C AND TARGET SPIN C CALL RDHMAT(ISP,MORE,HMAT,irowsize,icolsize) c **** parallel **** CALL DIAGHMAT(HMAT,irowsize,icolsize,ISP) c **** parallel **** deallocate(HMAT) !! parallel 6 CONTINUE 11 CONTINUE 10 CONTINUE C CLOSE (JBUFIR) CLOSE (JBUFFV) CLOSE (JBUFFZ) CLOSE (JBUFF1) CLOSE (JBUFF2) CLOSE (ITAPE1) C c **** parallel **** c call mpi_barrier(mpi_comm_world,ierr) c CLOSE (JBUFD,status='delete') !not if scratch if(.not.exs)then IFLAGMAIN=IFLAGMAIN+1 i1=iflagmain/100 i2=(iflagmain-100*(iflagmain/100))/10 i3=iflagmain-(100*(iflagmain/100))-i2*10 ffff='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) inquire (file=ffff,exist=ex) if(ex) goto 8448 endif c **** parallel **** C C if (iam.eq.0) then c CLOSE (JBUFD) !! not if parallel c ENDFILE ITAPE3 c REWIND ITAPE3 CLOSE (ITAPE3) endif C if (iam.eq.0)WRITE(IWRITE,1000) C c **** parallel **** if (iam.eq.0) then call cpu_time(timef) time=timef-timei TIME=TIME/60.0 write(33,2999) time 2999 format(//2x,' Total CPU time=:',1pg10.3,' min.') WRITE(IWRITE,999) TIME,nproc 999 FORMAT(//1X,'CPU TIME=',F9.3,' MIN -- processors:',i4) endif call mpi_barrier(mpi_comm_world,ierr) call blacs_gridexit(ictxt) c call blacs_exit(0) 9999 call mpi_finalize(ierr) c **** parallel **** C C STOP END ********************************************************************** BLOCK DATA IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS C C SET GLOBAL REAL CONSTANTS C DATA ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS/ 1 0.0D 00,0.1D 00,0.5D 00,1.0D 00,2.0D 00,3.0D 00,4.0D 00, 2 7.0D 00,1.1D 01,1.0D-09/ C END C*********************************************************************** SUBROUTINE CMATII(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT) C C TO CALCULATE THE CONTRIBUTION TO THE HAMILTONIAN FROM THE C INTERACTION OF A PAIR OF SIMILAR CONFIGURATIONS C THIS SUBMATRIX IS STORED IN CMAT C ZB IS THE ARRAY OF Z-COEFFICIENTS FOR CURRENT L,LP FOR C LAMBDA=LAMST TO LAMFIN C ICOUNT IS THE COUNTER TO LOCATE VSH FOR THIS CONFIG PAIR C AT EXIT ICOUNT HOLDS THE COUNT FOR THE NEXT CONFIGURATION PAIR C LAMCFG IS THE MAX. LAMBDA FOR WHICH THERE IS A VSH VALUE C FOR THIS CONFIG PAIR C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2) C DIMENSION ZB(ILDIM1) DIMENSION CMAT(IDIM3,IDIM3),VRMAT(IDIM3,IDIM3),RM(INDIM2) c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** C 1000 FORMAT(/,28X,'SUBROUTINE CMATII'/28X,'-----------------'/) 1001 FORMAT(' LAMBDA VRMAT') 1002 FORMAT(' LAMBDA CMAT') 1003 FORMAT(4X,I6,/2X,10E11.4,/2X,10E11.4) C IF (NBUG1 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1000) END IF C C ZEROISE THE MATRIX CMAT(N,NP) READY TO ACCUMULATE THE ANGULAR TERMS C DO 1 N=1,IDIM3 DO 2 NP=1,IDIM3 CMAT(NP,N)=ZERO 2 CONTINUE 1 CONTINUE C C SET LIMITS FOR N AND NP C NRA1=NVARY(L+1) NRA2=NVARY(LP+1) C C SET LAMBDA LIMITS CORRESPONDING TO VSHELL FILE C LAMLO=MOD (ABS (IRHSGL(ICOUNT)), 100) LAMUP=LAMCFG C DO 3 LAM=LAMLO,LAMUP,2 NSHSGL=ABS( IRHSGL(ICOUNT)) C C FIND THE NUMBER OF SHELLS CONTRIBUTING FOR THIS LAMBDA C NSH=NSHSGL/10000 IF(NSH .EQ. 0) THEN C C IF NO SHELLS CONTRIBUTE HIGHER VALUES OF LAMBDA ALSO C CANNOT CONTRIBUTE SO RETURN C ICOUNT=ICOUNT+1 RETURN END IF C IF(LAM .LT. LAMST .OR. LAM .GT. LAMFIN) THEN C C SKIP TO NEXT LAMBDA LOCATION OVER THE NSH VALUES C ICOUNT=ICOUNT+NSH GO TO 3 END IF C C ZEROISE VRMAT(N,NP) TO ACCUMULATE CONTRIBUTION FROM SHELLS C FOR THIS LAMBDA VALUE C DO 4 N=1,IDIM3 DO 5 NP=1,IDIM3 VRMAT(NP,N)=ZERO 5 CONTINUE 4 CONTINUE NSHHUN=NSH*100 C C SUM OVER SHELLS C DO 6 ISH=1,NSH ISIG=ABS( IRHSGL(ICOUNT))/100 - NSHHUN VS=VSH(ICOUNT) ICOUNT=ICOUNT+1 IF(ABS(VS) .LT. EPS) GO TO 6 LSIG=LJCOMP(ISIG) NSIG=NJCOMP(ISIG) C C LOCATE THE RADIAL INTEGRAL C CALL LOCCCD(LSIG+1, LSIG+1, NSIG, NSIG, LAM+1, RM) IF (NBUG3 .EQ. 1) THEN if (iam.eq.0) x WRITE(IWRITE,'(/''RK INTEGRAL FOR L='',I4,''LP='',I4, 1 ''LAM='',I4,''L1='',I4,''L1P='',I4,''N1='',I4,''N1P='', 2 I4)') L,LP,LAM,LSIG,LSIG,NSIG,NSIG if (iam.eq.0)WRITE(IWRITE,'(2X,8E14.7)') (RM(I),I=1,ICCLNG) END IF C C I=0 DO 7 N=1,NRA1 IF(L .EQ. LP) THEN DO 71 NP=1,N I=I+1 VRMAT(N,NP)=VRMAT(N,NP)+RM(I)*VS 71 CONTINUE C ELSE C DO 72 NP=1,NRA2 I=I+1 VRMAT(NP,N)=VRMAT(NP,N)+RM(I)*VS 72 CONTINUE C END IF C 7 CONTINUE C 6 CONTINUE IF (NBUG1 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1001) if (iam.eq.0)WRITE(IWRITE,1003)LAM, + ((VRMAT(NP,N),NP=1,NRA2),N=1,NRA1) END IF C C MOVE THE ZBUFF COUNTER TO LOCATE THE Z COEFFICIENT C FOR THIS LAMBDA VALUE C IZCONT=1+(LAM-LAMST)/2 Z=ZB(IZCONT) C C ADD Z*VR INTO CMAT(N,NP) C DO 9 N=1,NRA1 IF(L .EQ. LP) THEN NPLO=N ELSE NPLO=1 END IF C DO 10 NP=NPLO,NRA2 CMAT(NP,N)=CMAT(NP,N)+Z*VRMAT(NP,N) 10 CONTINUE C 9 CONTINUE IF (NBUG1 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1002) if (iam.eq.0)WRITE(IWRITE,1003)LAM, + ((CMAT(NP,N),NP=1,NRA2),N=1,NRA1) END IF C 3 CONTINUE C RETURN END ********************************************************************** SUBROUTINE CMATIJ( 1 L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT,IND) C C TO CALCULATE THE CONTRIBUTION TO THE HAMILTONIAN MATRIX C FROM THE INTERACTION OF A PAIR OF CONFIGURATIONS WHICH ARE C NOT IDENTICAL C THIS SUB MATRIX IS STORED IN CMAT C ZB IS THE ARRAY OF Z-COEFFICIENTS FOR CURRENT L,LP FOR C LAMCFG IS THE MAX. LAMBDA ALLOWED BY TARGET CONFIGURATION C ANGULAR MOMENTUM C ICOUNT IS THE COUNTER TO LOCATE VSHELL FOR THIS CONFIG PAIR C AT EXIT ICOUNT HOLDS THE COUNT FOR THE NEXT CONFIGURATION PAIR C IND RETURNS 1 IF IC,JC COUPLED C 0 IF THEY DID NOT C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** C IND INITIALLY SET ZERO C DIMENSION ZB(ILDIM1) DIMENSION CMAT(IDIM3,IDIM3),RM(INDIM2) C 1000 FORMAT(/,28X,'SUBROUTINE CMATIJ'/28X,'-----------------'/) 1002 FORMAT(' LAMBDA CMAT') 1003 FORMAT(4X,I6,/2X,10E11.4,/2X,10E11.4) C IF (NBUG1 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1000) END IF C IND=0 C IRSL=IRHSGL(ICOUNT) IF(IRSL .EQ. 9999) THEN ICOUNT=ICOUNT+1 RETURN END IF C C FIND THE INTERACTING SHELLS AND THE FIRST LAMBDA FROM IRSL C LAMLO=MOD(IRSL,100) IRS=IRSL/100 IRHO=IRS/100 ISIG=MOD(IRS,100) LRHO=LJCOMP(IRHO) LSIG=LJCOMP(ISIG) LAMUP=MIN(LRHO+LSIG,LAMCFG) C C TEST WHETHER THESE LIMITS OVERLAP THE LIMITS ON LAMBDA C FOR THE Z COEFFICIENTS C IF(LAMLO .GT. LAMFIN .OR. 1 LAMUP .LT. LAMST) THEN ICOUNT=ICOUNT+1+(LAMUP-LAMLO)/2 RETURN END IF C NRHO=NJCOMP(IRHO) NSIG=NJCOMP(ISIG) C C ZEROISE THE CMAT MATRIX READY TO SUM OVER LAMBDA C DO 1 N=1,IDIM3 DO 2 NP=1,IDIM3 CMAT(NP,N)=ZERO 2 CONTINUE 1 CONTINUE C C SET LIMITS ON N AND NP C NRA1=NVARY(L+1) NRA2=NVARY(LP+1) C C SET IND C IND=1 C DO 3 LAM=LAMLO,LAMUP,2 IF (LAM .LT. LAMST .OR. LAM .GT. LAMFIN) THEN C C MOVE ALONG VSHELL FILE C ICOUNT=ICOUNT+1 GO TO 3 END IF C VS=VSH(ICOUNT) ICOUNT=ICOUNT+1 IF(ABS(VS) .LE. EPS) THEN GO TO 3 END IF C C MOVE ZBUFF COUNTER TO LOCATE Z FOR THIS LAMBDA C IZCONT=1+(LAM-LAMST)/2 VZ=ZB(IZCONT)*VS C C LOCATE THE RADIAL INTEGRAL C IF (LRHO .LT. LSIG) THEN CALL LOCCCD(LRHO+1, LSIG+1, NRHO, NSIG, LAM+1, RM) ELSE CALL LOCCCD(LSIG+1, LRHO+1, NSIG, NRHO, LAM+1, RM) END IF C C RM(I) CONTAINS THE NRA1 X NRA2 MATRIX OF RADIAL C INTEGRALS. C WHEN L=LP ONLY THE LOWER HALF OF THE MATRIX IS STORED C RM(I) IS NOW MULTIPLIED BY THE FACTOR VZ AND STORED IN CMAT(N,NP) C WHEN L=LP IT IS STORED IN THE UPPER HALF OF CMAT C IF (NBUG3 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,'(/''RK INTEGRAL FOR L='',I4,''LP='',I4, 1 ''LAM='',I4,''L1='',I4,''L1P='',I4,''N1='',I4,''N1P='', 2 I4)') L,LP,LAM,LRHO,LSIG,NRHO,NSIG if (iam.eq.0)WRITE(IWRITE,'(2X,8E14.7)') (RM(I),I=1,ICCLNG) END IF C I=0 DO 4 N=1,NRA1 IF(L .EQ. LP) THEN NPUP=N ELSE NPUP=NRA2 END IF C DO 5 NP=1,NPUP I=I+1 IF(L .EQ. LP) THEN CMAT(N,NP)=CMAT(N,NP)+RM(I)*VZ ELSE CMAT(NP,N)=CMAT(NP,N)+RM(I)*VZ END IF C 5 CONTINUE C 4 CONTINUE C IF (NBUG1 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1002) if (iam.eq.0)WRITE(IWRITE,1003)LAM, + ((CMAT(NP,N),NP=1,NRA2),N=1,NRA1) END IF 3 CONTINUE C RETURN END ********************************************************************** SUBROUTINE CRLMAT C C CALCULATES THE REDUCED MATRIX WHOSE C ELEMENTS (MULTIPLIED BY A Z COEFFICIENT) C WILL BE USED TO CALCULATE THE ASYMPTOTIC C COEFFICIENTS C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM6=(IDIM4+1)/2) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2) C DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6), CRLADD(IDIM4) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 1000 FORMAT(/,28X,'SUBROUTINE CRLMAT'/28X,'-----------------', 1 //' IS JS ICS JCS LAM CRLADD') 1001 FORMAT(' ',5I5,E14.7) 1002 FORMAT(' CRL MATRIX LAMBDA=',I2) 1003 FORMAT(' ',8E14.7) C 3000 FORMAT(' **** STOP IN CRLMAT ****'/ 1 ' ERROR IN USING VSH FILE FOR SETS ', 2I4) 3001 FORMAT(' **** STOP IN CRLMAT ****'/ 1 ' WRONG VSH FILE LOADED') C IF (NBUG2 .GT. 1) THEN if (iam.eq.0)WRITE(IWRITE,1000) END IF LAMIND=(LAMAX+1)/2 C C READ FIRST RADIAL INTEGRAL POINTER BLOCK C READ (ITAPE1) (((IBBPOL(I,J,K),I=1,LRANG1),J=1,LRANG1), 1 K=1,LAMIND),NBBPOL,(IST3(I),I=1,LRANG2) C C ZEROISE CRL MATRIX C DO 100 IST=1,NAST DO 101 JST=1,NAST DO 102 LM=1,LAMAX CRL(IST,JST,LM)=ZERO 102 CONTINUE 101 CONTINUE 100 CONTINUE C NRECR=0 NB=0 C C READ BOUND BOUND MULTIPOLE INTEGRALS INTO RKCCD ARRAY C THIS SAVES SPACE C THESE ARE THE FIRST INTEGRALS STORED SO THEY START AT C THE BEGINNING OF THE SECOND RECORD C ISTART=LBUFF1*2 + 1 IFIN=ISTART + NBBPOL - 1 CALL RDRINT(ISTART,IFIN,NB,NRECR) C C READ ANGULAR INTEGRAL POINTERS AND SET C RECORD NUMBERS ZERO READY TO START READING C ANGULAR INTEGRALS C READ (JBUFIR, REC=2) NS,((LAMIJ(IS,JS),IS=1,NS),JS=1,NS), 1 (IVCONT(ISJS),ISJS=1,(NS*(NS+1))/2 + 1) IF (NS .NE. NSET) THEN if (iam.eq.0) WRITE (IWRITE,3001) STOP END IF NRECV1=0 NRECV2=0 C C LOOP OVER THE SETS C DO 1 IS=1,NSET C DO 2 JS=IS,NSET C LAMLU=LAMIJ(IS,JS) C C TEST WHETHER SETS IS,JS CAN COUPLE C IF THEY CAN LAMIJ CONTAINS THE LIMITS ON LAMBDA IMPOSED C BY SHELL COUPLING BETWEEN CONFIGURATIONS IN THE SETS C IF(LAMLU .EQ. -999) GO TO 2 IF(LAMLU .EQ. 0) GO TO 2 C C LAMLO=LAMLU/100 LAMUP=MOD(LAMLU,100) C C TRANSFER TARGET ANGULAR INTEGRALS FOR THIS SET C COMBINATION TO VSH C CALL RDVSH(IS,JS,NRECV1,NRECV2,NVSHEL) ICOUNT=1 C C LOOP OVER CONFIGURATIONS IN SET IS C DO 3 ICS=1,NCSET(IS) IF(IS .EQ. JS) THEN JCSLO=ICS ELSE JCSLO=1 END IF C C LOOP OVER CONFIGURATIONS IN SET JS C DO 4 JCS=JCSLO,NCSET(JS) C C THE CASE OF CONFIGURATIONS HAVING THE SAME SHELL STRUCTURE C IS TREATED SEPARATELY C IF (IRHSGL(ICOUNT) .LT. 0) THEN C C ZEROISE THE ARRAY CRLADD(LAM) READY TO ACCUMULATE TERMS C DO 40 LAM=1,LAMAX CRLADD(LAM)=ZERO 40 CONTINUE C C SET LAMBDA LIMITS CORRESPONDING TO VSHELL FILE C LAM1=MOD (ABS (IRHSGL(ICOUNT)), 100) LAM2=LAMUP C DO 41 LAM=LAM1,LAM2,2 NSHSGL=ABS( IRHSGL(ICOUNT)) C C FIND THE NUMBER OF SHELLS CONTRIBUTING FOR THIS LAMBDA C NSH=NSHSGL/10000 IF(NSH .EQ. 0) THEN C C IF NO SHELLS CONTRIBUTE HIGHER VALUES OF LAMBDA ALSO C CANNOT CONTRIBUTE SO JUMP OUT OF LOOP C LAM2=LAM-2 ICOUNT=ICOUNT+1 GO TO 44 END IF C IF(LAM .EQ. 0 .OR. LAM .GT. LAMAX) THEN C C SKIP TO NEXT LAMBDA LOCATION OVER THE NSH VALUES C ICOUNT=ICOUNT+NSH GO TO 41 END IF NSHHUN=NSH*100 C C SUM OVER SHELLS C DO 42 ISH=1,NSH ISIG=ABS( IRHSGL(ICOUNT))/100 - NSHHUN VS=VSH(ICOUNT) ICOUNT=ICOUNT+1 IF(ABS(VS) .LT. EPS) GO TO 42 LSIG=LJCOMP(ISIG) NSIG=NJCOMP(ISIG) C C LOCATE THE RADIAL INTEGRAL C CALL LOCMBB(IBBPOL, LSIG+1, LSIG+1, NSIG, NSIG, LAM, RVAL) C CRLADD(LAM)=CRLADD(LAM) + VS*RVAL IF (NBUG2 .GT. 1) THEN if (iam.eq.0)WRITE(IWRITE,1001)IS,JS,ICS,JCS,LAM,CRLADD(LAM) END IF C 42 CONTINUE 41 CONTINUE C ELSE C C NOW TREAT CASE OF CONFIGURATIONS HAVING DIFFERENT C SHELL STRUCTURE C IRSL=IRHSGL(ICOUNT) IF(IRSL .EQ. 9999) THEN ICOUNT=ICOUNT+1 GO TO 4 END IF C C FIND THE INTERACTING SHELLS AND THE FIRST LAMBDA FROM IRSL C LAM1=MOD(IRSL,100) IRS=IRSL/100 IRHO=IRS/100 ISIG=MOD(IRS,100) LRHO=LJCOMP(IRHO) LSIG=LJCOMP(ISIG) LAM2=MIN(LRHO+LSIG,LAMUP) NRHO=NJCOMP(IRHO) NSIG=NJCOMP(ISIG) C DO 43 LAM=LAM1,LAM2,2 C IF (LAM .EQ. 0 .OR. LAM .GT. LAMAX) THEN ICOUNT=ICOUNT+1 GO TO 43 END IF C VS=VSH(ICOUNT) ICOUNT=ICOUNT+1 IF(ABS(VS) .LE. EPS) THEN CRLADD(LAM)=ZERO GO TO 43 END IF C C LOCATE THE RADIAL INTEGRAL C CALL LOCMBB(IBBPOL, LRHO+1, LSIG+1, NRHO, NSIG, LAM, RVAL) C CRLADD(LAM)=VS*RVAL IF (NBUG2 .GT. 1) THEN if (iam.eq.0)WRITE(IWRITE,1001)IS,JS,ICS,JCS,LAM,CRLADD(LAM) END IF C 43 CONTINUE C END IF C C LOOP OVER ALL TARGET STATES COMPOSED FROM SET IS C 44 DO 5 IST=1,NSTAT(IS) C C FIND THE STATE NUMBER AND THUS FIND THE COEFFICIENT AIJ C ISTAT=NSETST(IS,IST) AI=AIJ(ISTAT,ICS) C C WHEN SETS IS, JS ARE THE SAME A PAIR OF CONFIGURATIONS C WILL CONTRIBUTE TWICE TO COMPENSATE FOR THE SAVING MADE C AS A RESULT OF SYMMETRY. JCSLO=KICS IN DO LOOP 7 C IF(IS .EQ. JS .AND. ICS .NE. JCS) THEN AI2=AIJ(ISTAT,JCS) END IF C C MODIFY LIMITS FOR LOOPING OVER STATES FROM SET JS WHEN IS=JS C IF(IS .EQ. JS) THEN JSTLO=IST ELSE JSTLO=1 END IF C C LOOP OVER STATES COMPOSED FROM THE SET JS C DO 6 JST=JSTLO,NSTAT(JS) JSTAT=NSETST(JS,JST) AJ=AIJ(JSTAT,JCS) AIAJ=AI*AJ IF(IS .EQ. JS .AND. ICS .NE. JCS) THEN AIAJ=AIAJ+AI2*AIJ(JSTAT,ICS) END IF C IF(ABS(AIAJ) .LE. EPS) GO TO 6 C C ADD IN CONTRIBUTION TO CRL MATRIX FOR EACH LAMBDA C IF (LAM1 .EQ. 0) LAM1=LAM1+2 LAM2=MIN(LAM2,LAMAX) C DO 61 LAM=LAM1,LAM2,2 IF (JSTAT .LT. ISTAT) THEN CRL(JSTAT,ISTAT,LAM)=CRLADD(LAM)*AIAJ + CRL(JSTAT,ISTAT, 1 LAM) ELSE CRL(ISTAT,JSTAT,LAM)=CRLADD(LAM)*AIAJ + CRL(ISTAT,JSTAT, 1 LAM) END IF 61 CONTINUE C 6 CONTINUE C 5 CONTINUE C 4 CONTINUE C 3 CONTINUE C C TEST THAT THE VSH ARRAY HAS BEEN COMPLETELY USED FOR THIS C SET COMBINATION C IF (ICOUNT .NE. NVSHEL + 1) THEN if (iam.eq.0) WRITE (IWRITE,3000) IS,JS STOP END IF C 2 CONTINUE C 1 CONTINUE C DO 7 LAM=1,LAMAX IF (NBUG2 .GT. 0) THEN if (iam.eq.0) WRITE (IWRITE,1002) LAM if (iam.eq.0) WRITE (IWRITE,1003) ! ((CRL(ISTAT,JSTAT,LAM),ISTAT=1,NAST),JSTAT=1,NAST) END IF 7 CONTINUE RETURN END ********************************************************************** SUBROUTINE CUSETL C C TO DETERMINE WHICH SETS COUPLE WITH EACH VALUE OF THE C CONTINUUM ANGULAR MOMENTUM WITHIN THE RANGE ALLOWED C BY LRGL. THE INFORMATION IS STORED IN /SETL/. C THE COMMON BLOCK /SETL/ WILL BE STORED ON DISC FOR C EACH LRGL,NPTY COMBINATION. IT IS INDEPENDENT OF THE C TOTAL SPIN IN THE DIRECT CASE. C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ILDIM1=MZLR3) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) C COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/SETL/LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/SETS/NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) C C ARRAYS INTERMEDIATE TO /SETL/ ARRAYS C DIMENSION ILVAL(ILDIM2),INSCOL(ILDIM2),INSETL(ISDIM1,ILDIM2) C 1000 FORMAT(////' TOTAL ANGULAR MOMENTUM ',I3,' PARITY ',I3/ 1 ' -------------------------------------') 1001 FORMAT(' KL L SETS COUPLED TO L') 1002 FORMAT(2X,15I6) 1003 FORMAT(/' CONTINUUM ANGULAR MOMENTUM ', 1 'TAKES VALUES BETWEEN ',I3,' AND ',I3) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** if (iam.eq.0)WRITE(IWRITE,1000) LRGL,NPTY IF (NBUG8 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1001) END IF 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 1 IS=1,NSET LCFG=LSET(IS) LMIN=MIN(LMIN,ABS(LCFG-LRGL)) LMAX=MAX(LMAX, LCFG) 1 CONTINUE LMAX=LMAX+LRGL ILTOT=LMAX-LMIN+1 C DO 2 KL=1,ILTOT ILVAL(KL)=LMIN+KL-1 2 CONTINUE C C LOOP OVER THIS RANGE OF CONTINUUM L C DO 3 KL=1,ILTOT L=ILVAL(KL) NL=0 C DO 4 IS=1,NSET LCFG=LSET(IS) IPI=IPISET(IS) C C PARITY TEST C IF(MOD(IPI+L+NPTY,2) .NE. 0) GO TO 4 C C TRIANGLE RULE C IF(LRGL .GT. LCFG+L .OR. 1 LRGL .LT. ABS(LCFG-L)) GO TO 4 C NL=NL+1 INSETL(NL,KL)=IS 4 CONTINUE C INSCOL(KL)=NL 3 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 5 IKL=1,ILTOT IF (INSCOL(IKL) .NE. 0) THEN KL=KL+1 NSCOL(KL)=INSCOL(IKL) LVAL(KL)=ILVAL(IKL) DO 51 NL=1,NSCOL(KL) NSETL(NL,KL)=INSETL(NL,IKL) 51 CONTINUE END IF 5 CONTINUE LTOT=KL IF (LTOT .EQ. 0) GO TO 7 LMIN=LVAL(1) LMAX=LVAL(LTOT) IF (NBUG8 .EQ. 1) THEN DO 6 KL=1,LTOT if (iam.eq.0)WRITE(IWRITE,1002)KL,LVAL(KL), + (NSETL(NL,KL),NL=1,NSCOL(KL)) 6 CONTINUE END IF if (iam.eq.0) WRITE (IWRITE,1003) LMIN,LMAX C 7 RETURN END ********************************************************************** SUBROUTINE LOC1CC(L,LASTL,NRECR1,R1) C C LOCATE THE ONE ELECTRON CONTINUUM CONTINUUM RADIAL INTEGRAL C L IS THE CONTINUUM ANGULAR MOMENTUM + 1 C NRECR1 IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE C R1 IS A SYMMETRIC MATRIX STORED AS A 1-DIMENSIONAL ARRAY C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG C DIMENSION R(IDIM3,IDIM3),R1(INDIM2) C IST3L=IST3(L) NBUFF1=(IST3L-1)/LBUFF1 !IST3L -> IST3L-1 IF (NBUFF1 .NE. NRECR1) THEN READ (JBUFF1, REC=NBUFF1)R1CC END IF C IBUFF1=MOD(IST3L-1,LBUFF1) !IST3L ->IST3L-1, )-1 -> ) C NRA=NVARY(L) DO 1 N=1,NRA C DO 2 NP=1,N IBUFF1=IBUFF1+1 IF (IBUFF1 .GT. LBUFF1) THEN NBUFF1=NBUFF1+1 READ (JBUFF1, REC=NBUFF1)R1CC IBUFF1=IBUFF1-LBUFF1 END IF C RVALUE=R1CC(IBUFF1) R(NP,N)=RVALUE R(N,NP)=RVALUE 2 CONTINUE C 1 CONTINUE C I=0 C DO 3 N=1,NRA C DO 4 NP=1,NRA I=I+1 R1(I)=R(NP,N) 4 CONTINUE C 3 CONTINUE C C IF THE END OF THE BUFFER HAS BEEN REACHED AND C THE LAST LOOP ON L NOT EXECUTED READ THE NEXT C BLOCK OF INTEGRALS IN READINESS C IF (IBUFF1+1 .GT. LBUFF1 .AND. L .LT. LASTL) THEN NBUFF1=NBUFF1+1 READ (JBUFF1, REC=NBUFF1)R1CC END IF C NRECR1=NBUFF1 C RETURN END ********************************************************************** SUBROUTINE LOCCCD(L1,L1P,N1,N1P,LAM,RM) C C LOCATE DIRECT CONTINUUM CONTINUUM RK INTEGRALS C SUBROUTINE RDRK(L,LP) MUST HAVE BEEN CALLED C PREVIOUSLY TO SET UP THE POINTER ARRAY ICTCCD C AND THE STORAGE ARRAY RKCCD C C IF L .NE. LP RM IS OF LENGTH NRANG2 * NRANG2 C IF L .EQ. LP RM CONTAINS LOWER TRIANGLE C THE INTEGRALS ARE STORED IN THE ORDER: C L1 .LE. L1P C IF L1 .EQ. L1P THEN N1 .GE. N1P C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG C DIMENSION RM(INDIM2) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 3001 FORMAT(' **** STOP IN LOCCCD ****'/ 1 ' NO POINTER FOUND IN ICTCCD FOR'/ 2 ' L1=',I4,'L1P=',I4,'LAM=',I4) MAXNC1=MAXNC(L1) MAXNCP=MAXNC(L1P) NI=N1-MAXNC1 NJ=N1P-MAXNCP C IF (L1 .LT. L1P) THEN ICCD=ICTCCD(L1,L1P,LAM) NN1P=MAXNHF(L1P)-MAXNCP NIJ=(NI-1)*NN1P + NJ ELSE IF (L1 .GT. L1P) THEN ICCD=ICTCCD(L1P,L1,LAM) NN1P=MAXNHF(L1) - MAXNC1 NIJ=(NJ-1)*NN1P + NI ELSE ICCD=ICTCCD(L1,L1P,LAM) IF (NI .GE. NJ) THEN NIJ=(NI*(NI-1))/2 + NJ ELSE NIJ=(NJ*(NJ-1))/2 + NI END IF END IF C IF (ICCD .EQ. -999) THEN if (iam.eq.0) WRITE (IWRITE,3001) L1,L1P,LAM END IF C NIJ=(NIJ-1)*ICCLNG + ICCD - ICCDIF C DO 1 I=1,ICCLNG NIJ=NIJ+1 RM(I)=RKCCD(NIJ) 1 CONTINUE C RETURN END ********************************************************************** SUBROUTINE LOCMBB(IBBPOL,L1,L1P,N1,N1P,LAM,RVAL) C C TO LOCATE THE BOUND BOUND MULTIPOLE INTEGRAL C C THE INTEGRALS ARE STORED WITH THE ORDER: C L1 .LE. L1P C IF L1 .EQ. L1P THEN N1 .GE. N1P C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM6=(IDIM4+1)/2) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG C DIMENSION IBBPOL(IDIM1,IDIM1,KDIM6) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 3001 FORMAT(' **** STOP IN LOCMBB ****'/ 1 ' NO POINTER FOUND IN IBBPOL FOR'/ 2 ' L1=',I4,'L1P=',I4,'LAM=',I4) C LM=(LAM+1)/2 MAXNC1=MAXNC(L1) MAXNCP=MAXNC(L1P) NI=N1-MAXNC1 NJ=N1P-MAXNCP C IF (L1 .LT. L1P) THEN IBB=IBBPOL(L1,L1P,LM) MAXN=MAXNHF(L1P)-MAXNCP NIJ=(NI-1)*MAXN + NJ ELSE IF (L1 .GT. L1P) THEN IBB=IBBPOL(L1P,L1,LM) MAXN=MAXNHF(L1)-MAXNC1 NIJ=(NJ-1)*MAXN + NI ELSE IBB=IBBPOL(L1,L1P,LM) IF (NI .GE. NJ) THEN NIJ=(NI*(NI-1))/2 + NJ ELSE NIJ=(NJ*(NJ-1))/2 + NI END IF END IF C IF (IBB .EQ. -999) THEN if (iam.eq.0) WRITE (IWRITE,3001) L1,L1P,LAM END IF C NIJ=NIJ+IBB-ICCDIF RVAL=RKCCD(NIJ) RETURN END C*********************************************************************** SUBROUTINE ORDER(EN,NCHAN,NORDER) IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C --- DEFINE NORDER, WHICH POINTS TO VALUES OF EN IN ASCENDING ORDER. C C --- INPUT PARAMETERS ... C (EN(M),M=1,NCHAN) = ENERGIES IN ARBITRARY ORDER; C --- OUTPUT PARAMETER ... C (NORDER(M),M=1,NCHAN) = THE POSITION OF THE M-TH BIGGEST ENERGY C IN THE EN ARRAY. C DIMENSION EN(NCHAN),NORDER(NCHAN) C NORDER(1)=1 IF(NCHAN.EQ.1) RETURN DO 9 M=2,NCHAN J=M X=EN(J)+1.D-7 J1=J-1 DO 7 I=1,J1 IF(J.LT.M) GO TO 6 I1=NORDER(I) IF(X.GT.EN(I1)) GO TO 7 6 J=J-1 NORDER(J+1)=NORDER(J) 7 CONTINUE 8 NORDER(J)=M 9 CONTINUE RETURN END ********************************************************************** SUBROUTINE RDHMAT(ISP,MORE,HMAT,irowsize,icolsize) C C READS THE HAMILTONIAN MATRIX ELEMENTS STORED BY STHMAT C AND ARRANGES THEM IN HMAT READY FOR THE C DIAGONALISATION ROUTINE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS) PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(INDIM3=INDIM2+IDIM4) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) C C CHARACTER*4 PARITY(0:1) C COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2), 1 ICONCH(ITDIM3),KCONCH(ITDIM3), 2 L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1), 3 NSCHAN(IPDIM1),NSPREC(IPDIM1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/SETL /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) COMMON/STG1 /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2) C c **** dario **** c PARAMETER (LWORK=2*IHDIM1*IHDIM1+6*IHDIM1+1) !DSYEVD c DIMENSION WORK(LWORK),B(IHDIM1,IHDIM1) !DSYEVD c EQUIVALENCE (WORK(1),A(1,1)) !DSYEVD c COMMON/INTHAM/HMAT(IHDIM4) c PARAMETER (LIWORK=10*IHDIM1) !LAPACK c COMMON /LAP1/A(IHDIM1,IHDIM1),XA(IHDIM1,IHDIM1),AUX(IHDIM1,6),DUM !LAPACK c DIMENSION ISUPP(2*IHDIM1),IWORK(LIWORK) !LAPACK include 'mpif.h' common/crashblock/istop common/parablock/iam,nproc common/procdat/ictxt,myrow,mycol dimension HMAT(irowsize,icolsize) c **** dario **** C DIMENSION BUF1D(INDIM3),H1MAT(IDIM3,IDIM3) C DIMENSION CF(ITDIM3,IDIM4) DIMENSION CF(ITDIM6,ITDIM6,IDIM4) DIMENSION ISTCHS(ITDIM6) C DATA PARITY/'EVEN','ODD'/ C 1000 FORMAT(/,28X,'SUBROUTINE RDHMAT'/28X,'-----------------'/) 1001 FORMAT(' ASYMPTOTIC COEFFS FOR LAMBDA = 1 TO LAMAX FOR ', 1 'TARGET SPIN (2S+1) =',I3/) 1002 FORMAT(' CHANNEL',I3,' (STATE',I3,', L=',I3, 1 '), WITH CHANNEL',I3,' (STATE',I3,', L=',I3,')') 1003 FORMAT (' SURFACE AMPLITUDE') 1004 FORMAT(' HAMILTONIAN MATRIX FOR TARGET SPIN (2S+1) =',I3/) 1005 FORMAT( 2X,8F14.7) 1006 FORMAT(16X,7F14.7) 1007 FORMAT(30X,6F14.7) 1008 FORMAT(44X,5F14.7) 1009 FORMAT(58X,4F14.7) 1010 FORMAT(72X,3F14.7) 1011 FORMAT(86X,2F14.7) 1012 FORMAT(100X,F14.7) 1013 FORMAT(' CONTRIBUTION FROM CHANNEL',I3,' (STATE',I3,', L=',I3, 1 '), WITH CHANNEL',I3,' (STATE',I3,', L=',I3,')') 1014 FORMAT (/' EIGENVALUE=',F12.7,' RYD - RELATIVE TO TARGET GROUND', 1 F12.7,' AU') 1015 FORMAT (' VECTOR') C 2001 FORMAT (5F14.7) C IF (NBUG4 .GT. 0 .OR. NBUG9 .GT. 0) THEN if (iam.eq.0)WRITE(IWRITE,1000) END IF IF (NBUG4 .GT. 0) THEN C C CONSTRUCT THE ARRAY ISTCHS(ICH), THE STATE NUMBER COUPLED TO C CHANNEL ICH C ICH=0 DO 10 IST=1,NAST DO 11 K=1,KCONAT(IST,ISP) ICH=ICH+1 ISTCHS(ICH)=IST 11 CONTINUE 10 CONTINUE END IF C IF (NBUG4 .GT. 1) THEN if (iam.eq.0)WRITE(IWRITE,1004) NSPVAL(ISP) END IF C C THE RECORD NUMBERS WERE CALCULATED IN STHMAT TO BE C IN THE CORRECT ORDER TO FILL THE HMAT ARRAY. C ONLY THE UPPER HALF OF THE SUBMATRICES ON THE C DIAGONAL OF THE HAMILTONIAN MUST BE STORED, C ALL OTHERS REMAIN RECTANGULAR. C NCH=NSCHAN(ISP) IF (NCH .EQ. 0 .AND. MORE .NE. 0) RETURN IREC=NSPREC(ISP) IREC2=NSPREC(ISP+1) IF (NCH .NE. 0) THEN READ (JBUFD, REC=IREC) BUF1D c call mpi_barrier(mpi_comm_world,ierr) c **** dario **** c call mpi_bcast(buf1d,indim3,mpi_real8,0,mpi_comm_world,ierr) c **** dario **** END IF NY=LENHAM(ISP) NY2= (NY*(NY+1))/2 C C WRITE ASYMPTOTIC FILE C if (iam.eq.0)WRITE(ITAPE3) LRGL,-NSPVAL(ISP),NPTY,NCH,NY,MORE if (iam.eq.0)WRITE(ITAPE3)(KCONAT(I,ISP),I=1,NAST) if (iam.eq.0)WRITE(ITAPE3)(L2PSPN(ICH,ISP),ICH=1,NCH) IF (NCH .EQ. 0) THEN if (iam.eq.0)WRITE(ITAPE3)ZERO if (iam.eq.0)WRITE(ITAPE3)ZERO if (iam.eq.0)WRITE(ITAPE3)ZERO RETURN END IF C C ZEROISE HAMILTONIAN MATRIX C c **** parallel **** c DO 1 I=1,NY2 c HMAT(I)=ZERO c 1 CONTINUE c **** parallel **** C C PREPARE TO FILL CF BUFFER WITH ASYMPTOTIC COEFFS C FROM THE START OF EACH HAMILTONIAN BLOCK. C THERE IS A COEFFICIENT FOR EACH VALUE OF LAMBDA FROM C ONE TO LAMAX C C IBUF=0 C C PREPARE TO FILL HMAT C SET INITIAL ROW LENGTH AND STARTING POSITION FOR C FIRST ROW OF SUB MATRICES C IROW = NY ISTART= 1 DO 2 ICH=1,NCH L=L2PSPN(ICH,ISP)+1 NRA1=NVARY(L) C C SET STARTING POSITION FOR FIRST ROW OF SUB MATRICES C JSTART=ISTART C DO 3 JCH=ICH,NCH LP=L2PSPN(JCH,ISP)+1 NRA2=NVARY(LP) DO 31 IB=1,LAMAX C CF(JCH,IB)=BUF1D(IB) CF(ICH,JCH,IB)=BUF1D(IB) CF(JCH,ICH,IB)=BUF1D(IB) 31 CONTINUE C C COPY THE SUB MATRIX FROM THE BUFFER TO THE C ARRAY H1MAT, DIMENSION NRA1 X NRA2 C IB=LAMAX DO 4 N=1,NRA1 C DO 5 NP=1,NRA2 IB=IB+1 H1MAT(NP,N)=BUF1D(IB) 5 CONTINUE C 4 CONTINUE C IREC = IREC+1 IF(IREC .LT. IREC2) THEN C C READ THE NEXT SUB MATRIX INTO THE BUFFER IN READINESS C READ (JBUFD, REC=IREC) BUF1D c call mpi_barrier(mpi_comm_world,ierr) c **** dario **** c call mpi_bcast(buf1d,indim3,mpi_real8,0,mpi_comm_world,ierr) c **** dario **** C END IF C IF(ICH .EQ. JCH) THEN C C ON THE DIAGONAL STORE UPPER HALF ONLY C NSTART = JSTART NROW = IROW C DO 61 N=1,NRA1 NPSTAR = NSTART C DO 71 NP=N,NRA2 c **** parallel **** c...........in serial stg3nx: c HMAT(NPSTAR) = H1MAT(NP,N) call locH(NY,NPSTAR,iprow,ipcol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HMAT(iprow,ipcol) = H1MAT(NP,N) c **** parallel **** NPSTAR = NPSTAR + 1 71 CONTINUE C NSTART = NSTART + NROW NROW = NROW - 1 61 CONTINUE JSTART=ISTART+NRA2 C ELSE C C OFF THE DIAGONAL STORE RECTANGULAR ARRAY C NROW = IROW - 1 NSTART = JSTART C DO 62 N=1,NRA1 NPSTAR = NSTART C DO 72 NP=1,NRA2 c **** parallel **** c...........in serial stg3nx: c HMAT(NPSTAR) = H1MAT(NP,N) call locH(NY,NPSTAR,iprow,ipcol,ipr,ipc) if (ipr.eq.myrow.and.ipc.eq.mycol) + HMAT(iprow,ipcol) = H1MAT(NP,N) c **** parallel **** NPSTAR = NPSTAR + 1 72 CONTINUE C NSTART = NSTART + NROW NROW = NROW - 1 62 CONTINUE JSTART=JSTART+NRA2 C END IF C IF (NBUG4 .GT. 1) THEN if (iam.eq.0)WRITE(IWRITE,1013) ICH,ISTCHS(ICH), 1 L2PSPN(ICH,ISP),JCH,ISTCHS(JCH),L2PSPN(JCH,ISP) IF (NBUG4 .EQ. 2) THEN NPR1=MIN(5,NRA1) NPR2=MIN(5,NRA2) ELSE NPR1=NRA1 NPR2=NRA2 END IF C IF (ICH .EQ. JCH) THEN C C PRINT TRIANGULAR MATRIX C NBLOCK=NPR1/8+1 NLO=1 NUP=MIN(8,NPR1) DO 12 NBLOC=1,NBLOCK IBLANK=0 DO 13 N=NLO,NUP NPLO=N NPUP=NUP DO 14 NBLOCP=NBLOC,NBLOCK IF (NBLOCP .EQ. NBLOC .AND. IBLANK .GT. 0) THEN IF (IBLANK .EQ. 1)THEN if (iam.eq.0)WRITE(IWRITE,1006)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE IF (IBLANK .EQ. 2) THEN if (iam.eq.0)WRITE(IWRITE,1007)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE IF (IBLANK .EQ. 3) THEN if (iam.eq.0)WRITE(IWRITE,1008)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE IF (IBLANK .EQ. 4) THEN if (iam.eq.0)WRITE(IWRITE,1009)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE IF (IBLANK .EQ. 5) THEN if (iam.eq.0)WRITE(IWRITE,1010)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE IF (IBLANK .EQ. 6) THEN if (iam.eq.0)WRITE(IWRITE,1011)(H1MAT(NP,N),NP=NPLO,NPUP) ELSE if (iam.eq.0)WRITE(IWRITE,1012)(H1MAT(NP,N),NP=NPLO,NPUP) END IF ELSE if (iam.eq.0)WRITE(IWRITE,1005)(H1MAT(NP,N),NP=NPLO,NPUP) END IF NPLO=NPUP+1 NPUP=MIN(NPUP+8,NPR2) 14 CONTINUE IBLANK=IBLANK+1 13 CONTINUE NLO=NUP+1 NUP=MIN(NUP+8,NPR1) 12 CONTINUE ELSE C C PRINT RECTANGULAR MATRIX C NBLOCK=NPR2/8+1 DO 15 N=1,NPR1 NPLO=1 NPUP=MIN(8,NPR2) DO 16 NBLOCP=1,NBLOCK if (iam.eq.0)WRITE(IWRITE,1005)(H1MAT(NP,N),NP=NPLO,NPUP) NPLO=NPLO+8 NPUP=MIN(NPUP+8,NPR2) 16 CONTINUE 15 CONTINUE END IF END IF C 3 CONTINUE C C WRITE CF MATRIX TO ASYMPTOTIC FILE C C WRITE (ITAPE3) ((CF(JCH,ILAM),ILAM=1,LAMAX),JCH=ICH,NCH) C ISTART = NPSTAR IROW = IROW - NRA1 C 2 CONTINUE C C WRITE CF MATRIX TO ASYMPTOTIC FILE C C if (iam.eq.0) + WRITE (ITAPE3) (((CF(ICH,JCH,ILAM),ICH=1,NCH),JCH=1,NCH), 1 ILAM=1,LAMAX) IF (NBUG4 .GT. 0) THEN if (iam.eq.0)WRITE(IWRITE,1001) NSPVAL(ISP) DO 17 ICH=1,NCH DO 18 JCH=ICH,NCH if (iam.eq.0)WRITE(IWRITE,1002)ICH,ISTCHS(ICH),L2PSPN(ICH,ISP), 1 JCH,ISTCHS(JCH),L2PSPN(JCH,ISP) if (iam.eq.0)WRITE(IWRITE,1005)(CF(ICH,JCH,ILAM),ILAM=1,LAMAX) 18 CONTINUE 17 CONTINUE END IF RETURN C END ********************************************************************** SUBROUTINE RDRINT(ISTART,IFIN,NB,NRECR) C C A DOUBLE BUFFERING READ ROUTINE C READS RADIAL INTEGRALS FROM FILE JBUFF1 INTO ARRAY C RKCCD STARTING AT FILE POSITION ISTART TO FILE POSITION C IFIN. CALLED FROM CRLMAT TO READ BOUND BOUND MULTIPOLE C INTEGRALS AND FROM RDRKCC TO READ RK INTEGRALS FOR GIVEN L,LP C C NRECR IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE C NB IS THE BUFFER INDICATOR IN THE DOUBLE BUFFERING C NRECR AND NB ARE SET ZERO ON FIRST ENTRY C THEY ARE NOT SET THEREAFTER BUT CONTAIN THE VALUES C LEFT BY THE PREVIOUS EXECUTION C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/RADBUF/BUFF1(LBUFF1),BUFF2(LBUFF1) COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 3000 FORMAT(' **** STOP IN RDRINT ****'/ 1 ' DIMENSION OF RKCCD IS NOT LARGE ENOUGH'/ 2 ' L=',I3,'LP=',I3,'IRDIM1=',I7,'ISTART=',I7,'IFIN=',I7/ 3 'TRY INCREASING MZMEG TO:',I3) C C C CHECK WHETHER DIMENSION OF RKCCD IS LARGE ENOUGH C IF ( IFIN-ISTART+1 .GT. IRDIM2 ) THEN MXMEM=IRDIM2/1000000+1 if (iam.eq.0)WRITE(IWRITE,3000) IRDIM2,ISTART,IFIN,MXMEM STOP END IF C C STORE ICCDIF IN COMMON READY TO BE SUBTRACTED FROM C STARTING VALUES IN POINTER ARRAY C ICCDIF = ISTART C NBLK1 = ISTART/LBUFF1 NBLK2 = IFIN/LBUFF1 IB1 = MOD(ISTART,LBUFF1) IB2 = MOD(IFIN,LBUFF1) IF (IB1 .EQ. 0) THEN NBLK1=NBLK1-1 IB1=LBUFF1 END IF IF (IB2 .EQ. 0) THEN ! CHECK IB2 AS WELL - NRB NBLK2=NBLK2-1 IB2=LBUFF1 END IF C C READ IN THE BLOCKS OF INTEGRALS USING TWO BUFFERS C BUT CHECK WHETHER FIRST BLOCK IS ALREADY IN THE BUFFER C IF (NRECR .NE. NBLK1) THEN READ (JBUFF1, REC=NBLK1) BUFF1 NB=0 END IF I=0 IBLO=IB1 C DO 1 NBLK=NBLK1+1, NBLK2 IF(NB .EQ. 0) THEN READ( JBUFF1, REC=NBLK) BUFF2 DO 21 IB=IBLO,LBUFF1 I=I+1 RKCCD(I)=BUFF1(IB) 21 CONTINUE NB=1 ELSE READ( JBUFF1, REC=NBLK) BUFF1 DO 22 IB=IBLO,LBUFF1 I=I+1 RKCCD(I)=BUFF2(IB) 22 CONTINUE NB=0 END IF IBLO=1 1 CONTINUE C C COPY THE LAST BUFFER INTO RKCCD C IF (NB .EQ. 0) THEN DO 31 IB=IBLO,IB2 I=I+1 RKCCD(I)=BUFF1(IB) 31 CONTINUE ELSE DO 32 IB=IBLO,IB2 I=I+1 RKCCD(I)=BUFF2(IB) 32 CONTINUE END IF NRECR=NBLK2 C RETURN END ********************************************************************** SUBROUTINE RDRKCC(L,LP,NRECIC,NB,NRECRK) C C FOR THIS L,LP COMBINATION STORE ALL CONTINUUM/CONTINUUM C RK INTEGRALS IN RKCCD C C NRECIC IS THE RECORD NUMBER OF THE BLOCK CURRENTLY IN STORE C NRECRK IS THE RECORD NUMBER OF THE RK BLOCK CURRENTLY IN STORE C NB IS THE BUFFER INDICATOR IN THE DOUBLE BUFFERING OF THE RK C NRECIC, NRECRK AND NB ARE SET ZERO ON FIRST ENTRY C THEY ARE NOT SET THEREAFTER BUT CONTAIN THE VALUES C LEFT BY THE PREVIOUS EXECUTION C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) C PARAMETER(LBUFF1=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 1001 FORMAT(' **** ERROR IN RDRKCC FOR L,LP=',2I4/ 1 ' CANNOT FIND STARTING RECORD FOR L1,L1P,LAMMIN=',3I4/ 2 ' TOTAL NUMBER OF RK INTEGRALS STORED FOR L,LP=',I4) C C IF TRIANGULAR RULES FOR LAMMIN,L,LP AND LAMMIN,L1,L1P ARE C NOT SATISFIED RETURN C STORE POINTER ARRAY C CALCULATE START AND END POSITIONS IN THE STORAGE OF THE C RADIAL INTEGRALS FOR THIS L, LP C IF (L .EQ. LP) THEN NRA=NVARY(L+1) ICCLNG = (NRA*(NRA+1))/2 ELSE NRA1=NVARY(L+1) NRA2=NVARY(LP+1) ICCLNG = NRA1 * NRA2 END IF C LDIF=LP-L+1 LAMMIN=LDIF IF (LAMMIN .GT. 2*LRANG1-1) RETURN C NRECIC=L*LLPDIM + LDIF READ (JBUFF2, REC=NRECIC)(((ICTCCD(I,J,K),I=1,LRANG1), 1 J=1,LRANG1),K=1,LRANG1+LRANG1-1) C C FIND START AND END POSITIONS OF RK INTEGRALS C L1=1 L1P=LAMMIN IF (L1P .GT. LRANG1) THEN L1=LAMMIN-LRANG1+1 L1P=LRANG1 END IF ICCD1=ICTCCD(L1,L1P,LAMMIN) NRKVAL=NCCD(LDIF,L+1) ICCD2=ICCD1+NRKVAL-1 C IF (ICCD1 .LT. LBUFF1) THEN if (iam.eq.0)WRITE(IWRITE,1001) L,LP,L1,L1P,LAMMIN,NRKVAL CALL EXIT (0) END IF C CALL RDRINT(ICCD1,ICCD2,NB,NRECRK) C RETURN END ********************************************************************** SUBROUTINE RDVSH(KIS,KJS,NRECV1,NRECV2,NVSHEL) C C TO READ IN VSH FILE FOR SETS KIS,KJS C NRECV1 IS THE RECORD IN BUFV1 C NRECV2 IS THE RECORD IN BUFV2 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM3=MZNR2) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF) C PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGBUF/BUFV1(LBUFFV),BUFV2(LBUFFV), 1 IBUFV1(LBUFFV),IBUFV2(LBUFFV) COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 1001 FORMAT(2I5,5X,I8,5X,I8,5X,E13.5) C 3000 FORMAT (' **** STOP IN RDVSH **** '/ 1 ' DIMENSION OF VSH NOT LARGE ENOUGH'/ 2 ' KIS=',I4,' KJS=',I4,' IVDIM1=',I7,' ICONT1=',I8, 3 ' ICONT2=',I8) C KISJS=NSET*(KIS-1) - (KIS*(KIS-1))/2 +KJS ICONT1=IVCONT(KISJS) ICONT2=IVCONT(KISJS+1) -1 NVSHEL=ICONT2-ICONT1+1 IF (ICONT2-ICONT1 .GE. IVDIM1) THEN if (iam.eq.0)WRITE(IWRITE, 3000) KIS,KJS,IVDIM1,ICONT1,ICONT2 STOP END IF C NBLK1=ICONT1/LBUFFV NBLK2=ICONT2/LBUFFV IB1=MOD(ICONT1,LBUFFV) IB2=MOD(ICONT2,LBUFFV) IF (IB1 .EQ. 0) THEN NBLK1=NBLK1-1 IB1=LBUFFV END IF C IF (NRECV1 .NE. NBLK1) THEN IF (NRECV2 .NE. NBLK1) THEN READ (JBUFFV, REC=NBLK1) BUFV1 READ (JBUFIR, REC=NBLK1) IBUFV1 NB=0 NRECV1=NBLK1 ELSE NB=1 END IF ELSE NB=0 END IF C I=0 IBLO=IB1 DO 1 NBLK=NBLK1+1,NBLK2 IF (NB .EQ. 0) THEN IF (NRECV2 .NE. NBLK) THEN READ (JBUFFV, REC=NBLK) BUFV2 READ (JBUFIR, REC=NBLK) IBUFV2 NRECV2=NBLK END IF NB=1 DO 21 IB=IBLO,LBUFFV I=I+1 VSH(I)=BUFV1(IB) IRHSGL(I)=IBUFV1(IB) 21 CONTINUE ELSE IF (NRECV1 .NE. NBLK) THEN READ (JBUFFV, REC=NBLK) BUFV1 READ (JBUFIR, REC=NBLK) IBUFV1 NRECV1=NBLK END IF NB=0 DO 22 IB=IBLO,LBUFFV I=I+1 VSH(I)=BUFV2(IB) IRHSGL(I)=IBUFV2(IB) 22 CONTINUE END IF IBLO=1 1 CONTINUE C C COPY THE LAST BUFFER C IF (NB .EQ. 0) THEN DO 31 IB=IBLO,IB2 I=I+1 VSH(I)=BUFV1(IB) IRHSGL(I)=IBUFV1(IB) 31 CONTINUE NB=1 ELSE DO 32 IB=IBLO,IB2 I=I+1 VSH(I)=BUFV2(IB) IRHSGL(I)=IBUFV2(IB) 32 CONTINUE NB=0 END IF IF (NBUG5 .EQ. 1) THEN IR0=0 !CAN BE RESET TO MATCH PARTICULAR STG1NX if (iam.eq.0)WRITE(IWRITE,1001) X (KIS,KJS,MOD(IR0+IR-1,LBUFFV)+1,IRHSGL(IR),VSH(IR),IR=1,I) ENDIF C RETURN END ********************************************************************** SUBROUTINE RECOV1(I,IDAT,ICURR) * * THIS ROUTINE IS CALLED ONLY IN THE CASE OF ARRAY OVERFLOW. * IF IPLACE=0 THE PROGRAM STOPS HERE, OTHERWISE THE PROGRAM RETURNS * TO THE CALLING ROUTINE AFTER SETTING IPLACE=I. * * I IS THE POSITION IN THE IDRAY AND IPRAY LISTS HOLDING * THE NAME OF THE RELEVANT DIMENSION PARAMETER AND * PREPROCESSOR PARAMETER. * IDAT IS THE ARRAY SIZE REQUIRED BY THE DATA * ICURR IS THE CURRENT ARRAY SIZE * IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' CHARACTER*6 IDRAY CHARACTER*3 IPRAY C C PARAMETER(IDIM1=MZLR1,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 3 LBUFF1=MZBUF) C PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) C COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint 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'/ c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** 1000 FORMAT(/' * ARRAY OVERFLOW *'/ 1 /1X,A6,' (MZ',A3,') SHOULD BE INCREASED FROM', 2 I7,' TO AT LEAST ',I7) * 1001 FORMAT(/' PROGRAM TERMINATES IN RECOV1'/) if (iam.eq.0)WRITE(IWRITE,1000)IDRAY(I),IPRAY(I),ICURR,IDAT 1 if (iam.eq.0)WRITE(IWRITE,1001) CALL EXIT (0) C RETURN END C*********************************************************************** SUBROUTINE RECSET(LLLO,LLUP,ICONTN) C C READS STGANG DATA AND RECONSTRUCTS COMMON/SETS/ 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 IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C CHARACTER*8 FILEA,FILEB,FILEC,FILED CHARACTER*4 TITLE,CONT CHARACTER*3 RELOP LOGICAL exs C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1) PARAMETER(ICDIM1=MZNCF,ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(ICDIM3=ICDIM2+ICDIM2-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) C PARAMETER(LBUFFV=MZBUF,LBUFFZ=MZBUF) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) COMMON/SHELLS/NJCOMP(ICDIM2),LJCOMP(ICDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) C COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN COMMON /NRB/MAXC COMMON /CPBHEADER/LRANG2ALL,MINIMALL,LR2 DATA IFLAGR/0/ C NAMELIST/STGNX/MINLT,MAXLT,LRGLLO,LRGLUP,MAXC,LFIXN,MJS,NPTYIN X ,RELOP,IWAVE X ,NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 C DIMENSION NOCCSH(ICDIM1),NOCORB(ICDIM2,ICDIM1), 1 NELCSH(ICDIM2,ICDIM1),J1QNRD(ICDIM3,3,ICDIM1) DIMENSION TITLE(18) DIMENSION LVAL(ILDIM2),ILVAL(ILDIM2),NSCOL(ILDIM2) DIMENSION IWAVE(1128) SAVE IFLAGR c **** parallel **** include 'mpif.h' common/parablock/iam,nproc COMMON /buff/status(MPI_STATUS_SIZE) c **** parallel **** 1000 FORMAT(' ',10X,'STGANG DATA'//) 1001 FORMAT(' LSET ISPSET IPISET NCSET NCFGSE') 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'/) 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 DATA CONT/'CONT'/ C C NRB NPTYMN=0 NPTYMX=1 NPTYIN=-1 C C IREAD=7 IWRITE=6 ISIZENX=8 inquire(file='RAD1.DAT',exist=exs) if(IFLAGR.eq.0)then OPEN(UNIT=IREAD,FILE='dstgnx',STATUS='UNKNOWN') OPEN(UNIT=IWRITE,FILE='rout3nx',STATUS='UNKNOWN') if(.not.exs) X OPEN(UNIT=ISIZENX,FILE='sizeNX.dat',STATUS='UNKNOWN') endif JBUFIR=20 JBUFFV=21 JBUFFZ=22 FILEA='ANG1.DAT' FILEB='ANG2.DAT' FILEC='ANG3.DAT' REWIND (IREAD) READ(IREAD,2001) (TITLE(I),I=1,18) IF (TITLE(1) .EQ. CONT) THEN INX2=36 OPEN (INX2,FILE='NX2.DAT',ACCESS='SEQUENTIAL', 1 STATUS='OLD',FORM='UNFORMATTED') C OPEN (INX2,FILE='tapenx2',ACCESS='SEQUENTIAL', C 1 STATUS='OLD',FORM='UNFORMATTED') ICONTN=1 READ (INX2)MAXN,MAXORB,LRANG3,NSETS1,MAXNCF,NAST,MAXNST, 1 NCHSUM,MAXNCH,ISRAN3,NCFG IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1) ELSE ICONTN=0 READ(IREAD, *) IBUG1,IBUG2,IBUG3,IBUG4,IBUG5,IBUG6,IBUG7, 1 IBUG8,IBUG9 END IF C C OPEN ANGULAR FILES C OPEN(UNIT=JBUFIR,FILE=FILEA,ACCESS='DIRECT',STATUS='OLD', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFFV) OPEN(UNIT=JBUFFV,FILE=FILEB,ACCESS='DIRECT',STATUS='OLD', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFFV) OPEN(UNIT=JBUFFZ,FILE=FILEC,ACCESS='DIRECT',STATUS='OLD', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFFZ) C CW if (iam.eq.0)WRITE(IWRITE,1000) C MAXC=1 IF(ICONTN.EQ.1)THEN MAXLT=0 MINLT=0 LRGLLO=0 LRGLUP=0 NBUG1=0 NBUG2=0 NBUG3=0 NBUG4=0 NBUG5=0 NBUG6=0 NBUG7=0 NBUG8=0 NBUG9=0 C C READ(IREAD,STGNX) C C ***** PARALLEL CPB ******* C MINLTALL=MINLT MAXLTALL=MAXLT C C ***** PARALLEL CPB ******* C C if(.not.exs)READ(ISIZENX,*)MINLT,MAXLT C c write(0,*)iam,minlt,maxlt C IF(NPTYIN.EQ.0)NPTYMX=0 IF(NPTYIN.EQ.1)NPTYMN=1 IF(MAXLT.GT.0)LRGLUP=MAXLT IF(MINLT.GT.0)LRGLLO=MINLT ELSE READ (IREAD,*) LRGLLO,LRGLUP ENDIF LLLO=LRGLLO LLUP=LRGLUP CW WRITE (IWRITE,1005) LRGLLO,LRGLUP IF (ILDIM3 .LE. LRGLUP) CALL RECOV1(16,LRGLUP+1,ILDIM3) CW if (iam.eq.0)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 if (iam.eq.0)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) C ELSE C 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 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,*) 16 CONTINUE IREAD1=36 FILED='CONF.DAT' OPEN(UNIT=36,FILE=FILED,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 if (iam.eq.0)WRITE(IWRITE,1008)NCFG IF (ICDIM1 .LT. NCFG) CALL RECOV1(24,NCFG,ICDIM1) READ(IREAD1,*) (NOCCSH(I),I=1,NCFG) CW if (iam.eq.0)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 if (iam.eq.0)WRITE(IWRITE,1010) I,(NOCORB(J,I),J=1,N) READ(IREAD1,*) (NELCSH(J,I),J=1,N) CW if (iam.eq.0)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 if (iam.eq.0)WRITE(IWRITE,1012) CW if (iam.eq.0)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 if (iam.eq.0)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 NCFGSE(IS,ICS)=IC 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 NCFGSE(IS,ICS)=IC C 2 CONTINUE C NSET=IS C DO 4 I=1,NSET NC=NCSET(I) CW if (iam.eq.0)WRITE(IWRITE,1002)LSET(I),ISPSET(I),IPISET(I),NCSET(I), CW 1 (NCFGSE(I,J),J=1,NC) 4 CONTINUE LRANG3=LRANG3+1 IF (ILDIM1 .LT. LRANG3) CALL RECOV1(15,LRANG3,ILDIM1) LLPDIM=LRANG3+LRANG3-1 C C FIND LRANG2 NEEDED FOR HIGHEST MAXLT ***CPB PARALLEL*** C IF(IFLAGR.eq.0)then MINIMALL=999 LRANG2ALL=1 DO 775 LRGL=MINTLTALL,MAXLTALL DO 776 NPTY0=NPTYMN,NPTYMX IF(NPTYIN.LT.0)THEN NPTY=NPTY0 ELSE NPTY=MOD(LRGL,2)+NPTY0 NPTY=MOD(NPTY,2) ENDIF LMIN=999 LMAX=0 C DO 777 IS=1,NSET LCFG=LSET(IS) LMIN=MIN(LMIN,ABS(LCFG-LRGL)) LMAX=MAX(LMAX, LCFG) 777 CONTINUE LMAX=LMAX+LRGL ILTOT=LMAX-LMIN+1 C DO 778 KL=1,ILTOT ILVAL(KL)=LMIN+KL-1 778 CONTINUE C C LOOP OVER THIS RANGE OF CONTINUUM L C DO 779 KL=1,ILTOT L=ILVAL(KL) C DO 7710 IS=1,NSET LCFG=LSET(IS) IPI=IPISET(IS) C C PARITY TEST C IF(MOD(IPI+L+NPTY,2) .NE. 0) GO TO 7710 C C TRIANGLE RULE C IF(LRGL .GT. LCFG+L .OR. 1 LRGL .LT. ABS(LCFG-L)) GO TO 7710 NSCOL(KL)=1 GO TO 779 7710 CONTINUE C 779 CONTINUE KL=0 DO 7711 IKL=1,ILTOT IF (NSCOL(IKL) .NE. 0) THEN KL=KL+1 LVAL(KL)=ILVAL(IKL) END IF 7711 CONTINUE LTOT=KL IF (LTOT .EQ. 0) GO TO 776 LMIN=LVAL(1) LMAX=LVAL(LTOT) LRANG2ALL=MAX(LMAX+1, LRANG2ALL) MINIMALL=MIN(LMIN+1, MINIMALL) 776 CONTINUE 775 CONTINUE LR2=LRANG2ALL-MAXLTALL endif C C write(0,*)' LR2= ',LR2 C write(0,*)' LRANG2ALL =',LRANG2ALL C write(0,*)' MINIMALL =',MINIMALL C C **** END CPB PARALLEL ******** C C Now work out LRANG2 for the range needed C C FIND LRANG2 NEEDED FOR REQUIRED RANGE OF LRGL 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 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) 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 NSCOL(KL)=1 GO TO 9 10 CONTINUE C 9 CONTINUE 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 IF (LTOT .EQ. 0) GO TO 6 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 IF (IDIM2 .LT. LRANG2) CALL RECOV1(2,LRANG2,IDIM2) C IFLAGR=IFLAGR+1 RETURN END ********************************************************************** SUBROUTINE RETAP1 C C READS BASIC INFORMATION FROM THE FILE RAD3 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/STG1 /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2) C 1000 FORMAT(' ',16I4) 1001 FORMAT(' ',6E14.7) REWIND ITAPE1 ITAPE=ITAPE1 READ(ITAPE) RA,BSTO READ(ITAPE) (NVARY(L),L=1,LRANG2) READ(ITAPE) ((ENDS(N,L),N=1,NVARY(L)),L=1,LRANG2) READ(ITAPE) ((COEFF(I,L),I=1,3),L=1,LRANG2) C RETURN END ********************************************************************** SUBROUTINE SETCHN(LRGLLO) C C FOR GIVEN LRGL AND PARITY TO CALCULATE C TOTAL NUMBER OF CHANNELS, ORDER OF CHANNELS FOR L AND STATE C AND THE FIRST CHANNEL NUMBER FOR L AND SET C SET UP THE ARRAY ICONCH TO CONVERT CHANNEL NUMBERS C IN THE NEW NO EXCHANGE CODE TO THOSE IN THE RMATRIX C CODE C SET UP NCONAT NO. OF CHANNELS COUPLED TO EACH STATE C AND L2P ARRAY LVALUE OF EACH CHANNEL ORDERED AS IN THE C RMATRIX CODE C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(ISDIM1=MZSET) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) C COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2), 1 ICONCH(ITDIM3),KCONCH(ITDIM3), 2 L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1), 3 NSCHAN(IPDIM1),NSPREC(IPDIM1) COMMON/HPOINT/KRECH(0:ILDIM4) COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/SETL /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) C COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN C DIMENSION ICHAN(ILDIM2,ITDIM1) c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** C 1000 FORMAT(/,28X,'SUBROUTINE SETCHN'/28X,'-----------------'/) 1001 FORMAT(' IS KL NCHSET IST ICHAN ISTCH') 1002 FORMAT(' ',6I6) 1003 FORMAT(' ICONCH(NCHAN)') 1004 FORMAT(/' TARGET SPIN (2S+1) =',I3/' -----------------------') 1005 FORMAT(' NO. OF CHANNELS=',I4) 1006 FORMAT(' L2P COUPLED TO EACH CHANNEL') 1007 FORMAT(' NO. OF CHANNELS COUPLED TO EACH STATE') 1008 FORMAT(' NSPREC(NSP)') C IF (NBUG8 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1000) if (iam.eq.0)WRITE(IWRITE,1001) END IF C C ZEROISE ICHAN ARRAY C DO 6 I=1,NAST DO 7 KL=1,LTOT ICHAN(KL,I)=0 7 CONTINUE 6 CONTINUE C C C COUNT THE CHANNELS IN NCHAN C NCHAN=0 C C LOOP OVER CONTINUUM ANGULAR MOMENTUM C DO 1 KL=1,LTOT C C LOOP OVER SETS COUPLED TO THIS L=LVAL(KL) C DO 2 I=1,NSCOL(KL) IS=NSETL(I,KL) C C IF THERE ARE NO STATES COMPOSED FROM THIS SET GO TO NEXT SET C IF(NSTAT(IS) .EQ. 0) GO TO 2 C C STORE FIRST CHANNEL NUMBER INVOLVING THIS SET AND L C NCHSET(IS,KL)=NCHAN+1 C C LOOP OVER STATES COMPOSED FROM THIS SET C DO 3 J=1,NSTAT(IS) IST=NSETST(IS,J) NCHAN=NCHAN+1 IF (ITDIM3 .LT. NCHAN) CALL RECOV1(21,NCHAN,ITDIM3) C C STORE THE CHANNEL NUMBER FOR EACH L AND STATE TEMPORARILY C ICHAN(KL,IST)=NCHAN C C STORE THE STATE NUMBER COUPLED TO EACH CHANNEL C ISTCH(NCHAN)=IST IF (NBUG8 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1002)IS,KL,NCHSET(IS,KL), 1 IST,ICHAN(KL,IST),ISTCH(NCHAN) END IF 3 CONTINUE C 2 CONTINUE C 1 CONTINUE C C CONSTRUCT THE CHANNEL CONVERSION ARRAY C FOR ORDERING THE HAMILTONIAN ON DISC C STORE L2P AND NCONAT ARRAYS FOR ASYMPTOTIC PROGRAM C C KRECH(2*LRGL+NPTY) WILL HOLD THE STARTING RECORD NO. C FOR EACH LRGL, NPTY COMBINATION C IN THIS VERSION THE STARTING RECORD IS 3 FOR EVERY C LRGL,NPTY COMBINATION TO SAVE DISC SPACE. C C LENHAM(ISP) WILL HOLD THE LENGTH OF THE HAMILTONIAN FOR C SPIN ISP C C NSPREC(ISP) WILL HOLD THE STARTING RECORD NO. FOR C THE HAMILTONIAN FOR SPIN ISP C IF(NPTYIN.LT.0)THEN IRECH=LRGL+LRGL+NPTY C IF (IRECH .EQ. LRGLLO+LRGLLO) KRECH(IRECH)=3 ELSE IRECH=LRGL C IF (IRECH .EQ. LRGLLO) KRECH(IRECH)=3 ENDIF KRECH(IRECH)=3 NSPREC(1)=KRECH(IRECH)+1 C KCOUNT=0 C DO 8 ISP=1,NSP ISPIN=NSPVAL(ISP) LENHAM(ISP)=0 ICOUNT=0 C DO 4 IST=1,NAST C C IF THE STATE DOES NOT HAVE CURRENT SPIN ENTER ZERO IN KCONAT C LCOUNT=0 IF (LSPN(IST) .NE. ISPIN) GO TO 41 C DO 5 KL=1,LTOT ICH=ICHAN(KL,IST) IF(ICH .NE. 0) THEN KCOUNT=KCOUNT+1 ICOUNT=ICOUNT+1 LCOUNT=LCOUNT+1 ICONCH(ICH)=ICOUNT KCONCH(ICH)=KCOUNT L2PSPN(ICOUNT,ISP)=LVAL(KL) L=LVAL(KL)+1 LENHAM(ISP)=LENHAM(ISP)+NVARY(L) END IF 5 CONTINUE C 41 KCONAT(IST,ISP)=LCOUNT 4 CONTINUE C NSCHAN(ISP)=ICOUNT C IF (ITDIM6 .LT. ICOUNT) CALL RECOV1(22,ICOUNT,ITDIM6) C NSPREC(ISP+1)=(ICOUNT*(ICOUNT+1))/2 + NSPREC(ISP) 8 CONTINUE C KRECH(IRECH+1)=NSPREC(NSP+1) C IF (NBUG8 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1003) if (iam.eq.0)WRITE(IWRITE,1002)(ICONCH(I),I=1,NCHAN) END IF DO 9 ISP=1,NSP if (iam.eq.0)WRITE(IWRITE,1004)NSPVAL(ISP) if (iam.eq.0)WRITE(IWRITE,1005)NSCHAN(ISP) if (iam.eq.0)WRITE(IWRITE,1006) if (iam.eq.0)WRITE(IWRITE,1002)(L2PSPN(I,ISP),I=1,NSCHAN(ISP)) if (iam.eq.0)WRITE(IWRITE,1007) if (iam.eq.0)WRITE(IWRITE,1002)(KCONAT(I,ISP),I=1,NAST) 9 CONTINUE IF (NBUG8 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1008) if (iam.eq.0)WRITE(IWRITE,1002)(NSPREC(ISP),ISP=1,NSP+1) END IF C RETURN END ********************************************************************** SUBROUTINE STGHRD(ICONTN) C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C C READ STATE DATA AND FORM COMMON/STATE/ C CHARACTER*8 FILAS,FILEA,FILEB,FILEC,FILED C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX) PARAMETER(IPDIM1=MZSPN) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C PARAMETER(LBUFD=IDIM3*IDIM3+IDIM4) C COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/NBUG /NBUGI(9) COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) COMMON/STATE1/NAST0,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) C DIMENSION NCST(ITDIM1), NTCON(ITDIM1) DIMENSION EN(ITDIM1),NORDER(ITDIM1),BIJ(ITDIM1,ISDIM2), 1 LL1(ITDIM1),LSPN1(ITDIM1),LPTY1(ITDIM1) C DIMENSION TITLE(18) CHARACTER*4 TITLE CHARACTER*11 FILE1,FILE2,FILE3 !,FILEH CHARACTER*1 NUM(0:9) C C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DATA IFLAG/0/ C SAVE IFLAG C NRB NAMELIST/STG3A/NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8 X,NBUG9,IPOLPH,RAD,ISORT,iprint NAMELIST/STG3B/NBUT,NCONT,IDIAG,NAST,INAST,TOLER,NEST C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc namelist/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest c **** parallel **** 2000 FORMAT (18A4) C2001 FORMAT (5F14.7) C2002 FORMAT (9I5) C 1000 FORMAT (TR31,'***********'/TR31,'** **'/TR31,'** NXHAM **'/ 1 TR31,'** **'/TR31,'***********'//) 1002 FORMAT (////TR28,' SUBROUTINE STGHRD'/TR28,' -----------------') 1003 FORMAT (//' OUTPUT CHANNEL UNIT NUMBER=',I4/ 1 ' INPUT CHANNEL UNIT NUMBER=',I4/ 2 ' TARGET ANGULAR FILE--INTEGER UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 3 ' REAL UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 4 ' RACAH ANGULAR FILE UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 5 ' RADIAL INTEGRAL FILE UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 6 ' RADIAL POINTER FILE UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 7 ' SEQENTIAL FROM NXRAD UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 8 ' UNDIAG. HAMILTONIAN FILE UNIT NUMBER=',I4, 3 ' NAME= ',A8/ 9 ' OUTPUT FOR ASYMPTOTIC PROGRAM UNIT NUMBER=',I4, 3 ' NAME= ',A8) 1004 FORMAT (' DEBUG PARAMETERS'/12I5) 1005 FORMAT (' NUMBER OF ATOMIC STATES ',I4) 1006 FORMAT (' STATE',I3,' L=',I3,' SPIN=',I3,' PARITY=',I1) 1007 FORMAT (' NUMBER OF STATES WITH EACH SYMMETRY'/' ',9 I4) 1008 FORMAT (' STATES INDEX FOR EACH SYMMETRY') 1009 FORMAT (' ',9 I4) 1010 FORMAT (' AIJ=',8F14.7,4(1X/6X,8F14.7)) 1011 FORMAT (' WARNING - THE NORMALIZATION OF THE STATE IS ', 1 F14.7, /' THE STATE IS BEING RENORMALISED') 1012 FORMAT (' ENAT=',2F14.7,I8) 1013 FORMAT (' EXPERIMENTAL ENERGIES',4(1X,/6X,8F14.7)) C 3001 FORMAT (' **** ERROR IN DATA ****'/ 1 ' STATE ',I3,' SHOULD CONTAIN ',I3,' CONFIGURATIONS '/ 2 ' IT EXPECTS ',I3,' COEFFICIENTS') C IREAD=7 IWRITE=6 JBUFD=26 ITAPE3=10 C C FILEA='ANG1.DAT' FILEB='ANG2.DAT' FILEC='ANG3.DAT' C i1=iflag/100 i2=(iflag-100*(iflag/100))/10 i3=iflag-(100*(iflag/100))-i2*10 C FILE1='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE2='RAD2.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE3='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) c c i1=iam/100 c i2=(iam-100*(iam/100))/10 c i3=iam-(100*(iam/100))-i2*10 c FILEH='HAM1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) c FILAS='H.DAT' C C ZEROISE KSPIN ARRAY USED IN SORTING AND STORING C STATE SPINS C DO 1 I=1,IPDIM1 KSPIN(I)=0 1 CONTINUE C C ZEROISE NSTAT ARRAY USED IN COUNTING STATES BELONGING C TO EACH SET C DO 2 IS=1,NSET NSTAT(IS)=0 2 CONTINUE C IF (ICONTN .EQ. 1) THEN INX2=36 END IF C C OPEN HAMILTONIAN AND ASYMPTOTIC FILES C c OPEN(UNIT=JBUFD,FILE=FILEH,ACCESS='DIRECT',STATUS='UNKNOWN', OPEN(UNIT=JBUFD,ACCESS='DIRECT',STATUS='SCRATCH', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFD) if(iflag.eq.0)then OPEN(UNIT=ITAPE3,FILE=FILAS,ACCESS='SEQUENTIAL',STATUS='UNKNOWN', 1 FORM='UNFORMATTED') endif C if(iflag.eq.0)then if (iam.eq.0) WRITE (IWRITE,1000) if (iam.eq.0) WRITE (IWRITE,1002) if (iam.eq.0) WRITE (IWRITE,1003) x IWRITE,IREAD,JBUFIR,FILEA,JBUFFV,FILEB, 1 JBUFFZ,FILEC,JBUFF1,FILE1,JBUFF2,FILE2, 2 ITAPE1,FILE3,JBUFD,FILEH,ITAPE3,FILAS endif IF (ICONTN .NE. 1) THEN READ (IREAD,*) (NBUGI(I),I=1,9) END IF C if(iflag.eq.0)then if (iam.eq.0) WRITE (IWRITE,1004) (NBUGI(I),I=1,9) endif IF (ICONTN .EQ. 1) THEN READ (INX2) (LL1(I),LSPN1(I),LPTY1(I),I=1,NAST0) ELSE READ (IREAD,*) NAST0 END IF if(iflag.eq.0)then if (iam.eq.0) WRITE (IWRITE,1005) NAST0 endif IF (ITDIM1 .LT. NAST0) CALL RECOV1(19,NAST0,ITDIM1) C DO 3 IST=1,NAST0 IF(ICONTN .EQ. 0)READ (IREAD,*)LL1(IST),LSPN1(IST),LPTY1(IST) if ((iam.eq.0).and.(iflag.eq.0)) WRITE (IWRITE,1006) + IST,LL1(IST),LSPN1(IST),LPTY1(IST) LST=LL1(IST) LSPNST=LSPN1(IST) LPTYST=LPTY1(IST) IF (IPDIM1 .LT. LSPNST) CALL RECOV1(23,LSPNST,IPDIM1) C C RECORD THAT THIS SPIN VALUE OCCURS C KSPIN(LSPNST)=1 C C FIND TO WHICH SET THIS STATE BELONGS C DO 31 IS=1,NSET IF (LST .EQ. LSET(IS) .AND. 1 LSPNST .EQ. ISPSET(IS) .AND. 2 LPTYST .EQ. IPISET(IS)) THEN NSTAT(IS)=NSTAT(IS)+1 IF (ITDIM2 .LT. NSTAT(IS)) CALL RECOV1(20,NSTAT(IS),ITDIM2) NSETST(IS,NSTAT(IS))=IST NCST(IST)=NCSET(IS) END IF 31 CONTINUE C 3 CONTINUE C if ((iam.eq.0).and.(iflag.eq.0))then WRITE (IWRITE,1007)(NSTAT(IS), IS=1,NSET) endif if ((iam.eq.0).and.(iflag.eq.0)) WRITE (IWRITE,1008) C DO 4 IS=1,NSET if ((iam.eq.0).and.(iflag.eq.0))then WRITE(IWRITE,1009)(NSETST(IS,JST),JST=1,NSTAT(IS)) endif 4 CONTINUE C C RECORD SPIN VALUES OCCURRING C ICOUNT=0 C DO 5 ISP=1,IPDIM1 IF (KSPIN(ISP) .NE. 0) THEN ICOUNT=ICOUNT+1 NSPVAL(ICOUNT)=ISP KSPIN(ISP)=ICOUNT END IF 5 CONTINUE C NSP=ICOUNT C C READ NO. OF CONFIGURATIONS IN EACH STATE AND COMPARE C WITH THE NO. OF CONFIGURATIONS IN THE RELEVANT SET C BEFORE READING THE CONFIGURATION COEFFICIENTS AND C ENERGY FOR EACH STATE. TEST NORMALISATION C IF (ICONTN .EQ. 1) THEN READ (INX2)(NTCON(IST), IST=1,NAST0) ELSE READ (IREAD,*)(NTCON(IST), IST=1,NAST0) END IF PT01=TENTH*TENTH C DO 6 IST=1,NAST0 NCONF=NCST(IST) IF (NCONF .NE. NTCON(IST)) THEN if ((iam.eq.0).and.(iflag.eq.0)) WRITE (IWRITE,3001) IST,NCONF,NTCON(IST) STOP END IF IF (ICONTN .EQ. 1) THEN READ(INX2)(BIJ(IST,IC), IC=1,NCONF),EN(IST) ELSE READ (IREAD,*) (BIJ(IST,IC), IC=1,NCONF),EN(IST) END IF if (iam.eq.0.and.iflag.eq.0)then WRITE (IWRITE,1010) (BIJ(IST,IC), IC=1,NCONF) endif ANORM=ZERO DO 61 IC=1,NCONF ANORM=ANORM+BIJ(IST,IC)**2 61 CONTINUE C IF (ABS(ANORM-ONE) .GE. PT01) THEN if (iam.eq.0.and.iflag.eq.0) WRITE (IWRITE,1011) ANORM DO 62 IC=1,NCONF BIJ(IST,IC)=BIJ(IST,IC)/SQRT(ANORM) 62 CONTINUE if ((iam.eq.0).and.(iflag.eq.0))then WRITE (IWRITE,1010) (BIJ(IST,IC), IC=1,NCONF) endif END IF C T=TWO*(EN(IST)-EN(1)) if ((iam.eq.0).and.(iflag.eq.0)) WRITE (IWRITE,1012) EN(IST),T 6 CONTINUE C IF (ICONTN .EQ. 1) THEN C FILED='DATAR3' C C OUR 'dstg3' DATASET, USUALLY REDEFINED ON INPUT TO UNIT 5 IN STG3 C c **** parallel **** c FILED='dstg3' OPEN (37,FILE='dstg3',STATUS='OLD',ACCESS='SEQUENTIAL') rewind(37) c **** parallel **** c READ (37,2000) (TITLE(I),I=1,18) C READ (37,2002) (IDUMMY(I),I=1,9) C READ (37,2002) (IDUMMY(I),I=1,9) C READ (37,2002) (IDUMMY(I),I=1,3) C READ (37,2002) N1,N2,N3,N4,NEST,N5 NEST=0 NAST=0 ISORT=0 iprint=0 C READ(37,STG3A) READ(37,STG3B) c **** parallel **** READ(37,matrixdat) c **** parallel **** IF(NEST.EQ.0)NEST=NAST ELSE ISORT=0 READ(IREAD,*) NEST END IF NEST0=IABS(NEST) C IF (NEST .NE. 0) THEN IF (ICONTN .EQ. 1) THEN READ (37,*)(ENAT(I),I=1,NEST0) ELSE READ(IREAD,*)(ENAT(IST),IST=1,NEST0) END IF if (iam.eq.0.and.iflag.eq.0)then WRITE(IWRITE,1013)(ENAT(IST),IST=1,NEST0) endif C IF(ISORT.GT.0)CALL ORDER(EN,NAST0,NORDER) C DO I=1,NEST0 C RYDBERGS RELATIVE TO GROUND IF(NEST.GT.0)THEN ENAT(I)=HALF*ENAT(I)+EN(1) ELSE if (iam.eq.0) + WRITE(6,*)' NAST/NEST.LT.0 (MERGE) NOT IN NX CODE YET' STOP' NAST/NEST.LT.0 (MERGE) NOT IN NX CODE YET' ENDIF C C B.P. CODES NOW USE NAST/NEST.LT.0 TO MERGE C AND EST/EN(1).NE.0 FOR ABSOLUTE A.U. C AND EST/EN(1).GT.0 RELATIVE RYD C AND EST/EN(1).LT.0 RELATIVE CM-1 COLD IF(NEST.LT.0)ENAT(I)=HALF*ENAT(I)+EN(1) I0=I IF(ISORT.GT.0)I0=NORDER(I) EN(I0)=ENAT(I) ENDDO END IF C C ORDER THE STATES IN ORDER OF ASCENDING ENERGY C CALL ORDER(EN,NAST0,NORDER) C if (iam.eq.0.and.iflag.eq.0)then WRITE(IWRITE,'(/'' REORDERED BY ENERGY'',/)') endif C DO 8 IST=1,NAST0 JST=NORDER(IST) ENAT(IST)=EN(JST) LL(IST)=LL1(JST) LSPN(IST)=LSPN1(JST) LPTY(IST)=LPTY1(JST) NCONF=NCST(JST) C DO 9 IC=1,NCONF AIJ(IST,IC)=BIJ(JST,IC) 9 CONTINUE C if (iam.eq.0.and.iflag.eq.0)then WRITE(IWRITE,1010)(AIJ(IST,IC), IC=1,NCONF) endif 8 CONTINUE C DO 10 IST=1,NAST0 if (iam.eq.0.and.iflag.eq.0)then WRITE(IWRITE,1006)IST,LL(IST),LSPN(IST),LPTY(IST) endif 10 CONTINUE C DO 11 IST=1,NAST0 T=TWO*(ENAT(IST)-ENAT(1)) if (iam.eq.0.and.iflag.eq.0) WRITE (IWRITE,1012) ENAT(IST),T,IST 11 CONTINUE C DO 12 IS=1,NSET C DO 13 IST=1,NSTAT(IS) JST=NSETST(IS,IST) C DO 14 KST=1,NAST0 IF (NORDER(KST) .EQ. JST) THEN NSETST(IS,IST)=KST GO TO 13 END IF 14 CONTINUE C 13 CONTINUE if (iam.eq.0.and.iflag.eq.0)then WRITE(IWRITE,1009)(NSETST(IS,JST),JST=1,NSTAT(IS)) endif 12 CONTINUE C IFLAG=IFLAG+1 IF (ICONTN .EQ. 1) THEN CLOSE (INX2) CLOSE (37) goto 3456 END IF CLOSE (37) C 3456 continue RETURN END ********************************************************************** SUBROUTINE STGRRD * THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES * CERTAIN VARIABLES IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' CHARACTER*1 NUM(0:9) CHARACTER*11 FILE1,FILE2,FILE3 CHARACTER*72 LVALUE(4)*1 DATA LVALUE/'S','P','D','F'/ DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DATA IFLAG2/0/ C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO) PARAMETER(IDIM11=MZSLT) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(LBUFF2=IDIM1*IDIM1*KDIM9) PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM C DIMENSION MAXNLG(IDIM2),MAXPN(IDIM2) DIMENSION IRAD(IDIM11),ZE(IDIM11),C(IDIM11) C SAVE IFLAG2 c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** * FORMAT STATEMENTS 2000 FORMAT(A) 1002 FORMAT(//' ',10X,'STGRAD DATA') 1053 FORMAT(//12X,'THIS IS A NO-EXCHANGE CALCULATION'//) 1054 FORMAT(//12X,'THIS CALCULATION INCLUDES EXCHANGE'//) 1003 FORMAT(//' OUTPUT CHANNEL UNIT NUMBER =',I3/ * ' INPUT CHANNEL UNIT NUMBER =',I3/ * ' INPUT SEQUENTIAL UNIT NUMBER =',I3/ * ' OUTPUT SEQUENTIAL UNIT NUMBER =',I3/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR RADIAL INT *EGRALS)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR CONTINUUM- *CONTINUUM POINTER ARRAYS)'/) 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,' LAM =',I2,' IBC =',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) 3001 FORMAT(' **** NOT YET ALLOWED IN THIS VERSION ****') * SET THE CHANNEL UNIT NUMBERS IREAD=7 IWRITE=6 JBUFF1=23+(iam+1000) JBUFF2=24+(iam+2000) ITAPE1=25+(iam+3000) C i1=iflag2/100 i2=(iflag2-100*(iflag2/100))/10 i3=iflag2-(100*(iflag2/100))-i2*10 C C FILE1='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE2='RAD2.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE3='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) * READ IN THE BASIC DATA C C OPEN RADIAL FILES C OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT',STATUS='OLD', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFF1) OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT',STATUS='OLD', 1 FORM='UNFORMATTED',RECL=MACDIM*LBUFF2) OPEN(UNIT=ITAPE1,FILE=FILE3,ACCESS='SEQUENTIAL',STATUS='OLD', 1 FORM='UNFORMATTED') CW if (iam.eq.0)WRITE(IWRITE,1002) CW if (iam.eq.0) WRITE (IWRITE,1053) READ (IREAD, *) NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7, * NBUG8,NBUG9 READ (IREAD, *) IPSEUD,NUMERC IF (IPSEUD.EQ.0) THEN CW if (iam.eq.0)WRITE(IWRITE,1005) ELSE CW if (iam.eq.0)WRITE(IWRITE,1006) ENDIF IF (NUMERC.EQ.0) THEN CW if (iam.eq.0)WRITE(IWRITE,1007) ELSE CW if (iam.eq.0)WRITE(IWRITE,1008) ENDIF CW if (iam.eq.0)WRITE(IWRITE,1009) READ (IREAD, *) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC CW if (iam.eq.0)WRITE(IWRITE, 1010)NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM,IBC IF (IDIM1 .LT. LRANG1) CALL RECOV1(1,LRANG1,IDIM1) IF (IDIM3 .LT. NRANG2) CALL RECOV1(3,NRANG2,IDIM3) IF (IDIM4 .LT. LAMAX) CALL RECOV1(4,LAMAX,IDIM4) READ (IREAD, *) (MAXNHF(I),I=1,LRANG1) CW if (iam.eq.0)WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1) READ (IREAD, *) (MAXNLG(I),I=1,LRANG1) CW if (iam.eq.0)WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1) * READ IN DATA DEFINING THE STATIC POTENTIAL OF THE GROUND STATE CW if (iam.eq.0)WRITE(IWRITE,1013) DO 10,I=1,LRANG1 READ (IREAD, *) ISTAT 10 CONTINUE C IF (NUMERC.EQ.0) THEN * ANALYTIC ORBITALS ORBITALS ARE TO BE USED * READ IN COEFFICIENTS DEFINING THE ANALYTIC ORBITALS DO 20,L=1,LRANG1 DO 30,N=L,MAXNHF(L) READ (IREAD, *) M READ (IREAD, *) (IRAD(IM),IM=1,M) READ (IREAD, *) (ZE(IM),IM=1,M) READ (IREAD, *) (C(IM),IM=1,M) 30 CONTINUE 20 CONTINUE IF (IBC.NE.0) THEN READ (IREAD, *) RA,BSTO ENDIF CW if (iam.eq.0)WRITE(IWRITE,1018) RA,BSTO ENDIF IF (IPSEUD.EQ.0) THEN * THIS IS A NON-MODEL POTENTIAL CALCULATION. IF (IBC.GT.1) THEN * READ IN INTEGRATION MESH READ (IREAD, *) NIX READ (IREAD, *) HINT READ (IREAD, *) IHX READ (IREAD, *) IRX ENDIF * INITIALISE MAXNC IN NON-MODEL POTENTIAL * CALCULATION. DO 150,I=1,LRANG1 MAXNC(I)=I-1 150 CONTINUE C WHAT ABOUT MODEL? ENDIF CW if (iam.eq.0)WRITE(IWRITE,1023) (MAXNC(I),I=1,LRANG1) C C READ THE MAXIMUM ENERGY FOR THE BUTTLE CORRECTION C READ(IREAD, *) EK2MAX,MJS 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 C IF (LRANG2.GT.LRANG1) THEN DO 170,I=LRANG1+1,LRANG2 MAXNHF(I)=I-1 MAXNLG(I)=I-1 170 CONTINUE ENDIF IFLAG2=IFLAG2+1 RETURN END ********************************************************************** SUBROUTINE STGRRX * THIS ROUTINE READS IN THE INPUT DATA AND INITIALISES * CERTAIN VARIABLES USING RMATRX OUTPUT FILE IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' CHARACTER*1 NUM(0:9) CHARACTER*11 FILE1,FILE2,FILE3 CHARACTER*72 LVALUE(4)*1 LOGICAL exs DATA LVALUE/'S','P','D','F'/ DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DATA IFLAG3/0/ C PARAMETER(MACDIM=MZMAC) C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX,IDIM5=MZNPO) PARAMETER(IDIM6=MZNIX) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(LBUFF2=IDIM1*IDIM1*KDIM9) PARAMETER(LBUFF1=MZBUF,LBUFFV=MZBUF,LBUFFZ=MZBUF) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON /NRB/MAXC c **** parallel **** include 'mpif.h' common/parablock/iam,nproc SAVE IFLAG3 c **** parallel **** C DIMENSION MAXNLG(IDIM2),MAXPN(IDIM2),IHX(IDIM6),IRX(IDIM6) * FORMAT STATEMENTS 2000 FORMAT(A) 1002 FORMAT(//' ',10X,'STGRAD DATA') 1053 FORMAT(//12X,'THIS IS A NO-EXCHANGE CALCULATION'//) 1054 FORMAT(//12X,'THIS CALCULATION INCLUDES EXCHANGE'//) 1003 FORMAT(//' OUTPUT CHANNEL UNIT NUMBER =',I3/ * ' INPUT CHANNEL UNIT NUMBER =',I3/ * ' INPUT SEQUENTIAL UNIT NUMBER =',I3/ * ' OUTPUT SEQUENTIAL UNIT NUMBER =',I3/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR RADIAL INT *EGRALS)'/ * ' OUTPUT DIRECT ACCESS UNIT NUMBER =',I3,'(FOR CONTINUUM- *CONTINUUM POINTER ARRAYS)'/) 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,' LAM =',I2,' IBC =',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) 3001 FORMAT(' **** NOT YET ALLOWED IN THIS VERSION ****') * SET THE CHANNEL UNIT NUMBERS IREAD=7 IWRITE=6 JBUFF1=23+(iam+1000) JBUFF2=24+(iam+2000) ITAPE1=25+(iam+3000) INX1=35 C inquire(file='RAD1.DAT',exist=exs) if(.not.exs)then i1=iflag3/100 i2=(iflag3-100*(iflag3/100))/10 i3=iflag3-(100*(iflag3/100))-i2*10 C FILE1='RAD1.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE2='RAD2.DAT'//NUM(i1)//NUM(i2)//NUM(i3) FILE3='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) else FILE1='RAD1.DAT' FILE2='RAD2.DAT' FILE3='RAD3.DAT' endif C * OPEN THE DIRECT ACCESS FILES OPEN(UNIT=JBUFF1,FILE=FILE1,ACCESS='DIRECT', * STATUS='OLD',FORM='UNFORMATTED',RECL=MACDIM*LBUFF1) OPEN(UNIT=JBUFF2,FILE=FILE2,ACCESS='DIRECT', * STATUS='OLD',FORM='UNFORMATTED',RECL=MACDIM*LBUFF2) OPEN(ITAPE1,FILE=FILE3,ACCESS='SEQUENTIAL', * STATUS='OLD',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 1 FORM='UNFORMATTED') C CW if (iam.eq.0)WRITE(IWRITE,1002) C CW if (iam.eq.0)WRITE(IWRITE,1009) LAM=0 READ (INX1) LRANG1,NRANG2,LAMAX,NIX,NPTS,NRKPTS READ (INX1) NELC,NZ,LFIXN,IPSEUD IF(MAXC.GT.1)NRANG2=MAXC CW if (iam.eq.0)WRITE(IWRITE,1010) NELC,NZ,LRANG1,NRANG2,LFIXN,LAMAX,LAM IF (IDIM1 .LT. LRANG1) CALL RECOV1(1,LRANG1,IDIM1) IF (IDIM3 .LT. NRANG2) CALL RECOV1(3,NRANG2,IDIM3) IF (IDIM4 .LT. LAMAX) CALL RECOV1(4,LAMAX,IDIM4) C READ (INX1) (MAXNHF(I),I=1,LRANG1), (MAXNLG(I),I=1,LRANG1) CW if (iam.eq.0)WRITE(IWRITE,1011) (MAXNHF(I),I=1,LRANG1) CW if (iam.eq.0)WRITE(IWRITE,1012) (MAXNLG(I),I=1,LRANG1) C DO 150 I=1,LRANG1 MAXNC(I)=I-1 150 CONTINUE C SHOULD MAXNC BE REDEFINED FOR MODEL POTENTIAL? IF(IPSEUD.GT.0)THEN READ(INX1)HINT,(IHX(I),I=1,NIX),(IRX(I),I=1,NIX) C ANSWER APPEARS TO BE NO C READ(INX1)(MAXNC(I),I=1,LRANG1) ENDIF * DEFINE THE MAXPN ARRAY . NOT USED? DO 160,I=1,LRANG1 MAXPN(I)=MAXNHF(I)-MAXNC(I)+I-1 160 CONTINUE * 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 * INITIALISE THE VARIABLES LAMBC AND LAMCC (NOT USED HERE) IF (LAM.GE.2) THEN LAMBC=LAMAX ELSE LAMBC=0 ENDIF IF (LAM.GE.3) THEN LAMCC=LAMAX ELSE LAMCC=0 ENDIF C CLOSE (INX1) C IFLAG3=IFLAG3+1 RETURN END *======================================================================= SUBROUTINE STHMAT(LRGLLO) C C CALCULATES DIRECT TERMS IN THE HAMILTONIAN MATRIX C AND STORES THEM ON DISC IN RECORDS NRANG2 X NRANG2 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2,IDIM4=MZLMX) PARAMETER(ICDIM2=MZORB) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) PARAMETER(IRDIM1=MZMEG*1000000+MZKIL*1000+1) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS,ITDIM3=MZCHF,ITDIM6=MZCHS) C PARAMETER(LBUFF1=MZBUF,LBUFFZ=MZBUF) C PARAMETER(IDIM2=ILDIM3+ILDIM1-1) PARAMETER(KDIM9=IDIM1+IDIM1-1) PARAMETER(ILDIM2=ILDIM1+ILDIM1-1) PARAMETER(ILDIM4=ILDIM3+ILDIM3) PARAMETER(INDIM2=IDIM3*IDIM3) PARAMETER(ISDIM3=(ISDIM1*(ISDIM1+1))/2 +1) PARAMETER(ITDIM4=ITDIM2*ITDIM2) PARAMETER(IVDIM1=ILDIM1*(ISDIM2*(ISDIM2+1))/2 + 1 ILDIM1*ISDIM2*(ICDIM2-1)) PARAMETER(IHDIM1=ITDIM6*IDIM3) PARAMETER(IHDIM2=(IHDIM1*(IHDIM1+1))/2) PARAMETER(IHDIM3=INDIM2*ITDIM4) C PARAMETER(IHDIM4=MAX(IHDIM2,IRDIM1+IVDIM1+LBUFFZ+LBUFF1)) PARAMETER(IADIM=IHDIM2,IBDIM=IRDIM1+IVDIM1+LBUFFZ+LBUFF1) PARAMETER(ICDIM=(IADIM+IBDIM+1)/2) PARAMETER(IHDIM4=IADIM*(IADIM/ICDIM)+IBDIM*(IBDIM/ICDIM) 1 -(IADIM/IBDIM)*(IBDIM/IADIM)*IADIM) PARAMETER(IRDIM2=IHDIM4-IVDIM1-LBUFFZ-LBUFF1) C COMMON/ANGPNT/LAMIJ(ISDIM1,ISDIM1),IVCONT(ISDIM3),KRECZ(0:ILDIM4), 1 IRHSGL(IVDIM1) COMMON/ASYMPR/CRL(ITDIM1,ITDIM1,IDIM4) COMMON/CHANN /NCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2), 1 ICONCH(ITDIM3),KCONCH(ITDIM3), 2 L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1), 3 NSCHAN(IPDIM1),NSPREC(IPDIM1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/HPOINT/KRECH(0:ILDIM4) COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/INTHAM/BUFFZ(LBUFFZ),VSH(IVDIM1), 1 R1CC(LBUFF1),RKCCD(IRDIM2) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/RPOINT/IST3(IDIM2),ICTCCD(IDIM1,IDIM1,KDIM9),NCCD(ILDIM2, 1 IDIM2),ICCDIF,ICCLNG COMMON/SETL /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/SETS /NSET,LSET(ISDIM1),ISPSET(ISDIM1),IPISET(ISDIM1), 1 NCSET(ISDIM1),NCFGSE(ISDIM1,ISDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) C COMMON/NRBPTY/NPTYMN,NPTYMX,NPTYIN C PARAMETER(INDIM3=INDIM2+IDIM4) C DIMENSION CMAT(IDIM3,IDIM3),HSMAT(IHDIM3) DIMENSION BUF1D(INDIM3),BUF2D(INDIM3) DIMENSION R1(INDIM2) DIMENSION ZB(ILDIM1) DIMENSION CF(ITDIM4,IDIM4) C c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** C 1000 FORMAT(/,28X,'SUBROUTINE STHMAT'/28X,'-----------------'/) 1001 FORMAT(' L LP IS JS IC JC IST JST') 1002 FORMAT(4X,8I4) 1003 FORMAT(2X,9E13.6) 1004 FORMAT(' ONE-ELECTRON INTEGRAL') 1005 FORMAT(' Z-COEFFICIENTS') 1006 FORMAT(' CF MATRIX OF ASYMPTOTIC COEFFS FOR SETS ',I3,I3) C 3000 FORMAT(' **** STOP IN STHMAT ****'/ 1 ' ERROR IN USING VSH FILE FOR SETS ', 2I4) C IF (NBUG6 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1000) if (iam.eq.0)WRITE(IWRITE,1001) END IF C IF(NPTYIN.LT.0)THEN LTEST=LRGLLO+LRGLLO LRGLPI=LRGL+LRGL+NPTY ELSE LTEST=LRGLLO LRGLPI=LRGL ENDIF C C AT FIRST ENTRY READ POINTER ARRAYS FOR ANGULAR AND RK INTEGRALS C IF (LRGLPI .EQ. LTEST) THEN READ (JBUFFZ, REC=2) L,KRECZ READ (ITAPE1)((NCCD(I,J),I=1,LLPDIM),J=1,LRANG2) END IF C C WRITE POINTER BLOCK FOR THE HAMILTONIAN MATRICES FOR C THIS LRGL AND PARITY C IRECH=LRGLPI NRECH=KRECH(IRECH) WRITE(JBUFD, REC=NRECH) NSP,(NSPVAL(I), 1 NSCHAN(I),NSPREC(I),LENHAM(I),I=1,NSP) IF (LTOT .EQ. 0) RETURN C C C READ FIRST BLOCK OF Z COEFFICIENTS FOR THIS LRGL AND PARITY C IRECZ=LRGLPI NRECZ=KRECZ(IRECZ) READ (JBUFFZ, REC=NRECZ)BUFFZ IZCON=1 C C SET RECORD NUMBERS TO ZERO READY FOR INITIAL ENTRY C INTO ROUTINES TO READ RADIAL INTEGRALS FROM DISC C NRECR1=0 NRECIC=0 NRECRK=0 NB=0 NRECV1=0 !NRB NRECV2=0 !NRB C LASTL=LVAL(LTOT)+1 C C LOOP OVER CONTINUUM ANGULAR MOMENTUM C DO 2 KL=1,LTOT L=LVAL(KL) NRA1=NVARY(L+1) C C LOCATE ONE ELECTRON INTEGRAL C CALL LOC1CC(L+1,LASTL,NRECR1,R1) C IF (NBUG6 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1004) NRA1SQ=NRA1*NRA1 if (iam.eq.0)WRITE(IWRITE,1003)(R1(NNP),NNP=1,NRA1SQ) END IF C DO 3 KLP=KL,LTOT LP=LVAL(KLP) NRA2=NVARY(LP+1) NRASQR=NRA1*NRA2 C C READ FROM DISC ALL CONTINUUM-CONTINUUM RK INTEGRALS C FOR THIS L,LP COMBINATION C CALL RDRKCC(L,LP,NRECIC,NB,NRECRK) C C LOOP OVER THE SETS COUPLED TO L AND LP C DO 4 I=1,NSCOL(KL) IS=NSETL(I,KL) ICHAN1=NCHSET(IS,KL) LCFG=LSET(IS) C C FIND THE SPIN, THE NO. OF CHANNELS IN THE C HAMILTONIAN FOR THIS SPIN AND THE RECORD C IN THE FILE AT WHICH IT STARTS C ISPIN=ISPSET(IS) ISP=KSPIN(ISPIN) NCH=NSCHAN(ISP) NSREC=NSPREC(ISP)-1 C JLO=1 IF(L .EQ. LP) JLO=I C DO 5 J=JLO,NSCOL(KLP) JS=NSETL(J,KLP) JCHAN1=NCHSET(JS,KLP) LCFGP=LSET(JS) LAMCFG=LCFG+LCFGP JSPIN=ISPSET(JS) C C LAMSIN IS THE INDICATOR TO SAY SETS IS, JS CAN INTERACT C IF THEY CANNOT IT WILL BE SET TO -999 AND USED TO SKIP C UNNECESSARY LOOPS AND ENSURE ZERO SUB MATRICES ARE C STORED IN THE HAMILTONIAN ON DISC IN THE APPROPRIATE C CHANNEL POSITIONS C LAMSIN=0 C C SUBROUTINE DIRANG ASSUMES JS .GE. IS C TO BE CONSISTENT IF IS .GT. JS THE SETS MUST BE REVERSED C IF(IS .GT. JS) THEN KIS=JS KJS=IS ELSE KIS=IS KJS=JS END IF C LAMLU=LAMIJ(KIS,KJS) C C TEST WHETHER SPIN ALLOWS SETS TO COUPLE. IF NOT C JUMP TO END OF LOOP C IF (ISPIN .NE. JSPIN) GO TO 5 C C TEST WHETHER SETS CANNOT COUPLE FOR SOME OTHER REASON C IF THEY CAN LAMIJ CONTAINS THE LIMITS ON LAMBDA IMPOSED C BY SHELL COUPLING BETWEEN CONFIGURATIONS IN THE SETS C IF(LAMLU .EQ. -999) THEN LAMSIN=-999 GO TO 50 END IF C LAMLO=LAMLU/100 LAMUP=MOD(LAMLU,100) C C CALCULATE LAMST,LAMFIN THE LIMITS ON LAMBDA IN THE STORAGE C OF THE Z COEFFICIENTS C LAMST=MAX(ABS(L-LP),LAMLO) LAMFIN=MIN( L+LP ,LAMUP) IF(LAMST .GT. LAMFIN) THEN LAMSIN=-999 GO TO 50 END IF C C TRANSFER TARGET ANGULAR INTEGRALS FOR THIS SET C COMBINATION TO VSH C CALL RDVSH(KIS,KJS,NRECV1,NRECV2,NVSHEL) ICOUNT=1 C C TRANSFER Z COEFFICIENTS FOR THIS SET COMBINATION TO ZB C IZ=1 DO 51 LAM=LAMST,LAMFIN,2 ZB(IZ)=BUFFZ(IZCON) IZ=IZ+1 IZCON=IZCON+1 IF (IZCON .GT. LBUFFZ) THEN NRECZ=NRECZ+1 READ (JBUFFZ, REC=NRECZ) BUFFZ IZCON=1 END IF 51 CONTINUE C IF (NBUG6 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1005) if (iam.eq.0)WRITE(IWRITE,1003)(ZB(IIZ),IIZ=1,IZ) END IF C C C CALCULATE THE ASYMPTOTIC COEFFICIENTS C CF(IJCHAN,LAMBDA) FOR CHANNELS IN THE C SETS IS,JS. LAMBDA .GE. 1 C C ZEROISE CF ARRAY C 50 DO 52 ICFS=1,ITDIM4 DO 53 LAM=1,LAMAX CF(ICFS,LAM)=ZERO 53 CONTINUE 52 CONTINUE ICFS=0 IF (LAMLU .GT. 0 .AND. LAMSIN .EQ. 0) THEN C C CALCULATE LIMITS ON LAMBDA CONSISTENT WITH Z-COEFFS C IF (LAMST .EQ. 0) THEN LAMS=2 IZCONT=1 ELSE LAMS=LAMST IZCONT=0 END IF LAMF=MIN(LAMAX,LAMFIN) DO 501 IST=1,NSTAT(IS) ISTAT=NSETST(IS,IST) IF (ICHAN1 .EQ. JCHAN1) THEN JSTLO=IST ELSE JSTLO=1 END IF DO 502 JST=JSTLO,NSTAT(JS) JSTAT=NSETST(JS,JST) IF (JSTAT .LT. ISTAT) THEN C C INTERCHANGE STATE NUMBERS SINCE ONLY UPPER HALF C OF THE CRL MATRIX IS STORED C KJSTAT=ISTAT KISTAT=JSTAT ELSE KJSTAT=JSTAT KISTAT=ISTAT END IF IZ=IZCONT ICFS=ICFS+1 DO 503 LAM=LAMS,LAMF,2 IZ=IZ+1 CF(ICFS,LAM)=2*CRL(KISTAT,KJSTAT,LAM)*ZB(IZ) 503 CONTINUE 502 CONTINUE 501 CONTINUE IF (NBUG7 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1006) IS,JS if (iam.eq.0)WRITE(IWRITE,1003) 1 ((CF(IJCHAN,LAM), LAM=1,LAMAX),IJCHAN=1,ICFS) END IF END IF C C C COLLECT TOGETHER THE TERMS FROM PAIRING OF CONFIGURATIONS C WITHIN THE SETS KIS,KJS C C ZEROISE HSMAT ARRAY C DO 1 IHS=1,IHDIM3 HSMAT(IHS)=ZERO 1 CONTINUE C C LOOP OVER CONFIGURATIONS IN SET KIS C IF THESE SETS CANNOT INTERACT PERFORM THE LOOPS ONLY C ONCE TO SET UP THE CHANNEL POSITIONS C FIRST SET THE UPPER LIMITS TO THE CONFIGURATION LOOPS C IF (LAMSIN .EQ. 0) THEN KICSUP=NCSET(KIS) KJCSUP=NCSET(KJS) ELSE KICSUP=1 KJCSUP=1 END IF C DO 6 KICS=1,KICSUP IF(KIS .EQ. KJS) THEN JCSLO=KICS ELSE JCSLO=1 END IF C C LOOP OVER CONFIGURATIONS IN SET KJS C DO 7 KJCS=JCSLO,KJCSUP IF (LAMSIN .EQ. 0) THEN C C THE CASE OF CONFIGURATIONS HAVING THE SAME SHELL STRUCTURE C IS TREATED SEPARATELY C IF (IRHSGL(ICOUNT) .LT. 0) THEN CALL CMATII(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT) ELSE CALL CMATIJ(L,LP,LAMST,LAMFIN,LAMCFG,ZB,ICOUNT,CMAT,IND) C C THE INDICATOR IND IS SET ZERO IF IC,JC DO NOT COUPLE C************* C THIS MAY CAUSE TROUBLE IF ALL PAIRS CANNOT COUPLE C EG. A ZERO SUB-MATRIX C************ C IF(IND .EQ. 0) GO TO 7 END IF C C ADD THE CMAT, MULTIPLIED BY ITS APPROPRIATE COEFFICIENT C AIJ INTO THE RELEVANT CHANNEL POSITIONS IN THE DIRECT C PART OF THE HAMILTONIAN MATRIX DMAT C WHEN IS .GT. JS THE SETS HAVE BEEN INTERCHANGED. C THE CONFIGURATIONS MUST THEN BE CHANGED BACK TO FIND AIJ C IF(IS .GT. JS) THEN ICS=KJCS JCS=KICS ELSE ICS=KICS JCS=KJCS END IF C END IF C C INITIALISE STORAGE PARAMETERS C IHS=1 ICHANN=ICHAN1 C C LOOP OVER ALL TARGET STATES COMPOSED FROM SET IS C DO 8 IST=1,NSTAT(IS) IF (LAMSIN .EQ. 0) THEN C C FIND THE STATE NUMBER AND THUS FIND THE COEFFICIENT AIJ C ISTAT=NSETST(IS,IST) AI=AIJ(ISTAT,ICS) C WHEN SETS IS, JS ARE THE SAME A PAIR OF CONFIGURATIONS C WILL CONTRIBUTE TWICE TO COMPENSATE FOR THE SAVING MADE C AS A RESULT OF SYMMETRY. JCSLO=KICS IN DO LOOP 7 C IF(IS .EQ. JS .AND. ICS .NE. JCS) THEN AI2=AIJ(ISTAT,JCS) END IF C END IF C C MODIFY LIMITS FOR LOOPING OVER STATES FROM SET JS WHEN IS=JS C AND L=LP. I.E. ICHAN1=JCHAN1 C IF(ICHAN1 .EQ. JCHAN1) THEN JSTLO=IST JCHANN=ICHANN ELSE JSTLO=1 JCHANN=JCHAN1 END IF C C LOOP OVER STATES COMPOSED FROM THE SET JS C DO 9 JST=JSTLO,NSTAT(JS) IF (LAMSIN .EQ. -999) GO TO 91 JSTAT=NSETST(JS,JST) AJ=AIJ(JSTAT,JCS) AIAJ=AI*AJ IF(IS .EQ. JS .AND. ICS .NE. JCS) THEN AIAJ=AIAJ+AI2*AIJ(JSTAT,ICS) END IF C IF(ABS(AIAJ) .LE. EPS) THEN IHS=IHS+NRASQR GO TO 91 END IF C C STORE MATRIX AIAJ*CMAT IN HSMAT C C IF THE CHANNELS ARE REVERSED WHEN CONVERTED TO THE RMATRIX C CONVENTION THE MATRIX MUST BE TRANSPOSED I.E. REVERSE THE C INDICES C IF (LP .NE. L .AND. 1 ICONCH(JCHANN) .LT. ICONCH(ICHANN)) THEN NUP=NRA2 NPUP=NRA1 ELSE NUP=NRA1 NPUP=NRA2 END IF C DO 10 N=1,NUP C DO 11 NP=1,NPUP C C WHEN THE CONTINUUM ANGULAR MOMENTA ARE EQUAL CMAT MUST C BE FILLED UP TO A SQUARE C IF(LP .EQ. L .AND. NP .LT. N) THEN CM=CMAT(N,NP) C C IF THE CHANNELS ARE REVERSED WHEN CONVERTED TO THE C RMATRX CONVENTION THE MATRIX MUST BE TRANSPOSED C ELSE IF (LP .NE. L .AND. 1 ICONCH(JCHANN) .LT. ICONCH(ICHANN)) THEN CM=CMAT(N,NP) C ELSE CM=CMAT(NP,N) END IF HSMAT(IHS)=HSMAT(IHS)+CM*AIAJ IHS=IHS+1 11 CONTINUE C 10 CONTINUE C C 91 JCHANN=JCHANN+1 9 CONTINUE C ICHANN=ICHANN+1 8 CONTINUE IF (NBUG6 .EQ. 1) THEN C if (iam.eq.0)WRITE(IWRITE,1002)L,LP,IS,JS,IC,JC,IST,JST if (iam.eq.0)WRITE(IWRITE,1002)L,LP,IS,JS,ICS,JCS,NSTAT(IS), + NSTAT(JS) if (iam.eq.0)WRITE(IWRITE,1003)(HSMAT(IHSM),IHSM=1,IHS-1) END IF C 7 CONTINUE C 6 CONTINUE IF (NBUG6 .EQ. 1) THEN if (iam.eq.0)WRITE(IWRITE,1002)L,LP,IS,JS if (iam.eq.0)WRITE(IWRITE,1003)(HSMAT(IHSM),IHSM=1,IHS-1) END IF C C TEST THAT THE VSH ARRAY HAS BEEN COMPLETELY USED FOR THIS C SET COMBINATION IF THE SETS INTERACTED C IF (LAMSIN .EQ. 0) THEN IF (ICOUNT .NE. NVSHEL + 1) THEN if (iam.eq.0) WRITE (IWRITE,3000) KIS,KJS STOP END IF END IF C C TRANSFER HSMAT TO A BUFFER MATRIX BY MATRIX C INCLUDE ASYMPTOTIC COEFFICIENTS FOR EACH LAMBDA C AT THE BEGINNING OF THE BUFFER FOR EACH CHANNEL C COMBINATION C ICHAN2=ICHANN-1 JCHAN2=JCHANN-1 IHS=0 ICFS=0 C C SET BUFFER NUMBER FOR DOUBLE BUFFERING C NBUF=1 C DO 12 ICH=ICHAN1,ICHAN2 C C CONVERT CHANNEL NUMBER TO RMATRIX CONVENTION C ICHCON=ICONCH(ICH) C IF (ICHAN1 .EQ. JCHAN1) THEN JCHLO=ICH ELSE JCHLO=JCHAN1 END IF C DO 13 JCH=JCHLO,JCHAN2 C IF(JCH .EQ. ICH) THEN C C ADD IN THE ONE-ELECTRON RK INTEGRAL AND THEN C ADD THE STATE ENERGY INTO THE DIAGONAL ELEMENTS C IHS1=IHS DO 131 IR=1,NRASQR IHS1 =IHS1+1 HSMAT(IHS1)=HSMAT(IHS1)+R1(IR) 131 CONTINUE C IHS1=IHS+1 IST=ISTCH(ICH) E=ENAT(IST) DO 132 IR=1,NRA1 HSMAT(IHS1)=HSMAT(IHS1)+E IHS1 = IHS1+NRA1+1 132 CONTINUE C END IF C C CALCULATE THE RECORD NUMBER SO THAT THE MATRICES ARE C STORED IN SEQUENCE IN THE RMATRIX CONVENTION C JCHCON=ICONCH(JCH) IF(JCHCON .LT. ICHCON) THEN IJREC=NSREC+NCH*(JCHCON-1)-(JCHCON*(JCHCON-1))/2+ICHCON ELSE IJREC=NSREC+NCH*(ICHCON-1)-(ICHCON*(ICHCON-1))/2+JCHCON END IF C IF (NBUF .EQ. 1) THEN C C FILL BUF1D C ICFS=ICFS+1 DO 141 IB=1,LAMAX BUF1D(IB)=CF(ICFS,IB) 141 CONTINUE C DO 142 IB=1+LAMAX,NRASQR+LAMAX IHS=IHS+1 BUF1D(IB)=HSMAT(IHS) 142 CONTINUE C WRITE(JBUFD,REC=IJREC)BUF1D NBUF=2 C ELSE C C FILL BUF2D C ICFS=ICFS+1 DO 143 IB=1,LAMAX BUF2D(IB)=CF(ICFS,IB) 143 CONTINUE DO 144 IB=1+LAMAX,NRASQR+LAMAX IHS=IHS+1 BUF2D(IB)=HSMAT(IHS) 144 CONTINUE C WRITE(JBUFD,REC=IJREC) BUF2D NBUF=1 C END IF C 13 CONTINUE C 12 CONTINUE C 5 CONTINUE C 4 CONTINUE C 3 CONTINUE C 2 CONTINUE C RETURN END C*********************************************************************** SUBROUTINE WRITE1 C C WRITE SECTION OF ASYMPTOTIC FILE WHICH IS INDEPENDENT OF C LRGL, SPIN AND PARITY C C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(IPDIM1=MZSPN) PARAMETER(ISDIM1=MZSET,ISDIM2=MZNCS) PARAMETER(ITDIM1=MZTAR,ITDIM2=MZNSS) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) COMMON/STG1 /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2) COMMON/CPBHEADER/LRANG2ALL,MINIMALL,LR2 COMMON /CPBCOEFFALL/COEFFALL(3,IDIM2) DATA IFLAG4/0/ SAVE IFLAG4 c **** parallel **** include 'mpif.h' common/parablock/iam,nproc c **** parallel **** if (iflag4.ne.0) return if (iam.ne.0) return !! parallel C C WRITE(ITAPE3) NELC,NZ,LRANG2,LAMAX,NAST,RA,BSTO,LRANG1, C 1 (MAXNHF(I),I=1,LRANG2),(MAXNC(I),I=1,LRANG1),NRANG2 WRITE(ITAPE3) NELC,NZ,LRANG2ALL,LAMAX,NAST,RA,BSTO C C**** C NRANG2 SHOULD BE REPLACED BY NOTERM C WRITE(ITAPE3)(ENAT(I),I=1,NAST) WRITE(ITAPE3)(LL(I),I=1,NAST) WRITE(ITAPE3)(LSPN(I),I=1,NAST) C WRITE(ITAPE3)(LPTY(I),I=1,NAST) WRITE(ITAPE3)((COEFFALL(I,L),I=1,3),L=1,LRANG2ALL) C iflag4=iflag4+1 RETURN END c **** parallel **** c---------------------------------------------------------------------- c ************ subroutines for parallel code ***************** c*********************************************************************** subroutine read3p implicit real*8(a-h,o-z) c.......read parameters for the Diagonalization routine include 'mpif.h' character*1 filec,filed,fileu logical ex COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint common/crashblock/istop common/parablock/iam,nproc common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc common /cases/more,mskip,ipwinit,ipwfinal,inast namelist/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest ireadp = 37 !! check if unit 37 is allowed inquire (file='dstg3',exist=ex) if (.not.ex) then if (iam.eq.0) write (iwrite,7777) if (iam.eq.0) write (0,7777) istop=1 go to 8888 endif 7777 format (/' dstg3 does not exist....stopping.') open(unit=ireadp,file='dstg3',status='old',form='formatted') itype=1 irsrc=0 icsrc=0 MB=16 NB=16 NPROW=-1 NPCOL=-1 read(ireadp,matrixdat,end=7778) 7778 IF(NPROW.LE.0.OR.NPCOL.LE.0)THEN T=NPROC NPROW=NINT(SQRT(T)) NPCOL=NPROW ENDIF if (nprow*npcol.ne.nproc) then if (iam.eq.0) then write(0,*)'USAGE: nprow*npcol must be equal to # of processors' write(iwrite,*) + 'USAGE: nprow*npcol must be equal to # of processors' endif STOP c istop=1 c return endif if (MB.ne.NB) MB=NB c....... open timeXX.dat file (XX is the number of processors) if(iam.eq.0)then ipwinit=1 !NX fixed ic = nproc/100 id = (nproc - 100*ic)/10 iu = (nproc - 100*ic - 10*id) ich0 = ichar('0') ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) if (ipwinit.eq.1) then open(unit=33,file='time'//filec//filed//fileu//'.dat', + status='unknown') write(33,2500) nproc,NB else open(unit=33,file='time'//filec//filed//fileu//'.dat', + status='unknown',position='append') write(33,2505) endif 2500 format(15x,'Timing'/3x,i3,' processors - (block size=:',i3,')'/) 2505 format(/5x,'New Run'/) endif c **** parallel **** 8888 return end c*********************************************************************** subroutine initctxt implicit real*8(a-h,o-z) c.......Initialize a single BLACS context include 'mpif.h' common/procdat/ictxt,myrow,mycol common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc c.......Initialize a single BLACS context CALL BLACS_GET( -1, 0, ictxt ) CALL BLACS_GRIDINIT( ictxt, 'R', NPROW, NPCOL ) CALL BLACS_GRIDINFO( ictxt, NPROW, NPCOL, MYROW, MYCOL ) return end c*********************************************************************** subroutine Hsize(NN,irowsize,icolsize) implicit real*8(a-h,o-z) c.......calculate the size of the local Hamiltonian include 'mpif.h' common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc common/procdat/ictxt,myrow,mycol c.......dimension of the local matrix HMAT irowsize = numroc(NN,nb,myrow,irsrc,nprow) icolsize = numroc(NN,mb,mycol,icsrc,npcol) return end ********************************************************************** SUBROUTINE DIAGHMAT(HMAT,irowsize,icolsize,ISP) c-------------------------------------------------------------- c sets up and diagonalizes the hamiltonian matrix c and determines the surface amplitudes from the eigenvectors. c-------------------------------------------------------------- c **** parallel **** IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' c.....variables needed: l2pspn, nvary, ends, wmat, lrgl, and nspval PARAMETER(ITDIM1=MZTAR,ITDIM3=MZCHF,ITDIM2=MZNSS,ITDIM6=MZCHS) PARAMETER(ISDIM1=MZSET) PARAMETER(ILDIM1=MZLR3,ILDIM2=ILDIM1+ILDIM1-1,ILDIM3=MZLR4) PARAMETER(IPDIM1=MZSPN) COMMON/CHANN /NnCHAN,ISTCH(ITDIM3),NCHSET(ISDIM1,ILDIM2), 1 ICONCH(ITDIM3),KCONCH(ITDIM3), 2 L2PSPN(ITDIM6,IPDIM1),KCONAT(ITDIM1,IPDIM1), 3 NSCHAN(IPDIM1),NSPREC(IPDIM1) PARAMETER(IDIM1=MZLR1,IDIM2=ILDIM3+ILDIM1-1,IDIM3=MZNR2) COMMON/RAD /RA,BSTO,LRANG1,LRANG2,MAXNHF(IDIM2),MAXNC(IDIM1), 1 NRANG2,NVARY(IDIM2),LAMAX,NELC,NZ,LLPDIM COMMON/STG1 /COEFF(3,IDIM2),ENDS(IDIM3,IDIM2) PARAMETER(IHDIM1=ITDIM6*IDIM3) DIMENSION WMAT(ITDIM6,IHDIM1) COMMON/SETL /LRGL,NPTY,LTOT,LVAL(ILDIM2),NSCOL(ILDIM2), 1 NSETL(ISDIM1,ILDIM2) COMMON/STATE1/NAST,NSTAT(ISDIM1),NSETST(ISDIM1,ITDIM2), 1 LL(ITDIM1),LSPN(ITDIM1),LPTY(ITDIM1), 2 NSP,NSPVAL(IPDIM1),KSPIN(IPDIM1),LENHAM(IPDIM1) COMMON/INFORM/IREAD,IWRITE,IPUNCH,iprint COMMON/DISC /JBUFF1,JBUFF2,JBUFIR,JBUFFV,JBUFFZ,JBUFD, 1 ITAPE1,ITAPE3 c......for debugger needed ground=enat(1) COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS PARAMETER(ISDIM2=MZNCS) COMMON/NBUG /NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8, 1 NBUG9 COMMON/STATE2/AIJ(ITDIM1,ISDIM2),ENAT(ITDIM1) CHARACTER*4 PARITY(0:1) DATA PARITY/'EVEN','ODD'/ c **** parallel **** include 'mpif.h' common/crashblock/istop common/parablock/iam,nproc allocatable :: value(:),vector(: , :) real*8,allocatable :: wmatarray(:) integer :: status(MPI_STATUS_SIZE) integer, allocatable :: npointer(:,:) integer :: kcount, nksize(0:nproc-1) integer :: itagwmat,itagpoint1,itagpoint2 parameter (itagwmat=42,itagpoint1=43,itagpoint2=44) common /matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc,nproctest dimension HMAT(irowsize,icolsize) DIMENSION wmatl(ITDIM6,IHDIM1) C c....... initialization of surface amplitudes do igv=1,ihdim1 do ich=1,itdim6 c wmat(ich,igv)=0.0D0 wmatl(ich,igv)=0.0D0 enddo enddo c.......size of the hamiltonian matrix, channels, etc. NY = LENHAM(ISP) nchan = NSCHAN(ISP) c.......allocate eigenvalues and eigenvectors arrays c....... eigenvalues are global, eigenvectors local allocate(value(NY),vector(irowsize,icolsize),stat=ierr) if (ierr.ne.0) then print*,'program could not allocate space for eigenvectors' stop endif if (iam.eq.0) then if(iprint.ge.0)write(0,*) ' Symmetry:',LRGL,' ',PARITY(NPTY), + ' size=:',NY write(33,*) ' Symmetry:',LRGL,' ',PARITY(NPTY), + ' size=:',NY endif call cpu_time(timei) c DO NOT use PDSYEVX, does NOT ensure orthogonality of e-vectors across processors. c call DIAG_PDSYEVX(HMAT,irowsize,icolsize,vector,value,NY) c Use PDSYEVD call DIAG_PDSYEVD(HMAT,irowsize,icolsize,vector,value,NY) call cpu_time(timef) if (iam.eq.0) then time = timef-timei if(iprint.ge.0)write(0,2999) time/60. write(33,2999) time/60. 2999 format(2x,'digonalization time=: ',1pg10.3,' min.') endif IF (NBUG9 .GT. 0) THEN GROUND=ENAT(1) if (iam.eq.0) WRITE(IWRITE,'(/'' EIGENVALUES AND SURFACE'' 1 ,'' AMPLITUDES FOR TARGET SPIN (2S+1) ='',I3)') 2 NSPVAL(ISP) do NO=NY,1,-1 ENLEVL=(VALUE(NO)-GROUND)*TWO if (iam.eq.0) WRITE (IWRITE,1014)VALUE(NO),ENLEVL C if (iam.eq.0) WRITE (0,1014)VALUE(NO),ENLEVL enddo END IF c ... initialization of local surface amplitudes c allocate(wmatarray(irowsize*icolsize), * npointer(irowsize*icolsize,2),stat=ierr) if (ierr.ne.0) stop * 'allocation fails for wmatl, wmatarray and npointer' do I=1,ny do J=1,NCHAN wmatl(J,I)=0.0d0 enddo enddo do I=1,irowsize*icolsize wmatarray(I)=0.0d0 npointer(I,1)=0 npointer(I,2)=0 enddo c.......determine the surface amplitudes, wmatl(channel,eigenvalue) c....... wmat(channel,eigvn) = sum{ vector(icont,eigvn)*ends(icont,channel)} isizemax = 0 do 197 ixch=1,nchan lp = l2pspn(ixch,isp)+1 isizemax = isizemax + nvary(lp) 197 continue do 205 irow=1,irowsize do 200 icol=1,icolsize c............... identifying irow,icol --> I,J call demapping(irow,icol,I,J) c................identifying the coupled channel number JCHAN jmx=0 jmn=0 do 198 ixch=1,nchan lp = l2pspn(ixch,isp)+1 noterm = nvary(lp) jmn = jmx jmx = jmx + noterm if (I.le.jmx.and.I.gt.jmn) then jchan = ixch jchinit = jmx - noterm icont = I - jchinit go to 199 endif 198 continue go to 200 c............... wmat is constructed only from isizemax terms 199 if (I.gt.isizemax) go to 200 wmatl(jchan,J) = wmatl(jchan,J) + + vector(irow,icol)*ends(icont,lp) 200 continue 205 continue c c new gathering of surface amplitudes c c.......Store non-zero elements, and their co-ords, in buffers to be sent to c.......the first processor kcount=0 do J=1,ny do jchan=1,nchan if(wmatl(jchan,J).NE.(0.0d0))then kcount=kcount+1 wmatarray(KCOUNT) = wmatl(jchan,J) npointer(KCOUNT,1)=jchan npointer(KCOUNT,2)=J endif enddo enddo 7778 deallocate(vector) c....... initialization of global surface amplitudes if (nproctest.ne.0) go to 4000 c....... Sum the local wmatl to the global wmat: c....... First, gather all values of kcount to the first processor CALL MPI_GATHER(KCOUNT,1,MPI_INTEGER,nksize(0),1,MPI_INTEGER,0, * MPI_COMM_WORLD,ierr) c....... Send all buffers to the first processor. Must make sure the c....... mpi tag on messages are different, especially the 2 integer sends. c....... Tag values are set in parameter statement above. if (iam.ne.0) then CALL MPI_SEND(wmatarray,KCOUNT,MPI_REAL8,0,itagwmat, * MPI_COMM_WORLD,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr CALL MPI_SEND(npointer(1,1),KCOUNT,MPI_INTEGER,0,itagpoint1, * MPI_COMM_WORLD,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr CALL MPI_SEND(npointer(1,2),KCOUNT,MPI_INTEGER,0,itagpoint2, * MPI_COMM_WORLD,ierr) if (ierr.ne.0) write(0,*)'MPI_SEND: iam, ierr = ', iam, ierr else c....... Receive the buffers from other processors, and add to global wmat do J = 1, nproc-1 CALL MPI_RECV(wmatarray,nksize(J),MPI_REAL8,J,itagwmat, * MPI_COMM_WORLD,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr CALL MPI_RECV(npointer(1,1),nksize(J),MPI_INTEGER,J, * itagpoint1,MPI_COMM_WORLD,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr CALL MPI_RECV(npointer(1,2),nksize(J),MPI_INTEGER,J, * itagpoint2,MPI_COMM_WORLD,status,ierr) if (ierr.ne.0) write(0,*)'MPI_RECV: iam, ierr = ', iam, ierr do L=1,nksize(J) wmatl(npointer(L,1),npointer(L,2)) = * wmatl(npointer(L,1),npointer(L,2)) + * wmatarray(L) enddo enddo endif deallocate(wmatarray,npointer) c c c....... sum the local wmatl to the global wmat c call mpi_reduce(wmatl,wmat,itdim6*ihdim1,mpi_real8,mpi_sum,0, c + mpi_comm_world,ierr) if (iam.ne.0) go to 4000 c.......write eigenvalues and surface amplitudes to itape3 write (itape3) (value(i),i=NY,1,-1) write (itape3) ((wmatl(k,i),k=1,nchan),i=NY,1,-1) IF (NBUG9 .GT. 0) THEN if (iam.eq.0) WRITE (IWRITE,1003) if (iam.eq.0) WRITE (IWRITE,1005) + ((wmatl(k,i),k=1,nchan),i=NY,1,-1) END IF C WRITE(IWRITE,'(//I7,'': '',''HAMILTONIAN DIAGONALISED FOR L='', 1 I3,'', PARITY='',A4,'', TARGET SPIN (2S+1)='',I2//)') 2 NY,LRGL,PARITY(NPTY),NSPVAL(ISP) C 1003 FORMAT (' SURFACE AMPLITUDE') 1005 FORMAT( 2X,8F14.7) 1014 FORMAT (/' EIGENVALUE=',F12.7,' RYD - RELATIVE TO TARGET GROUND', 1 F12.7,' AU') c........ VERY IMPORTANT !! (leave place for next symmetry) 4000 deallocate(value) c deallocate(vector) c deallocate(wmatl) c if(iam.eq.0)call flush_(iwrite) c RETURN END c*********************************************************************** subroutine diag_pdsyevx(HNP1,irowsize,icolsize,Z,W,NN) implicit real*8(a-h,o-z) c.......Matrix diagonalization for parallel machine c.......driver for the SCAlapack subroutine pdsyevx include 'mpif.h' common/parablock/iam,nproc common/procdat/ictxt,myrow,mycol dimension HNP1(irowsize,icolsize) dimension Z(irowsize,icolsize) ! local eigenvectors dimension W(NN) ! global eigenvalues c.......parameters and variables used in pdsyevx parameter (rminusone=-1.0d0,zero=0.0d0) integer desca(9),descz(9) integer ifail,iwork,iclustr allocatable :: ifail(:),iwork(:),iclustr(:) allocatable :: work(:),gap(:) integer ione,itwo,ifour,ifive,isix data ione,itwo,ifour,ifive,isix/1,2,4,5,6/ common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc c.......other parameters needed in pssyevx abstol = rminusone ! error tolerance of eigenvalues orfac = rminusone ! reorthogonalization lda = max(irowsize,ione) allocate(ifail(NN),stat=ierr) if(ierr.ne.0)stop 'allocation fails for ifail' c.......space required to guarantee that all eigenvectors are computed neig = NN NNN = max(NN,nb,itwo) np0 = numroc(NNN,nb,0,0,nprow) mq0 = numroc(max(neig,nb,itwo),nb,0,0,npcol) lwork = ifive*NN + max(ifive*NNN,np0*mq0+itwo*nb*nb) + + iceil(neig,nprow*npcol)*NNN allocate(work(lwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for work' c.......working array liwork = isix*(max(NN,nprow*npcol + ione, ifour)) allocate(iwork(liwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for iwork' c.......gap and iclustr allocate(iclustr(nproc*2),stat=ierr) if (ierr.ne.0)stop 'allocation fails for iclustr' allocate(gap(nproc),stat=ierr) if (ierr.ne.0)stop 'allocation fails for gap' c.......basic array descriptors call descinit(desca,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call descinit(descz,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call mpi_barrier(mpi_comm_world,ierr) c......compute the entire eigendecomposition in PDSYEVX c pdsyevx for SUN, SGI and IBM --- pssyevx for CRAY call pdsyevx('V','A','U',NN,HNP1,1,1,desca,zero,zero,13, + -13,abstol,m,nz,W,orfac,Z,1,1,descz,work, + lwork,iwork,liwork,ifail,iclustr,gap,info) c if(info.ne.0) then c if (iam.eq.0) write(0,*) ' pdsyevx error - info=:',info c stop 'pdsyevx error' c endif deallocate(ifail) deallocate(work) deallocate(iwork) deallocate(gap) deallocate(iclustr) return end c*********************************************************************** subroutine diag_pdsyevd(HNP1,irowsize,icolsize,Z,W,NN) implicit real*8(a-h,o-z) c.......Matrix diagonalization for parallel machine c.......driver for the SCAlapack subroutine pdsyevd include 'mpif.h' common/parablock/iam,nproc common/procdat/ictxt,myrow,mycol dimension HNP1(irowsize,icolsize) dimension Z(irowsize,icolsize) ! local eigenvectors dimension W(NN) ! global eigenvalues c.......parameters and variables used in pssyevd integer desca(9),descz(9) integer iwork allocatable :: iwork(:) allocatable :: work(:) integer ione,itwo,ithree,isix,iseven,ieight data ione,itwo,ithree,isix,iseven,ieight/1,2,3,6,7,8/ common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc c.......other parameters needed in pdsyevd lda = max(irowsize,ione) c.......space required to guarantee that all eigenvectors are computed np0 = numroc(NN,nb,myrow,irsrc,nprow) mq0 = numroc(NN,nb,mycol,icsrc,npcol) itril = ithree*NN + max(nb*(np0+ione),ithree*nb) lwork = itwo*NN + max(ione+isix*NN+itwo*np0*mq0,itril) allocate(work(lwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for work' c.......working array liwork = iseven*NN + ieight*npcol + itwo allocate(iwork(liwork),stat=ierr) if(ierr.ne.0)stop 'allocation fails for iwork' c.......basic array descriptors call descinit(desca,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call descinit(descz,NN,NN,nb,mb,irsrc,icsrc,ictxt,lda,info) call mpi_barrier(mpi_comm_world,ierr) c......compute the entire eigendecomposition in PDSYEVD c pdsyevd SUN, SGI and IBM, pssyevd CRAY call pdsyevd('V','U',NN,HNP1,1,1,desca,W,Z,1,1,descz,work, + lwork,iwork,liwork,info) c if(info.ne.0) then c if (iam.eq.0) write(0,*) ' pdsyevd error - info=:',info c stop 'pdsyevd error' c endif deallocate(work) deallocate(iwork) return end c*********************************************************************** subroutine demapping(irow,icol,I,J) implicit real*8(a-h,o-z) c.......demapping the local (irow,icol) in the global (I,J) position common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc common/procdat/ictxt,myrow,mycol ix = mod(irow,nb) if (ix.eq.0) ix=nb iy = mod(icol,mb) if (iy.eq.0) iy=mb c....... l and m: local block coordinate l = (irow-ix)/nb m = (icol-iy)/mb I = (l*nprow + myrow) * nb + ix J = (m*npcol + mycol) * mb + iy return end c*********************************************************************** subroutine locH(NN,ijgbindex,irow,icol,ipr,ipc) implicit real*8(a-h,o-z) c.......Mapping the global matrix into the local matrix c....... IJGBINDEX ----> (I,J) ------> (irow,icol) include 'PARAM' common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc common/procdat/ictxt,myrow,mycol c........identifying I,J from ijgbindex call ntoij(NN,ijgbindex,i,j) c........mapping the GLOBAL matrix A(I,J) into LOCAL A(irow,icol) c........ipr and ipc: row and column processes call mapping(I,J,irow,icol,ipr,ipc) return end c*********************************************************************** subroutine ntoij(NN,IJ,i,j) implicit real*8(a-h,o-z) parameter(two=2.0d0) c........identifying i,j from index IJ b = (2*NN+3)/two bb=b*b q=two*(NN+IJ) i = b - sqrt(bb-q) j = NN + IJ - i*(2*NN - i + 1)/two return end c*********************************************************************** subroutine ijton(NN,IJ,i,j) c.......identifying index IJ from i,j IJ = (NN*(NN+1))/2-((NN-i+1)*(NN-i+2))/2+j-i+1 return end c*********************************************************************** subroutine mapping(I,J,irow,icol,ipr,ipc) c.......mapping the GLOBAL matrix A(I,J) into LOCAL A(irow,icol) include 'PARAM' common/matrixdat/MB,NB,nprow,npcol,itype,irsrc,icsrc common/procdat/ictxt,myrow,mycol c..........ipr and ipc: row and column processes ipr = mod(irsrc + (I-1)/MB , nprow) ipc = mod(icsrc + (J-1)/NB , npcol) c.......column identification (J --> icol) mblock = (J-1)/(npcol*NB) iypos = mod(J-1,NB) + 1 icol = mblock*NB + iypos c.......row identification (I --> irow) lblock = (I-1)/(nprow*MB) ixpos = mod(I-1,MB) + 1 irow = lblock*MB + ixpos return end c*********************************************************************** SUBROUTINE HEADER C INCLUDE 'PARAM' C implicit real*8(a-h,o-z) CHARACTER*1 NUM(0:9) CHARACTER*11 abcd LOGICAL EX,exs DATA IFLAGHEAD/0/ SAVE IFLAGHEAD c c DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ c NAMELIST/STGNX/MINLT,MAXLT,LRGLLO,LRGLUP,MAXC,LFIXN,MJS,NPTYIN X ,RELOP,IWAVE X ,NBUG1,NBUG2,NBUG3,NBUG4,NBUG5,NBUG6,NBUG7,NBUG8,NBUG9 c include 'mpif.h' common/parablock/iam,nproc common/CPBCOEFF/LRANG2UP(1128),NNFILE c IF(IFLAGHEAD.ne.0)return c DO I=0,MZLR3*MZLR4 ! HIGHEST SPIN* HIGHEST ANGULAR MOMENTUM i1=I/100 i2=(I-100*(I/100))/10 i3=I-(100*(I/100))-i2*10 abcd='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) inquire (file=abcd,exist=ex) if(.NOT.ex)then II=I goto 2222 endif ENDDO c c 2222 continue c inquire (file='sizeNX.dat',exist=exs) if(exs)then OPEN(UNIT=2228,FILE='sizeNX.dat',STATUS='UNKNOWN') REWIND(2228) DO I=1,II READ(2228,*) NDUM,LRANG2UP(I) ENDDO NNFILE=II else inquire (file='RAD1.DAT',exist=exs) if(exs)then OPEN(UNIT=2228,FILE='dstgnx',STATUS='UNKNOWN') READ(2228,*) !Continuation line READ(2228,STGNX) NNFILE=1 LRANG2UP(1)=MAXLT else WRITE(6,*)' no sizeNX.dat,RAD1.DAT or RAD1.DATXXX' stop ' no sizeNX.dat,RAD1.DAT or RAD1.DATXXX' endif endif c close(2228) iflaghead=1 return end c*********************************************************************** SUBROUTINE BUTTLECOEFF C C READS BASIC INFORMATION FROM THE FILE RAD3 C IMPLICIT REAL*8 (A-H,O-Z) C INCLUDE 'PARAM' C PARAMETER(IDIM1=MZLR1,IDIM3=MZNR2) PARAMETER(ILDIM1=MZLR3,ILDIM3=MZLR4) PARAMETER(IDIM2=ILDIM3+ILDIM1-1) C CHARACTER*1 NUM(0:9) CHARACTER*11 abcd C LOGICAL ex C COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS COMMON /CPBCOEFF/LRANG2UP(1128),NNFILE COMMON /CPBHEADER/LRANG2ALL,MINIMALL,LR2 COMMON /CPBCOEFFALL/COEFFALL(3,IDIM2) C C DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ DATA IFLAGBUTT/0/ SAVE IFLAGBUTT C DIMENSION NVARY(IDIM2),COEFFD(3,IDIM2),ENDS(IDIM3,IDIM2) C if(IFLAGBUTT.ne.0) return C DO I=1,IDIM2 DO J=1,3 COEFFALL(J,I)=ZERO ENDDO ENDDO DO I=1,NNFILE LRANG2UP(I)=LRANG2UP(I)+LR2 ENDDO C DO III=0,NNFILE-1,1 K=III+1 i1=III/100 i2=(III-100*(III/100))/10 i3=III-(100*(III/100))-i2*10 abcd='RAD3.DAT'//NUM(i1)//NUM(i2)//NUM(i3) inquire (file=abcd,exist=ex) if(ex)then OPEN(UNIT=2228,FILE=abcd,ACCESS='SEQUENTIAL',STATUS='OLD', 1 FORM='UNFORMATTED') else OPEN(UNIT=2228,FILE='RAD3.DAT',ACCESS='SEQUENTIAL', 1 STATUS='OLD',FORM='UNFORMATTED') endif C REWIND(2228) !old ITAPE=ITAPE1 READ(2228) RA,BSTO READ(2228) (NVARY(L),L=1,LRANG2UP(K)) READ(2228) ((ENDS(N,L),N=1,NVARY(L)),L=1,LRANG2UP(K)) READ(2228) ((COEFFD(I,L),I=1,3),L=1,LRANG2UP(K)) C DO J=1,LRANG2UP(K) IF(COEFFD(3,J).ne.ZERO)THEN COEFFALL(1,J)=COEFFD(1,J) COEFFALL(2,J)=COEFFD(2,J) COEFFALL(3,J)=COEFFD(3,J) ENDIF ENDDO CLOSE(2228) ENDDO C IFLAGBUTT=1 RETURN END c **** parallel ****