C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c KRA - substitution coordinates from single isotopic substitution in c a linear/symmetric/asymmetric top from Kraitchman's equations c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c - The program uses single isotopic equations straight out of c Gordy&Cook c - Experimental errors in rotational constants are propagated into c errors in coordinates c - The much more realistic errors according to the Costain criterion c are then added to the errors above, c - For planar molecules suitable planarity conditions can be built in c to use any two out of the three measured constants c c C Version 4a.IV.2017 ----- Zbigniew KISIEL ----- C 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 ca 95: Creation c 23.10.97: Addition of reference data and modifications to output c 30.05.99: Debugging of input c 1.04.01: Merging with symmetric top equations from another version c 14.10.02: Some output mods c 27.12.11: More explicit output C 28.08.12: Change to columnwise A,B,C input more convenient for copy/paste C 21.10.12: Automatic generation of .EVA type lines and general code cleanup C 27.11.13: More comments and symmetric case output bug eliminated C 4.04.17: Planar calculation output also echoed to the .EVA file C C_____________________________________________________________________________ C C Input: C C MOLNAM.KRA - The input file with name specified in response to the prompt C from the program. The use of the extension .KRA is recommended. C C Output: C C MOLNAM_KRA.OUT - main output file C C MOLNAM_KRA.EVA - stub file for the EVAL program containing the coordinates C and their uncertainties evaluated by KRA. C It is still necessary to assign signs, add two preceding C lines containing the descriptive comment and the number C of atoms, and append lines with declarations of internals C that are to be evaluated C_____________________________________________________________________________ c c The syntax of the input file: c c line0: comment line c line1: i j c line1a: k l c line2: A, dA \ c line2a: B, dB |- for the parent species c line2b: C, dC | c line3: Mass / c line4: A, dA \ c line4a: B, dB |- for a substituted species c line4b: C, dC | c line5: dMass / c . c . c lastline: end c c c - Line1 selects the type of top: c c i=0: defines that the parent is a diatomic, linear or a symmetric top c (off-axis asymmetric substitution is allowed for a symmetric top) c If i=0 then the value of j is irrelevant but is still necessary on c input so it is best to also set it to zero. c c i>0: defines that the parent is an asymmetric top, in which case i and j c i<0 can take values out of 1,2,3 (i.e. a,b,c). i cannot be equal to j, c and (i,j) selects the inertial plane for planar coordinate c calculation. c One of i and j can be negative in which case the respective moment c of inertia is calculated using the planar relation defined in line1a c c - Line1a is only necessary if one of i,j is negative and defines the c planar relation for calculation of that moment. For the most common case c of a planar molecule in the ab inertial plane there are three alternatives: c c line1 i,j: 1 2 -1 2 1 -2 c line1a k,l: none 3 -2 3 -1 c c and, for example, k,l = 3 -2 define I.a = I.c - I.b c For a non-zero inertia defect each these three cases gives a slightly c different result in the 'PLANAR' section, while the results in c the 'NONPLANAR' section will be identical. c c - Line 2: rotational constant A and its error dA for the parent species (MHz). c If the species is not an asymmetric top then set A to a value c large enough to pretend that it is infinity, say to 9999999.0 c c - Line 2a: B, dB for the parent (MHz) c - Line 2b: C, dC for the parent (MHz) c c - Line 3: mass of the parent isotopic species in u c c - Lines 4: rotational constant A and its error dA for the substituted c species (MHz) c c - Line 4a: B, dB for the substituted species (MHz) c - Line 4b: C, dC for the substituted species (MHz) c c - Line 5: the mass difference between the mass of substituted species and c that of the parent species (u), see the table below for some c common values. c c - Lines 4-5 can be repeated as many times as necessary, and are best c terminated with a line beginning with the word 'end' c (without the inverted commas) c c c c NOTE: c c - Lines that begin with a recognised comment flag i.e. c with either one of c,C,#,!,% are ignored so that c comment lines can be inserted anywhere in the data file. c The first line is the only exception as it is always treated as c comment describing the whole data set. c c - The X,Y,Z coordinates in the output are consistent with the c x,y,z indices in the equations in Gordy&Cook c c - The program tries to extract a six character descriptor from the comment lines C preceding each isotopic species. This descriptor will begin with the first c non blank character after the character defining the comment line. c The descriptor will be carried over to the .EVA file. C c - Remember that all substitution coordinates determined with Kraitchman's c equations arise from a square root operation so that their sign is c indeterminate. For convenience KRA prints all coordinates positive, and c not with +- sign in front of each. c It is, of course, necessary to assign signs to coordinates for their c use in determination of internal structural parameters, as for example c with EVAL. c C_____________________________________________________________________________ C C References: C C [1] W.Gordy and R.L.Cook, "Microwave Molecular Spectra", Wiley, C New York, 1984. c The relevant substitution relations are: c linear: Eq.13.46 c off axis symm.: Eq.13.56-13.57 c planar asymm: Eq.13.64-13.65 c general asymm. Eq.13.71-13.73 c Costain rule Eq.13.183 c C [2] C.C.Costain, "Further comments on the accuracy of r_s substitution structures", C Trans.Am.Crystallogr.Assoc. 2,157-164(1966). C C NOTE: what became the 'Costain rule' started off in the Costain paper C as an estimate of 0.0012/r based on NNO data (p.159). C The generalisation to the form 0.0015/a as found in Gordy&Cook C was made later. A preceding version of the G&C text can be found, C for example, in: C C [3] M.D.Harmony, V.W.Laurie, R.L.Kuczkowski, R.H.Schwendeman, C D.A.Ramsay, F.J.Lovas, W.J.Lafferty, A.G.Maki, C "Molecular Structures of Gas-Phase Polyatomic Molecules C Determined by Spectroscopic Methods", J. Phys. Chemn. Ref. Data, C 8, 619-721 (1979), p.624. C PDF version (bitmap scan) of that paper is freely available from C the NIST Standard Reference Data under the direct link: C http://www.nist.gov/data/PDFfiles/jpcrd146.pdf C_____________________________________________________________________________ C C Mass differences for use in KRA: C C Updated from "Atomic Weights and Isotopic Compositions" at: C http://www.nist.gov/pml/data/comp.cfm based on primary sources up to 2007 C C M.parent (u) M.isot (u) Delta.M (u) C c 1.H = 1.00782503207 2.H = 2.0141017778 1.0062767 c 3.H = 3.0160492777 2.0082242 C 12.C = 12. 13.C = 13.0033548378 1.0033548 c c 14.N = 14.0030740048 15.N = 15.0001088982 0.9970349 c c 16.O = 15.99491461956 17.O = 16.99913170 1.0042171 c 18.O = 17.9991610 2.0042464 c c 20.Ne = 19.9924401754 21.Ne = 20.99384668 1.0014065 c 22.Ne = 21.991385114 1.9989449 c c 28.Si = 27.9769265325 29.Si = 28.976494700 0.9995682 c 30.Si = 29.97377017 1.9968436 c c 32.S = 31.97207100 33.S = 32.97145876 0.9993878 c 34.S = 33.96786690 1.9957959 c 36.S = 35.96708076 3.9950098 c c 35.Cl = 34.96885268 37.Cl = 36.96590259 1.9970499 c c 40.Ar = 39.9623831225 36.Ar = 35.967545106 -3.9948380 c 38.Ar = 37.9627324 -1.9996507 c c 79.Br = 78.9183371 81.Br = 80.9162906 1.9979535 c C_____________________________________________________________________________ C c IMPLICIT REAL*8(A-H,O-Z) PARAMETER (CONV=505379.0089D0, maxat=100) ! CODATA2014 c integer iflagx,iflagy,ifirst,isecnd real*8 isot character planout(maxat)*100,plandesc*100 character*2 ixx,iyy,izz,ixxp,iyyp,irecal character filnam*50,coment*78,fsign*1,ssign*1,errmes*50, * isodes*50,filout*50,fileva*50,filgen*50 CHARACTER*3 AX(3) c COMMON RIX,RIY,RIZ, RIX1,RIY1,RIZ1, RM, * erix,eriy,eriz, erix1,eriy1,eriz1 COMMON /general/dpx,dpy,dpz,DIX,DIY,DIZ,xy,xz,yz common /moment/parent(3,2),isot(3,2),iflagx,iflagy,ifirst,isecnd common /isocom/isolen,isodes c c equivalence * (parent(1,1),ria),(parent(2,1),rib),(parent(3,1),ric), * (parent(1,2),eria),(parent(2,2),erib),(parent(3,2),eric), * (isot(1,1),ria1),(isot(2,1),rib1),(isot(3,1),ric1), * (isot(1,2),eria1),(isot(2,2),erib1),(isot(3,2),eric1) DATA AX/'a','b','c'/ isym=0 iflagx=2 iflagy=3 C 2 CONTINUE write(*,'(1x//)') WRITE(*,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | KRA - SINGLE ISOTOPIC SUBSTITUTION - Various ', * 'permutations',T79,'|'/ * ' | of Kraitchman''s equat', * 'ions for symmetric/asymmetric tops',T79,'|'/ * ' |',76(1H_),'|'/' version 4a.IV.2017',T64,'Zbigniew KISIEL'/) C c...establish and open input (unit 4) and output (units 3 and 7) files c 120 write(*,'(1x/'' Name of data file: '',$)') read(*,'(a)',err=120)filnam open(4,file=filnam,status='old',err=120) read(4,'(a)',err=120)coment c filgen=filnam do 314 nc=1,len_trim(filnam) if(filgen(nc:nc).eq.' '.or.filgen(nc:nc).eq.'.')goto 315 314 continue 315 nc=nc-1 filout=filgen(1:nc)//'_kra.out' fileva=filgen(1:nc)//'_kra.eva' c open(3,file=filout,status='unknown') write(3,3344) write(3,'(1x,78(1h-)/1x,a/1x,78(1h-)/)')coment write(*,'(/1x,78(1h-)/1x,a/1x,78(1h-)/)')coment c OPEN(7,FILE=fileva,STATUS='UNKNOWN') write(7,'(2a)')'! ',coment(1:len_trim(coment)) write(7,'(''!''/)') close(7) c c...read and process the control lines c 1121 errmes='Problems reading i,j (internally IFLAGX,IFLAGY)' 2002 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 2002 c read(4,*,err=2000,end=99)iflagx,iflagy c c...diatomic, linear or symmetric top c if(iflagx.eq.0)then isym=1 goto 121 endif c c...error trapping in control indices c write(errmes,'(2i7)')iflagx,iflagy errmes='Illegal i,j = '//errmes(1:14) if(iflagx.lt.0.and.iflagy.lt.0)goto 2000 if(iabs(iflagx).lt.1.or.iabs(iflagx).gt.3)goto 2000 if(iabs(iflagy).lt.1.or.iabs(iflagy).gt.3)goto 2000 if(iflagx.eq.iflagy)goto 2000 c c...choice of recalculated I for an asymmetric planar molecule c if(iflagx.lt.0.or.iflagy.lt.0)then errmes='Problems reading k,l (internally IFIRST,ISECND)' 2003 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 2003 c read(4,*,err=2000,end=99)ifirst,isecnd write(errmes,'(2i7)')ifirst,isecnd errmes='Illegal k,l = '//errmes(1:14) if(iabs(ifirst).lt.1.or.iabs(ifirst).gt.3)goto 2000 if(iabs(isecnd).lt.1.or.iabs(isecnd).gt.3)goto 2000 if(iabs(ifirst).eq.iabs(isecnd))goto 2000 c fsign=' ' ssign='+' if(ifirst.lt.0)fsign='-' if(isecnd.lt.0)ssign='-' c write(3,'(1x)') write(*,'(1x)') c if(iflagx.lt.0)then write(3,37)ax(iabs(iflagx)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)), * ax(iabs(iflagy)) write(*,37)ax(iabs(iflagx)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)), * ax(iabs(iflagy)) write(plandesc,37)ax(iabs(iflagx)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)), * ax(iabs(iflagy)) endif c if(iflagy.lt.0)then write(3,38)ax(iabs(iflagx)), * ax(iabs(iflagy)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)) write(*,38)ax(iabs(iflagx)), * ax(iabs(iflagy)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)) write(plandesc,38)ax(iabs(iflagx)), * ax(iabs(iflagy)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)) endif else write(*,36)ax(iabs(iflagx)),ax(iabs(iflagy)) write(3,36)ax(iabs(iflagx)),ax(iabs(iflagy)) write(plandesc,36) * ax(iabs(iflagx)),ax(iabs(iflagy)) endif c write(3,'(1x)') write(*,'(1x)') c 36 format(' Planar calculation will be made from I.',a1, * ' and I.',a1) 37 format(' Planar calculation will be made from I.',a1, * ' = ',a1,'I.',a1,' ',a1,' I.',a1,' and I.',a1) 38 format(' Planar calculation will be made from I.',a1, * ' and I.',a1,' = ',a1,'I.',a1,' ',a1,' I.',a1) c c...data for the parent species c errmes='Problem reading one of (A,dA), (B,dB), (C,dC)' 121 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 121 read(4,*,err=2000,end=99)a,era read(4,*,err=2000,end=99)b,erb read(4,*,err=2000,end=99)c,erc c errmes='Problem reading mass of parent' 123 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 123 read(4,*,err=2000,end=99)tm c RIA=CONV/A RIB=CONV/B RIC=CONV/C eria=conv*era/a**2 erib=conv*erb/b**2 eric=conv*erc/c**2 defect=RIC-RIA-RIB rkappa=(2.d0*b-a-c)/(a-c) APLUSB=(A+B)/C PA=0.5D0*(RIB+RIC-RIA) PB=0.5D0*(RIA+RIC-RIB) PC=0.5D0*(RIA+RIB-RIC) write(3,'(1x,78(1h-))') C C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C c...data for the isotopic species (looped over from label 124 as many times c as available in the input file) c errmes='Problem reading one of isotopic (A,dA), (B,dB), (C,dC)' isolen=0 natom=0 c 124 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 124 if(ierr.eq.2.or.ierr.eq.3)then write(*,'(/'' end of data''/)') goto 99 endif c READ(4,*,err=2000,end=99)A1,era1 if(a1.le.0.0)goto 2000 READ(4,*,err=2000,end=99)B1,erb1 READ(4,*,err=2000,end=99)C1,erc1 c errmes='Problem reading isotopic mass change' 126 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 126 READ(4,*,err=2000,end=99)DM if(dm.eq.0.0.or.dabs(dm).gt.100.)goto 2000 RM=TM*DM/(TM+DM) C RIA1=CONV/A1 RIB1=CONV/B1 RIC1=CONV/C1 eria1=conv*era1/a1**2 erib1=conv*erb1/b1**2 eric1=conv*erc1/c1**2 defec1=RIC1-RIA1-RIB1 rkapp1=(2.d0*b1-a1-c1)/(a1-c1) APLUS1=(A1+B1)/C1 PA1=0.5D0*(RIB1+RIC1-RIA1) PB1=0.5D0*(RIA1+RIC1-RIB1) PC1=0.5D0*(RIA1+RIB1-RIC1) C c c...Linear+symmetric top KRAITCHMAN calculation c if(isym.eq.1)then if(rib.eq.ric)then xx=c yy=b zz=a erxx=erc eryy=erb erzz=era RIX=RIc RIY=RIb RIZ=RIa erix=eric eriy=erib eriz=eria xx1=c1 yy1=b1 zz1=a1 erxx1=erc1 eryy1=erb1 erzz1=era1 RIX1=RIc1 RIY1=RIb1 RIZ1=RIa1 erix1=eric1 eriy1=erib1 eriz1=eria1 iflagz=1 iflagy=2 else xx=a yy=b zz=c erxx=era eryy=erb erzz=erc RIX=RIa RIY=RIb RIZ=RIc erix=eria eriy=erib eriz=eric xx1=a1 yy1=b1 zz1=c1 erxx1=era1 eryy1=erb1 erzz1=erc1 RIX1=RIa1 RIY1=RIb1 RIZ1=RIc1 erix1=eria1 eriy1=erib1 eriz1=eric1 iflagz=3 iflagy=2 endif c CALL SUBSTS(ZS,YS,erzs,erys,z2,y2) <--- SUBSTS c if(zs.ne.0.d0)then ERCOSZ=DSQRT((0.0015D0/ZS)**2+ERZS**2) else ercosz=0.d0 endif if(ys.ne.0.d0)then ERCOSY=DSQRT((0.0015D0/YS)**2+ERYS**2) else ercosy=0.d0 endif iyy=' ' izz=' ' if(y2.lt.0.d0)iyy='*i' if(z2.lt.0.d0)izz='*i' c write(3,'(1x,78(1h-))') WRITE(3,10) 10 FORMAT(1X/' The parent species:'/) WRITE(3,33)' X, Y, Z = ',XX,YY,ZZ write(3,33)' eX, eY, eZ = ',erxx,eryy,erzz WRITE(3,33)' IX, IY, IZ = ',RIx,RIy,RIz WRITE(3,33)'eIX,eIY,eIZ = ',erIx,eRIy,eRIz WRITE(3,'(1x)') WRITE(3,35)' Mass = ',TM WRITE(3,15) 15 FORMAT(1X/' The isotopic species:'/) WRITE(3,33)' X, Y, Z = ',xx1,yy1,zz1 write(3,33)' eX, eY, eZ = ',erxx1,eryy1,erzz1 WRITE(3,33)' IX, IY, IZ = ',RIx1,RIy1,RIz1 WRITE(3,33)'eIX,eIY,eIZ = ',eRIx1,eRIy1,eRIz1 33 FORMAT(1X,A,3F20.8) WRITE(3,'(1x)') WRITE(3,35)' Mass change = ',DM WRITE(3,35)' Total mass = ',TM+DM WRITE(3,35)' M DM/(M+DM) = ',RM 35 format(1x,a,f20.8) WRITE(3,43) 43 FORMAT(/' KRAITCHMAN RESULTS:'/) c WRITE(3,'(28x,a1,22x,a1/)')AX(iabs(IFLAgz)),ax(iabs(iflagy)) WRITE(3,334)' SYMMETRIC: ',zs,izz ,erzs, * Ys,iyy ,erys WRITE(3,334)'+Costain err.',zs,izz ,ercosz, * Ys,iyy ,ercosy c write(3,'(1x/1x,78(1h-))') c c...lines for EVA c write(3,'(1x/1x,78(1h-))') C OPEN(7,FILE=fileva,STATUS='UNKNOWN',ACCESS='APPEND') write(7,44)isodes(1:6),zs,izz,ercosz, * Ys,iyy,ercosy, * 0.0d0,' ',0.0d0 close(7) c44 format(a6,f12.5,a2,f8.5,2x, f10.5,a2,f8.5,2x, f10.5,a2,f8.5) goto 124 endif C C...Planar KRAITCHMAN calculation C CALL SUBSTP(XP,YP,erxp,eryp,x2,y2) <--- SUBSTP ERCOSX=DSQRT((0.0015D0/XP)**2+ERXP**2) ERCOSY=DSQRT((0.0015D0/YP)**2+ERYP**2) ixxp=' ' iyyp=' ' if(x2.lt.0.d0)ixxp='*i' if(y2.lt.0.d0)iyyp='*i' c if(iflagx.lt.0.or.iflagy.lt.0)then if(iflagx.lt.0)irecal='I'//ax(iabs(iflagx)) if(iflagy.lt.0)irecal='I'//ax(iabs(iflagy)) srix=rix sriy=riy srix1=rix1 sriy1=riy1 endif C c...General KRAITCHMAN calculation c xx=a yy=b zz=c erxx=era eryy=erb erzz=erc RIX=RIA RIY=RIB RIZ=RIC erix=eria eriy=erib eriz=eric xx1=a1 yy1=b1 zz1=c1 erxx1=era1 eryy1=erb1 erzz1=erc1 RIX1=RIA1 RIY1=RIB1 RIZ1=RIC1 erix1=eria1 eriy1=erib1 eriz1=eric1 CALL SUBST(X,Y,Z,erx,ery,erz,x2,y2,z2,R) <--- SUBST ERCX=DSQRT((0.0015D0/X)**2+ERX**2) ERCY=DSQRT((0.0015D0/Y)**2+ERY**2) ERCZ=DSQRT((0.0015D0/Z)**2+ERZ**2) ERR=dsqrt( (X*ercx)**2+(Y*ercy)**2+(Z*ercz)**2 )/R c ixx=' ' iyy=' ' izz=' ' if(x2.lt.0.d0)ixx='*i' if(y2.lt.0.d0)iyy='*i' if(z2.lt.0.d0)izz='*i' c write(3,'(1x,78(1h-))') WRITE(3,10) WRITE(3,33)' X, Y, Z = ',XX,YY,ZZ write(3,33)' eX, eY, eZ = ',erxx,eryy,erzz WRITE(3,33)' IX, IY, IZ = ',RIx,RIy,RIz WRITE(3,33)'eIX,eIY,eIZ = ',erIx,eRIy,eRIz WRITE(3,33)' PX, PY, PZ = ',PA,PB,PC c WRITE(3,33)'In.Defect= ',defect c WRITE(3,33)'Kappa = ',rkappa c WRITE(3,33)'A+B /C = ',APLUSB WRITE(3,'(1x)') WRITE(3,35)' Mass = ',TM WRITE(3,15) WRITE(3,33)' X, Y, Z = ',xx1,yy1,zz1 write(3,33)' eX, eY, eZ = ',erxx1,eryy1,erzz1 WRITE(3,33)' IX, IY, IZ = ',RIx1,RIy1,RIz1 WRITE(3,33)'eIX,eIY,eIZ = ',eRIx1,eRIy1,eRIz1 WRITE(3,33)' PX, PY, PZ = ',PA1,PB1,PC1 c WRITE(3,33)'In.Defect= ',defec1 c WRITE(3,33)'Kappa = ',rkapp1 c WRITE(3,33)'A+B /C = ',APLUS1 WRITE(3,'(1x)') WRITE(3,35)' Mass change = ',DM WRITE(3,35)' Total mass = ',TM+DM WRITE(3,35)' M DM/(M+DM) = ',RM WRITE(3,43) c if(iflagx.lt.0.or.iflagy.lt.0) * write(3,336)srix, sriy, irecal, * srix1,sriy1,irecal 336 FORMAT(' Modified values used for PLANAR calculation:'/ * ' parent IX,IY =',2f20.8,' <-- ',a,' recalculated'/ * ' isotopicIX,IY =',2f20.8,' <-- ',a,' recalculated'/) c WRITE(3,'(28x,a1,22x,a1/)')AX(iabs(IFLAgx)),ax(iabs(iflagy)) WRITE(3,334)' PLANAR: ',XP,ixxp,erxp, * YP,iyyp,eryp WRITE(3,334)'+Costain err.',XP,ixxp,ercosx, * YP,iyyp,ercosy 334 FORMAT(1X,A,3X,F12.5,a2,'+-',F8.5, * F11.5,a2,'+-',F8.5) c WRITE(3,'(1x/28x,a1,22x,a1,22x,a1/)')'a','b','c' WRITE(3,34)' NONPLANAR: ',X,ixx,erx,Y,iyy,ery, * Z,izz,erz WRITE(3,34)'+Costain err.',X,ixx,ercx,Y,iyy,ercy, * Z,izz,ercz 34 FORMAT(1X,A,3X,F12.5,a2,'+-',F8.5, * F11.5,a2,'+-',F8.5, * F11.5,a2,'+-',F8.5) c write(3,134)r,err 134 format(1x/64x,'R=',f9.5,' +-',F8.5) c WRITE(3,'(1x)') WRITE(3,33)'DIX,DIY,DIZ = ',DIX,DIY,DIZ WRITE(3,33)'DPX,DPY,DPZ = ',DPX,DPY,DPZ WRITE(3,33)'IXY,IXZ,IYZ = ',XY ,XZ ,YZ c c...lines for EVA c write(3,'(1x/1x,78(1h-))') C OPEN(7,FILE=fileva,STATUS='UNKNOWN',ACCESS='APPEND') write(7,44)isodes(1:6),X,ixx,ercx, * Y,iyy,ercy, * Z,izz,ercz 44 format(a6,f12.5,a2,f8.5,2x, f10.5,a2,f8.5,2x, f10.5,a2,f8.5) close(7) c natom=natom+1 izz=' ' write(planout(natom),44) * isodes(1:6),XP,ixxp,ercosx, * YP,iyyp,ercosy, * 0.0,izz,0.0 C isolen=0 GOTO 124 c C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c 99 close(3) close(4) c if(iflagx.ne.0)then OPEN(7,FILE=fileva,STATUS='UNKNOWN',ACCESS='APPEND') write(7,'(1x//''!''/''!'',2a/''!''/)')plandesc(1:20), * plandesc(29:len_trim(plandesc)) do 210 n=1,natom write(7,'(a)')planout(n)(1:len_trim(planout(n))) 210 continue close(7) endif c write(*,'(1x/'' Output has been written to:''/5x,a/5x,a)') * filout(1:len_trim(filout)),fileva(1:len_trim(fileva)) c write(*,199) 199 format(1x// * ' NOTE: The coordinate values in the .EVA file still need to be'/ * ' assigned suitable signs.'/ * ' These lines also have to be preceded by a line with'/ * ' a single number specifying the number of coordinate'/ * ' lines and followed by internal coordinate definitions.' * //) c STOP c 2000 write(*,2001)errmes 2001 Format(1x//'**** INPUT ERROR: ',a//) c stop END c c------------------------------------------------------------------------ c SUBROUTINE SUBST(X,Y,Z,erx,ery,erz,x2,y2,z2,R) c C...ASYMMETRIC TOP - general case c [Gordy&Cook, eq. 13.71-13.74 p.662] c IMPLICIT REAL*8(A-Z) COMMON IX,IY,IZ,IX1,IY1,IZ1,RM, * erix,eriy,eriz,erix1,eriy1,eriz1 COMMON /general/dpx,dpy,dpz,DIX,DIY,DIZ,xy,xz,yz c DIX=IX1-IX DIY=IY1-IY DIZ=IZ1-IZ XY=IX-IY XZ=IX-IZ YZ=IY-IZ DPX=0.5D0*(-DIX+DIY+DIZ) DPY=0.5D0*(-DIY+DIZ+DIX) DPZ=0.5D0*(-DIZ+DIX+DIY) c X2=DPX*(1+DPY/XY)*(1+DPZ/XZ)/RM X=DSQRT(DABS(X2)) c y2=DPY*(1+DPZ/YZ)*(1-DPX/XY)/RM Y=DSQRT(DABS(y2)) c z2=DPZ*(1-DPX/XZ)*(1-DPY/YZ)/RM Z=DSQRT(DABS(z2)) c A=(DIX+DIY+DIZ)/(2.D0*RM) IF(A.LT.0.D0)WRITE(*,'(1x,a)')' WARNING (negative value) ' R=DSQRT(DABS(A)) c c c...Errors in Kraitchman quantities Delta(P.x), Delta(P.y) and Delta(P.z) c erdpx=0.5d0*dsqrt(erix**2+eriy**2+eriz**2+erix1**2+eriy1**2+ $ eriz1**2) erdpy=erdpx erdpz=erdpx c c 2 2 2 c....Errors in squares of Kraitchman coordinates ie. in X , Y and Z c c The following expressions for the error in the square of the Kraitchman c coordinate have been evaluated with DERIVE by starting from the two c expressions commented out for each error c c X2:=dpx*(dpz+ix-iz)*(dpy+ix-iy)/(rm*(ix-iy)*(ix-iz)) c c DIF(X2,dpx)**2*erdpx**2+DIF(X2,dpy)**2*erdpy**2+DIF(X2,dpz)**2*er c $dpz**2+DIF(X2,ix)**2*erix**2+DIF(X2,iy)**2*eriy**2+DIF(X2,iz)**2*e c $riz**2 c erx2= $dpx**2*(dpy**2*(dpz**2*(erix**2*(2*ix-iy-iz)**2+eriy**2*(ix-iz)**2 $+eriz**2*(ix-iy)**2)+2*dpz*(ix-iz)**2*(erix**2*(2*ix-iy-iz)+eriy** $2*(ix-iz))+(ix-iz)**2*(erdpz**2*(ix-iy)**2+(ix-iz)**2*(erix**2+eri $y**2)))+2*dpy*(ix-iy)**2*(dpz**2*(erix**2*(2*ix-iy-iz)+eriz**2*(ix $-iy))+dpz*erix**2*(ix-iz)**2+erdpz**2*(ix-iz)**2*(ix-iy))+(ix-iy)* $*2*(dpz**2*(erdpy**2*(ix-iz)**2+(ix-iy)**2*(erix**2+eriz**2))+2*dp $z*erdpy**2*(ix-iz)**3+(ix-iz)**2*(erdpy**2*(ix**2-2*ix*iz+iz**2)+e $rdpz**2*(ix**2-2*ix*iy+iy**2))))/(rm**2*(ix-iy)**4*(ix-iz)**4)+erd $px**2*(dpz+ix-iz)**2*(dpy+ix-iy)**2/(rm**2*(ix-iy)**2*(ix-iz)**2) erx2=dsqrt(erx2) c c Y2:=dpy/rm*(1+dpz/(iy-iz))*(1+dpx/(iy-ix)) c c DIF(Y2,dpx)**2*erdpx**2+DIF(Y2,dpy)**2*erdpy**2+DIF(Y2,dpz)**2*er c $dpz**2+DIF(Y2,ix)**2*erix**2+DIF(Y2,iy)**2*eriy**2+DIF(Y2,iz)**2*e c $riz**2 c ery2= $(dpx**2*(dpy**2*(dpz**2*(erix**2*(iy-iz)**2+eriy**2*(ix-2*iy+iz)** $2+eriz**2*(ix-iy)**2)+2*dpz*(iy-iz)**2*(erix**2*(iy-iz)+eriy**2*(- $ix+2*iy-iz))+(iy-iz)**2*(erdpz**2*(ix-iy)**2+(iy-iz)**2*(erix**2+e $riy**2)))+erdpy**2*(iy-iz)**2*(ix-iy)**2*(dpz+iy-iz)**2)-2*dpx*(ix $-iy)**2*(dpy**2*(dpz**2*(eriy**2*(ix-2*iy+iz)+eriz**2*(ix-iy))-dpz $*eriy**2*(iy-iz)**2+erdpz**2*(iy-iz)**2*(ix-iy))+erdpy**2*(iy-iz)* $*2*(ix-iy)*(dpz+iy-iz)**2)+(ix-iy)**2*(dpy**2*(dpz**2*(erdpx**2*( $iy-iz)**2+(ix-iy)**2*(eriy**2+eriz**2))+2*dpz*erdpx**2*(iy-iz)**3 $+(iy-iz)**2*(erdpz**2*(ix-iy)**2+erdpx**2*(iy**2-2*iy*iz+iz**2))) $+erdpy**2*(iy-iz)**2*(ix-iy)**2*(dpz+iy-iz)**2))/(rm**2*(ix-iy)**4 $*(iy-iz)**4) ery2=dsqrt(ery2) c c Z2:=dpz/rm*(1+dpx/(iz-ix))*(1+dpy/(iz-iy)) c c DIF(Z2,dpx)**2*erdpx**2+DIF(Z2,dpy)**2*erdpy**2+DIF(Z2,dpz)**2*er c $dpz**2+DIF(Z2,ix)**2*erix**2+DIF(Z2,iy)**2*eriy**2+DIF(Z2,iz)**2*e c $riz**2 c erz2= $(dpx**2*(dpy**2*(dpz**2*(erix**2*(iy-iz)**2+eriy**2*(ix-iz)**2+eri $z**2*(ix+iy-2*iz)**2)+erdpz**2*(iy-iz)**2*(ix-iz)**2)-2*dpy*(iy-iz $)**2*(dpz**2*(erix**2*(iy-iz)+eriz**2*(ix+iy-2*iz))+erdpz**2*(iy-i $z)*(ix-iz)**2)+(iy-iz)**2*(dpz**2*(erdpy**2*(ix-iz)**2+(iy-iz)**2* $(erix**2+eriz**2))+erdpz**2*(iy-iz)**2*(ix-iz)**2))-2*dpx*(ix-iz)* $*2*(dpy**2*(dpz**2*(eriy**2*(ix-iz)+eriz**2*(ix+iy-2*iz))+erdpz**2 $*(iy-iz)**2*(ix-iz))-dpy*(iy-iz)**2*(dpz**2*eriz**2+2*erdpz**2*(iy $-iz)*(ix-iz))+(iy-iz)**2*(ix-iz)*(dpz**2*erdpy**2+erdpz**2*(iy-iz) $**2))+(ix-iz)**2*(dpy**2*(dpz**2*(erdpx**2*(iy-iz)**2+(ix-iz)**2* $(eriy**2+eriz**2))+erdpz**2*(iy-iz)**2*(ix-iz)**2)-2*dpy*(iy-iz)** $3*(dpz**2*erdpx**2+erdpz**2*(ix-iz)**2)+(iy-iz)**2*(dpz**2*(erdpy $**2*(ix-iz)**2+erdpx**2*(iy-iz)**2)+erdpz**2*(iy-iz)**2*(ix-iz)** $2)))/(rm**2*(ix-iz)**4*(iy-iz)**4) erz2=dsqrt(erz2) c c...Errors in Kraitchman coordinates X, Y and Z c erx=0.5d0*erx2/x ery=0.5d0*ery2/y erz=0.5d0*erz2/z c c RETURN END c c------------------------------------------------------------------------ c SUBROUTINE SUBSTS(Z,Y,erz,ery,z2,y2) c C...SYMMETRIC TOP [on axis: Gordy&Cook, eq. 13.46, p.658] c [off-axis: Gordy&Cook, eq. 13.56,13.57, p.660] c IMPLICIT REAL*8(A-Z) COMMON IX,IY,IZ,IX1,IY1,IZ1,RM, * erix,eriy,eriz,erix1,eriy1,eriz1 c c...substitution on the symmetry axis c if((ix1.eq.iy1).or.(iy1.eq.iz1))then DIy=Iy1-Iy z2=diy/rm z=dsqrt(abs(z2)) ERZ2=DSQRT(ERIY1**2+ERIY**2) erz=erz2/(2.d0*z) y=0 y2=0 ery=0 return endif c c...substitution off the symmetry axis c y2=((ix1-iy1)*(ix1-iz1)) / (rm * (ix1-iy1-iz1+iy)) z2=((iy1-iy )*(iz1-iy )) / (rm * (iy1+iz1-ix1-iy)) z=dsqrt(abs(z2)) y=dsqrt(abs(y2)) C ERY2=DSQRT( $erix1**2*(ix1**2+2*ix1*(iy-iy1-iz1)+iy*(-iy1-iz1)+iy1**2+iy1*iz1+i $z1**2)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy**2*(ix1-iz1)**2*(ix1-iy $1)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy1**2*(iy-iz1)**2*(ix1-iz1)** $2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriz1**2*(iy-iy1)**2*(ix1-iy1)**2/(r $m**2*(ix1+iy-iy1-iz1)**4) ) ery=erz2/(2.d0*y) ERZ2=DSQRT( $erix1**2*(iy-iz1)**2*(iy-iy1)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy* $*2*(ix1*(2*iy-iy1-iz1)+iy**2-2*iy*(iy1+iz1)+iy1**2+iy1*iz1+iz1**2) $**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy1**2*(iy-iz1)**2*(ix1-iz1)**2/ $(rm**2*(ix1+iy-iy1-iz1)**4)+eriz1**2*(iy-iy1)**2*(ix1-iy1)**2/(rm* $*2*(ix1+iy-iy1-iz1)**4) ) erz=erz2/(2.d0*z) c c RETURN END c c------------------------------------------------------------------------ c SUBROUTINE SUBSTP(X,Y,erx,ery,x2,y2) c C..ASYMMETRIC TOP P L A N A R C [Gordy&Cook, eq. 13.64,13.65 p.661] C C All input parameters are passed through the COMMON blocks c c iflagx, iflagy - inertial coordinates to be treated as x and y c 1,2,3 for a,b,c. A negative value can be specified for c one of iflagx,iflagy, in which case the necessary quantities c for this axis are evaluated from moments for the remaining c two axes using the planarity relations: c I.b+ I.a= I.c c DI.b+DI.a=DI.c c as defined for the case in question by ifirst,isecnd, ie. if c iflagx is negative then c I.iflagx = I.ifirst + I.isecnd c c X, Y - r.s coordinates from planar equations c ERX, ERY - errors on r.s coordinates from error propagation c X2, Y2 - squares of r.s coordinates as calculated with the equations, c these can be negative which can be used for marking imaginary c coordinates c IMPLICIT REAL*8(A-Z) integer iflagx,iflagy,ifirst,isecnd,i1,i2 COMMON IX,IY,IZ, IX1,IY1,IZ1, RM, * erix,eriy,eriz, erix1,eriy1,eriz1 common /moment/parent(3,2),isot(3,2),iflagx,iflagy,ifirst,isecnd c if(iflagx.gt.0)then ix=parent(iflagx,1) ix1=isot(iflagx,1) DIX=IX1-IX ERIX=PARENT(IFLAGX,2) ERDIX=DSQRT(PARENT(IFLAGX,2)**2+ISOT(IFLAGX,2)**2) else i1=iabs(ifirst) i2=iabs(isecnd) s1=isign(1,ifirst) s2=isign(1,isecnd) dix=s1*(isot(i1,1)-parent(i1,1)) * +s2*(isot(i2,1)-parent(i2,1)) c ix=s1*parent(i1,1)+s2*parent(i2,1) ix1=s1*isot(i1,1)+s2*isot(i2,1) ERDIX=DSQRT(ISOT(i1,2)**2+PARENT(i1,2)**2 * +ISOT(i2,2)**2+PARENT(i2,2)**2) ERIX=DSQRT(PARENT(i1,2)**2+PARENT(i2,2)**2) endif c if(iflagy.gt.0)then iy=parent(iflagy,1) iy1=isot(iflagy,1) DIY=IY1-IY ERIY=PARENT(IFLAGY,2) ERDIY=DSQRT(PARENT(IFLAGY,2)**2+ISOT(IFLAGY,2)**2) else i1=iabs(ifirst) i2=iabs(isecnd) s1=isign(1,ifirst) s2=isign(1,isecnd) diy=s1*(isot(i1,1)-parent(i1,1)) * +s2*(isot(i2,1)-parent(i2,1)) iy=s1*parent(i1,1)+s2*parent(i2,1) iy1=s1*isot(i1,1)+s2*isot(i2,1) ERDIy=DSQRT(ISOT(i1,2)**2+PARENT(i1,2)**2 * +ISOT(i2,2)**2+PARENT(i2,2)**2) ERIy=DSQRT(PARENT(i1,2)**2+PARENT(i2,2)**2) endif XY=IX-IY ERXY=DSQRT(ERIX**2+ERIY**2) c WRITE(*,33)'DIX,DIY,XY = ',DIX,DIY,XY c WRITE(3,33)'DIX,DIY,XY = ',DIX,DIY,XY c33 FORMAT(1X,A,3F20.5) C x2=DIY*(1.d0+DIX/XY)/RM X=DSQRT(DABS(x2)) y2=DIX*(1.d0-DIY/XY)/RM Y=DSQRT(DABS(y2)) c erx2= *DSQRT(dix**2*(diy**2*erxy**2+erdiy**2*xy**2)+2*dix*erdiy**2*xy**3 * +diy**2*erdix**2*xy**2+erdiy**2*xy**4)/(xy**2*ABS(rm)) erx=erx2/(2.d0*x) ery2= *DSQRT(dix**2*(diy**2*erxy**2+erdiy**2*xy**2)+diy**2*erdix**2*xy**2 * -2*diy*erdix**2*xy**3+erdix**2*xy**4)/(xy**2*ABS(rm)) ery=ery2/(2.d0*y) c RETURN END c c------------------------------------------------------------------------ c subroutine test(iodev,ierr) character line*80,isodes*50 common /isocom/isolen,isodes c c...read a line from the input file (on IODEV), test it and c return IERR as follows: c c IERR=0 the line is a bona fide data line, in which case the input c file is backspaced for numerical input c IERR=1 the line is a comment line, in which case it is echoed c directly to the output c IERR=2 end of data, as specified by the string 'end' has been encountered c IERR=3 end of file has been encountered c read(iodev,'(a)',err=3,end=3)line c if(line(1:3).eq.'end')then ierr=3 return endif c if(line(1:1).eq.'c'.or.line(1:1).eq.'C'.or.line(1:1).eq.'!' * .or.line(1:1).eq.'#'.or.line(1:1).eq.'%')then ierr=1 if(len_trim(line).gt.1)then write(3,'(1x,a)')line(2:) if(len_trim(line).gt.isolen)then isolen=len_trim(line) isodes=line(2:len_trim(line)) do 1 n=1,len_trim(line)-1 if(isodes(1:1).ne.' ')then write(*,'(1x,2a)')'processing: ', * isodes(1:len_trim(isodes)) return endif isodes=isodes(2:) 1 continue endif endif else backspace(iodev) ierr=0 endif return c 3 ierr=2 return end c c------------------------------------------------------------------------ c------------------------------------------------------------------------