C------------------------------------------------------------------------------- C P I C O N T R - To extract parameter contributions to measured C transition frequencies from parameter derivatives C printed in the .EGY file generated by SPCAT by setting C the EGYFLG in the .int file to 2 C------------------------------------------------------------------------------- C C C Information on transitions to use for contributions comes from the C SPFIT results file as reformatted with PIFORM (usually given the extension .RES) C C Information on parameters of fit comes from the .PAR file for the fit. C C The output file is basically the .RES file with parameter contributions C appended to each transition line. C C C NOTE: It is crucial to ensure that the .PAR, .RES and .EGY files correspond C to the same fit C C Only MAXVIB different vibrational states are allowed. The change C to two digit state identifier is accounted for but derivatives seem to be C picked incorrectly for V>9 C C Some declaration schemes involving dependent parameters could give C unreliable results (needs to be checked). C C Declaration scheme with V2V1=99 at the end of a parameter identifier C also needs checking. C C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C Input/output file names are declared in file PICONTR.INP, which should C consist of four active lines: C c name of the file for output: molnam_contrib.res c name of .par file for the fit: molnam.par c name of .res file from PIFORM: molnam.res c name of .egy file with derivatives: molnam.egy C_______________________________________ C d e s c r i p t o r |_____ d a t a (starting on column 41) C C C C ver. 7.III.2023 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c 21.12.21: creation, with some adaptations from VICONTR c 15.02.22: tweaking c 21.02.22: further development c 2.03.22: modified parameter selection c 24.06.22: debugging input of 3rd line of .PAR c 26.01.23: adding transfer of end of line annotations to the output c 7.03.23: MAXLEV warning C C_____________________________________________________________________________ C C The number of energy levels for J from 0 to n is: n(n+1)+n+1 C MAXLEV needs to account for the potential number of energy levels C C Key counters: C C nlevels - the number of read energy levels with their derivatives C npars - the number of parameters in the .PAR file C ncontr - the number of parameters for which frequency contributions C are to be worked out C C pardesc - alphanumeric parameter descriptors C numpar - numbers of parameters for working out derivatives c program picontr c IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*8 (I-N) parameter (maxlev=300000, maxcon=500, maxvib=20) c character line*300,parnote*2(maxcon),annot*250, * fileout*50,filepar*50,fileres*50,fileegy*50, * pardesc*10(maxcon),descdummy*10 integer*8 indlev(maxlev),indpar(maxcon),iqn(6), * indup,indlw,numpar(11), * ivlow(maxcon),ivup(maxcon) real*8 deriv(maxlev,maxcon),param(maxcon),contrib(11), * deru(11),derl(11) c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Deal with the control file (unit 4) c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c WRITE(*,5500) 5500 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | PICONTR - Contribution extractor from EGY file from SPCAT', * T79,'|'/ * ' | obtained by setting EGYFLG=2 in the .INT file', * T79,'|'/ * ' |',76(1H_),'|'/' version 7.III.2023',T64,'Zbigniew KISIEL'/) c open(4,file='picontr.inp',status='old',err=1027) read(4,'(40x,a)',err=1027)fileout read(4,'(40x,a)',err=1027)filepar read(4,'(40x,a)',err=1027)fileres read(4,'(40x,a)',err=1027)fileegy close(4) c c open(3,file=filepar,status='old',err=1030) ! .par file = unit 3 write(*,'(1x,2a)')'Reading from: ',filepar(1:len_trim(filepar)) c open(1,file=fileres,status='old',err=1029) ! .res file = unit 1 write(*,'(1x,2a)')'Reading from: ',fileres(1:len_trim(fileres)) c open(4,file=fileegy,status='old',err=1028) ! .egy file = unit 4 write(*,'(1x,2a)')'Reading from: ',fileegy(1:len_trim(fileegy)) c open(2,file=fileout,status='unknown') ! .log output = unit 2 write(*,'(1x,2a)')'Writing to: ',fileout(1:len_trim(fileout)) c open(7,file='picontr_identified.out',status='unknown') ! identified levels = unit 7 write(7,50)'List of lines with identified energy levels ' * //'and their derivatives' c open(8,file='picontr_unidentified.out',status='unknown') ! unidentified levels = unit 8 write(8,50)'List of lines with unidentified energy levels' 50 format(75(1H-)/1x,a/75(1H-)//) c c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Read the .PAR file (unit 3) and fish out the parameter names and their c consecutive numbers c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...read the number of vibrational states: c 1/use NFIRST to skip the obsolete alphanumeric descriptor at the beginning of c third .PAR line c 2/use NSHIFT to account for >9 states c read(3,'(a)',end=19,err=19)line read(3,'(a)',end=19,err=19)line read(3,'(a)',end=19,err=19)line c nfirst=2 1143 read(line(nfirst:),*,err=1144)ndummy,nvib if(nvib.lt.0)nvib=-nvib goto 1145 c 1144 nfirst=nfirst+1 if(nfirst.gt.20)goto 1142 goto 1143 1145 write(*,'(1x)') write(*,20)nvib,' vibrational state(s) declared in .PAR' nshift=0 if(nvib.ge.10)nshift=2 c c...skip to the first line with a parameter declaration c 16 read(3,'(a)',end=19,err=19)line if(line(14+nshift:14+nshift).ne.' '.or. * line(33+nshift:33+nshift).ne.'E'.or. * line(49+nshift:49+nshift).ne.'E')goto 16 backspace (3) c c...read the parameter indices c npars=0 nonnot=0 18 read(3,'(a)',err=19,end=19)line read(line(1:13+nshift),*,err=19)index if(index.lt.0)goto 18 ! skip a dependent parameter read(line(14+nshift:37+nshift),*,err=19)parval npars=npars+1 param(npars)=parval indpar(npars)=index c c...read the paramater names (if not available then use the parameter index) c if(len_trim(line).lt.56+nshift)then write(pardesc(npars),'(i10)')indpar(npars) noannot=noannot+1 89 if(pardesc(npars)(1:1).eq.' ')then pardesc(npars)(1:)=pardesc(npars)(2:) goto 89 endif else read(line(56+nshift:len_trim(line)),*,err=19)pardesc(npars) endif parnote(npars)=' ' if(nshift.eq.0)then read(line(12:12),'(I1)')ivup(npars) read(line(13:13),'(I1)')ivlow(npars) c if(ivup(npars).ne.ivlow(npars))parnote(npars)='* ' else read(line(12:13),'(I2)')ivup(npars) read(line(14:15),'(I2)')ivlow(npars) c if(ivup(npars).ne.ivlow(npars))parnote(npars)='* ' endif if(dabs(param(npars)).lt.1.d-25)parnote(npars)='x ' goto 18 c 19 close(3) write(*,20)npars,' declared parameters in .PAR' 20 format(i8,a) c if(nvib.gt.1)then endif c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Read the .EGY file (unit 4): c 1/ fish out energy labels and convert to one index c 2/ fish out derivatives (in EGY six per line) c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c nfull=npars/6 ! full derivative rows nleft=npars-6*nfull ! completing derivative row c write(*,*)nfull,nleft c nlevels=0 read(4,'(a)',err=24,end=24)line nqn=(len_trim(line)-64)/3 write(*,20)nqn,' quantum numbers describing energy levels' rewind(4) c open(10,file='edump',status='unknown') ! debug c c . . . . . c 23 read(4,'(a)',err=24,end=24)line read(line(65:),'(6i3)')(iqn(i),i=1,nqn) nlevels=nlevels+1 c if(nlevels.gt.maxlev)then write(*,1101)nlevels,maxlev stop endif 1101 format(1x//1x,60(1H-)// * ' The number of input energy levels reached:',i7/ * ' while the maximum allowed is MAXLEV=',i7// * ' Increase MAXLEV manually and recompile the program.'//) c index=iqn(1) do 101 n=2,nqn index=index+iqn(n)*1000**(n-1) 101 continue c indlev(nlevels)=index c write(10,'(1x,i20,3x,6i4)')index,(iqn(i),i=1,nqn) ! debug c jstart=1 do 25 n=1,nfull read(4,'(a)',err=24,end=24)line lstart=9 do 26 j=jstart,jstart+5 read(line(lstart:lstart+19),27)jj,deriv(nlevels,j) lstart=lstart+20 26 continue 27 format(i6,e14.3) jstart=jstart+6 25 continue c if(nleft.ne.0)then read(4,'(a)',err=24,end=24)line lstart=9 do 28 j=jstart,npars read(line(lstart:lstart+19),27)jj,deriv(nlevels,j) lstart=lstart+20 28 continue endif c goto 23 c c . . . . . c 24 close(4) write(*,20)nlevels,' energy levels in .EGY' ncontr=0 c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Decide on parameters for which frequency contributions are to be worked out c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c 30 write(*,45) 45 format(1x/ * ' List of parameters declared in the .PAR file ', * '(frequency contributions for '/ * ' up to ten of these can be selected, ', * '0 or just ENTER ends selection):'/) nn=npars/4 write(*,'(a)')' N vv Parname N vv Parname '// * 'N vv Parname N vv Parname' write(*,135) 135 format(' |',19x,'|',19x,'|',19x,'|') do 31 n=1,nn write(*,132) * n, ivup(n), ivlow(n), * parnote(n), pardesc(n), * n+nn, ivup(n+nn), ivlow(n+nn), * parnote(n+nn), pardesc(n+nn), * n+2*nn, ivup(n+2*nn),ivlow(n+2*nn), * parnote(n+2*nn), pardesc(n+2*nn), * n+3*nn, ivup(n+3*nn),ivlow(n+3*nn), * parnote(n+3*nn), pardesc(n+3*nn) c c write(*,32)n, parnote(n), pardesc(n), c * n+nn, parnote(n+nn), pardesc(n+nn), c * n+2*nn, parnote(n+2*nn), pardesc(n+2*nn), c * n+3*nn, parnote(n+3*nn), pardesc(n+3*nn) c 31 continue 32 format(i6,a,a,i8,a,a,i8,a,a,i8,a,a) 132 format(i4,1x,i1,i1,a,a,i5,1x,i1,i1,a,a, * i5,1x,i1,i1,a,a, * i5,1x,i1,i1,a,a) c if(4*nn.ne.npars)then nleft=npars-4*nn do 49 nnn=npars-nleft+1,npars write(*,133)nnn,ivup(nnn),ivlow(nnn), * parnote(nnn),pardesc(nnn) c write(*,33)nnn,parnote(nnn),pardesc(nnn) 49 continue endif 33 format(61x,i5,a,a) 133 format(59x,i5,1x,i1,i1,a,a) c 34 write(*,74) 74 format(1x/' x = zero valued parameter') if(noannot.gt.0)then write(*,77) 77 format(' ------> some alphanumeric ', * 'parameter descriptors are missing in .PAR,'/ * ' numerical parameter identifiers are', * ' used instead') endif c if(ncontr.ge.1)then write(*,70)(numpar(iii),iii=1,ncontr) endif 70 format(1x/ * ' Parameters declared so far: ',11i4) c write(*,72) 72 format(1x/ * ' Type value of N to select parameter for contributions, '/ * ' or just ENTER to terminate .... ',$) read(*,'(i5)',err=30)ipar if(ipar.lt.-(npars-numshift).or.ipar.gt.npars)goto 30 if(ipar.eq.0)goto 71 ncontr=ncontr+1 if(ipar.lt.0)then numpar(ncontr)=-ipar else numpar(ncontr)=ipar endif c if(ncontr.eq.10)then goto 71 else goto 30 endif c 71 write(*,'(1x)') c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Read the .RES file (unit 1): c Write copy of the .RES with appended contributions to the log file (unit 2) c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...deal with the header c read(1,'(a)',err=36,end=36)line write(2,'(a)')line(1:len_trim(line)) read(1,'(a)',err=36,end=36)line write(2,'(a)')line(1:len_trim(line)) do 75 n=1,5 read(1,'(a)',err=36,end=36)line write(2,'(a)')line(1:72) 75 continue c c...right justify parameter identifiers and write derivatives header c do 81 i=1,npars j=len_trim(pardesc(i)) if(j.eq.10)goto 81 descdummy=pardesc(i) pardesc(i)=' ' pardesc(i)(11-j:10)=descdummy(1:j) 81 continue c write(2,80)'Parameter descriptor:',(pardesc(numpar(i)),i=1,ncontr) * ,' SUM' write(2,85)'Parameter number: ',( numpar(i) ,i=1,ncontr) write(2,85)'Parameter index: ',( indpar(numpar(i)),i=1,ncontr) 80 format(a,52x,11(2x,a10)) 85 format(a,52x,11( i12)) write(2,'(1x)') c write(7,82)(pardesc(numpar(i)),i=1,ncontr) write(7,87)( numpar(i) ,i=1,ncontr) write(7,87)( indpar(numpar(i)),i=1,ncontr) 82 format('n_transition level index',2x,11(4x,a10)) 87 format( 34x,11( i14)) write(7,'(1x)') c c- - - top of loop for reading transitions - - - - - - - - - - - - - - - - - - - c ntread=0 nident=0 nunident=0 nout=0 35 read(1,'(a)',err=36,end=36)line c c...not a transition line, so just copy over to the output c if(line(25:25).eq.',')goto 126 if(line(51:51).eq.'.')goto 125 ! MMW line if(line(50:50).eq.'.')goto 125 ! IR line c 126 write(2,'(a)')line(1:len_trim(line)) goto 35 c c...deal with a transition line c 125 read(line(6:),'(6i3)',err=126)(iqn(i),i=1,nqn) indup=iqn(1) do 102 j=2,nqn indup=indup+iqn(j)*1000**(j-1) 102 continue c nst=6+nqn*3+2 read(line(nst:),'(6i3)',err=126)(iqn(i),i=1,nqn) indlw=iqn(1) do 103 j=2,nqn indlw=indlw+iqn(j)*1000**(j-1) 103 continue c c...loop over all parameters for which contributions are desired c do 76 j=1,ncontr ipard=numpar(j) c c c...find the upper level for the transitions among the energy levels from .EGY c do 37 i=1,nlevels if(indup.eq.indlev(i))then derup=deriv(i,ipard) ! upper level goto 38 endif 37 continue write(2,'(a)')line(1:len_trim(line)) ! unmodified output for missing derivative nunident=nunident+1 write(8,'(a)')line(1:len_trim(line))// ! output for missing derivative * ' <--- upper level not found in .EGY' goto 35 c c...find the lower level for the transitions among the energy levels from .EGY c 38 do 39 i=1,nlevels if(indlw.eq.indlev(i))then derlw=deriv(i,ipard) ! lower level goto 40 endif 39 continue c write(2,'(a)')line(1:len_trim(line)) ! unmodified output for missing derivative nunident=nunident+1 write(8,'(a)')line(1:len_trim(line))// ! output for missing derivative * ' <-- lower level not found in .EGY' goto 35 c c...work out parameter contribution c 40 contrib(j)=(derup-derlw)*param(ipard) deru(j)=derup derl(j)=derlw 76 continue ! end of contribution loop c c...output the .RES line modified by appended parameter contributions c (first convert the contribution to cm-1 if IR line) c sumcon=0.d0 do 88 j=1,ncontr sumcon=sumcon+contrib(j) 88 continue c colannot=0 ! extract end of line annotation if(len_trim(line).ge.86)then annot=line(86:len_trim(line)) colannot=74+(ncontr+1)*12 endif c if(line(50:50).eq.'.')then ! IR line do 86 j=1,ncontr contrib(j)=contrib(j)/29979.2458d0 86 continue sumcon=sumcon/29979.2458d0 write(line(74:),'(11f12.4)')(contrib(j),j=1,ncontr),sumcon else write(line(74:),'(11f12.3)')(contrib(j),j=1,ncontr),sumcon ! MW line endif c nident=nident+1 if(colannot.ne.0)then line(colannot:)=' '//annot(1:len_trim(annot)) endif c write(2,'(a)')line(1:len_trim(line)) ! main OUTPUT c if(40*(nident/40).eq.nident)then ! intermediate header write(2,90) * (pardesc(numpar(i)),i=1,ncontr),' SUM' write(2,'(1x)') endif 90 format(1x/73(1h-),11(2x,a10)) c write(7,153)line(1:5),' upper',indup,(deru(j),j=1,ncontr) ! identified levels write(7,153)' ', ' lower',indlw,(derl(j),j=1,ncontr) 153 format(1x,a5,a,i20,2x,1PE14.5,10e14.5) c goto 35 c c- - - end of transition reading loop- - - - - - - - - - - - - - - - - - - - - - c 36 close(1) close(2) close(7) close(8) c write(*,151) * ' Main output written to: ', * fileout(1:len_trim(fileout)) write(*,151)' Echo files also written:' write(*,'(1x)') write(*,152)nident, * ' identified transitions (elevels+derivs): ', * 'picontr_identified.out' write(*,152)nunident, * ' unidentified transitions: ', * 'picontr_unidentified.out' 151 format(1x/a,a) 152 format(1x,i7,a,a) c stop c C_____________________________________________________________________________ c c...Error conditions c c 1027 write(*,'(1x//'' Missing or unreadable PICONTR.INP file''//1x)') stop 1028 write(*,'(1x//'' Missing or unreadable .EGY file: '',a//1x)') * fileegy(1:len_trim(fileegy)) stop 1029 write(*,'(1x//'' Missing or unreadable .RES file: '',a//1x)') * fileres(1:len_trim(fileres)) stop 1030 write(*,'(1x//'' Missing or unreadable .PAR file: '',a//1x)') * filepar(1:len_trim(filepar)) stop 1142 write(*,'(1x// * '' Problem reading the third line of the .PAR file:''// * 1x,a//1x)') * line(1:len_trim(line)) c end c C_____________________________________________________________________________ C_____________________________________________________________________________