program grasprad C C Written by T G Lee and C P Ballance C C CPB 21/2/11 .... follow the steps below C C Warning : you may be tempted to capitalise this code C but you will break the search strings !!! C C Extract the radiative rates from a GRASP.OUT file C Default extraction : Dirac Fock C C other options : Breit, Breit+QED C C ****** In the GRASP.INP, we need 2 runs ******* C C MCT 1 0 ! creates a radmatrixE1M1 with E1, M1 transitions C C MCT 2 0 ! creates a radmatrixE2M2 with E2 , M2 transitions C C Step 1. run GRASP with MCT 1 0 C C Step 2. run grasprad with ie1m1=1 =====> radmatrixE1M1 C C Step 3. run GRASP with MCT 2 0 C C Step 4. run grasprad.f iwth ie1m1=2 =====> radmatrixE2M2 C C In radmatrixE1M1 -E1(positive) , M1 (negative) C In radmatrixE2M2 -E2(positive) , M2 (negative) C C In fact the second run picks up the radmatrixE1M1 file C and creates adf04rad with everything contained. This enables C a simple column 'cut and paste' into a standard adf04 file. C implicit real*8(a-h,o-z) character*4 char4 CHARACTER*1 NUM(0:9) character*37 marker1 character*13 fileaa,fileab,fileac character*11 marker2,marker3 character*8 markercoulomb !will be Breit character*72 char1 character*20 blankline character*3 expon character*4 emant logical EX,EX2 real*8,allocatable :: ENAT(:),TRANSMATRIX(:,:) integer,allocatable :: JLEVEL(:),JPAR(:), X NUMCSF(:),IPOINTER(:,:) DATA NUM/'0','1','2','3','4','5','6','7','8','9'/ C C C NAMELIST/adf04/nosoftemps,noselevs,noswitch,ie1m1,maxline,isub C OPEN(9,file='ELEV.DAT',form='formatted',status='unknown') OPEN(10,file='drad',form='formatted',status='unknown') OPEN(11,file='GRASP.OUT',form='formatted',status='unknown') C IBIGNUMBER=10000 MAXLINE=10000000 C C nosoftemps=0 noselevs=0 noswitch=0 ie1m1=1 isub=0 C READ(10,adf04) C if(isub.eq.0)isub=noselevs C if(noselevs.ne.0)ibignumber=noselevs C allocate(ENAT(IBIGNUMBER)) allocate(TRANSMATRIX(IBIGNUMBER,IBIGNUMBER)) allocate(JLEVEL(IBIGNUMBER)) allocate(JPAR(IBIGNUMBER)) allocate(NUMCSF(IBIGNUMBER)) allocate(IPOINTER(IBIGNUMBER,2)) C fileaa='radmatrixE'//NUM(ie1m1)//'M'//NUM(ie1m1) OPEN(12,file=fileaa,form='formatted',status='unknown') C C Initialise C DO I=1,IBIGNUMBER DO J=1,IBIGNUMBER TRANSMATRIX(J,I)=0.0d0 ENDDO ENAT(I)=0.0d0 JLEVEL(I)=0 JPAR(I)=0 NUMCSF(I)=0 IPOINTER(I,1)=I IPOINTER(I,2)=I ENDDO C do I=1,noswitch READ(10,*)III,IIJ IPOINTER(III,1)=III IPOINTER(III,2)=IIJ enddo C do II=1,MAXLINE read(11,2222)marker1 if(marker1.eq.' eigenenergies relative to the lowest')then write(0,*)' energy levels found at line ',II,' .... reading ', X noselevs do III=1,4 read(11,2223)char1 enddo ! should be at top of energy levels goto 50 endif if(II.eq.MAXLINE)then write(0,*)'over ', X MAXLINE,'lines of input and cannot energy levels' stop endif enddo C 50 continue C C do II=1,noselevs read(11,2224)idum1,idum2,char4,idum3,tmp1,tmp2,tmp3 write(9,2224)idum1,idum2,char4,idum3,tmp1,tmp2,tmp3 NUMCSF(II)=idum3 C C if(char4.eq.'even')then JPAR(II)=0 else JPAR(II)=1 endif C JLEVEL(II)=idum2 C ENAT(II)=tmp3 C write(0,*)NUMCSF(II),JLEVEL(II),JPAR(II),ENAT(II) enddo C C C C do JJ=1,MAXLINE read(11,2225)marker2 if(marker2.eq.' Electric 2')then write(0,*)'FOUND E' do III=1,2 read(11,2223)char1 enddo read(11,2226)markercoulomb write(0,*)'Using Dirac-Fock values' if(markercoulomb.eq.' using C')then ! C==> B will give Breit continue C C start building the TRANSMATRIX C do III=1,8 read(11,2223)char1 enddo C C 75 read(11,2227)blankline C write(0,*)blankline if(blankline.eq.' ')then goto 100 else backspace(11) read(11,3334)idum1,idum2,tmp1,tmp2,tmp3,tmp4 C write(0,3334)idum1,idum2,tmp1,tmp2,tmp3,tmp4 TRANSMATRIX(idum1,idum2)=tmp4 goto 75 endif C 3334 FORMAT(1X,2I5,1P,4E12.4) C else write(0,*)'through the Coulomb transitions' goto 99 endif endif C if(JJ.eq.MAXLINE)then write(0,*)'over',maxline,' lines processed ... somethings up' stop endif C 99 continue enddo C 100 continue do KK=1,MAXLINE read(11,2228)marker3 if(marker3.eq.' Magnetic 2')then write(0,*)'FOUND MAG ' do III=1,11 read(11,2223)char1 enddo C 175 read(11,2227)blankline C write(0,*)blankline if(blankline.eq.' ')then goto 200 else backspace(11) read(11,*)idum1,idum2,tmp1,tmp2,tmp3 C write(0,*)idum1,idum2,tmp1,tmp2,tmp3 TRANSMATRIX(idum1,idum2)=-tmp3 goto 175 endif else continue endif enddo 200 continue C C C Start of the transmatrix dump C do I=1,noselevs do J=1,noselevs write(12,*)IPOINTER(I,2),IPOINTER(J,2), X TRANSMATRIX(IPOINTER(I,2),IPOINTER(J,2)) enddo enddo CLOSE(12) C fileaa='radmatrixE'//NUM(1)//'M'//NUM(1) fileab='radmatrixE'//NUM(2)//'M'//NUM(2) INQUIRE (FILE=fileaa,EXIST=EX) INQUIRE (FILE=fileab,EXIST=EX2) write(0,*)EX,EX2 IF((EX).AND.(EX2))then write(0,*)'both e1m1 and e2m2 files exist ... creating adf04 rad' OPEN(12,file=fileaa,form='formatted',status='unknown') OPEN(13,file=fileab,form='formatted',status='unknown') OPEN(14,file='adf04rad',form='formatted',status='unknown') OPEN(15,file='SCRATCHW',form='formatted',status='unknown') REWIND(12) REWIND(13) do II=1,noselevs do JJ=1,noselevs TRANSMATRIX(JJ,II)=0.0d0 enddo enddo C do II=1,noselevs do JJ=1,noselevs READ(12,*)III,JJJ,tmpa TRANSMATRIX(III,JJJ)=TRANSMATRIX(III,JJJ)+abs(tmpa) READ(13,*)III,JJJ,tmpa TRANSMATRIX(III,JJJ)=TRANSMATRIX(III,JJJ)+abs(tmpa) enddo enddo C do II=1,isub do JJ=II+1,isub tmpaa=TRANSMATRIX(II,JJ) TRANSMATRIX(II,JJ)=max(1.0d-30,tmpaa) rewind(15) write(15,730)TRANSMATRIX(II,JJ) rewind(15) read(15,740)emant,expon WRITE(14,750)JJ,II,emant,expon enddo enddo C 730 format(1pe9.2) 740 format(1x,a4,1x,a3) 750 format(2i5,(1x,a4,a3)) C ENDIF C close(12) close(13) close(14) C 2222 FORMAT(A37) 2223 FORMAT(A72) 2224 FORMAT(I5,1X,I5,4X,A4,I7,2X,F5.3,3X,E16.9,E16.9) 2225 FORMAT(A11) 2226 FORMAT(A8) 2227 FORMAT(A20) 2228 FORMAT(A11) C end