c ****************************************************************** c SHLX2TEP - SHELX to ORTEP-III Data Conversion Program c c Martin Kroeker c c Version 1.0 - February 2004 fix for serious bug that converted all +1 c translation elements resulting from c combinations of (+0.5) (+0.5), e.g. c centering and mirror plane to +1/2 symbol c on symmetry card c Version 0.99 - September 2002 accept LATT cards with whitespace between c LATT and number (as written by PLATON) c Version 0.98 - August 2002 ignore REM among the atom cards c Version 0.97 - February 2002 (handle P-1 given as LATT 1 without SYMM) c Version 0.96 - August 2001 (handle symmetry cards in PLATON-generated c shelx files - these have lots of whitespace c automatically retry reading with old shelx c format if new format cannot be read) c Version 0.95 - April 2000 (now ignores more shelx commands ) c Version 0.94 - December 1999 based on c Version 0.9 BETA - August 27, 1996 (by Mike Burnett,ORNL) c c Send comments, questions, problems, etc. to c martin@ruby.chemie.uni-freiburg.de c c but *NOT* to ORNL staff c c ****************************************************************** common /CELL/ A, B, C, ALPHA, BETA, GAMMA common /Uij/ U11, U22, U33, U12, U13, U23 common /SYMM/ fs(3,3,96), ts(3,96), numsymm common /FVRS/ fvar(10),ueq character line*80, lab*4, char1 CHARACTER*2 ATSYM(40),SYMBLS(103) INTEGER NSFAC,ORDNL(40) logical stereo,oldshx dimension fsym(3,4) character*36 maksym DATA SYMBLS/'H ', 'HE', 1 'LI','BE','B ','C ','N ','O ','F ','NE', 2 'NA','MG','AL','SI','P ','S ','CL','AR', 3 'K ','CA', A 'SC','TI','V ','CR','MN','FE','CO','NI','CU','ZN', X 'GA','GE','AS','SE','BR','KR', 4 'RB','SR', A 'Y ','ZR','NB','MO','TC','RU','RH','PD','AG','CD', X 'IN','SN','SB','TE','I ','XE', 5 'CS','BA', B 'LA','CE','PR','ND','PM','SM','EU','GD','TB','DY', C 'HO','ER','TM','YB','LU', A 'HF','TA','W ','RE','OS','IR','PT','AU','HG', X 'TL','PB','BI','PO','AT','RN', 6 'FR','RA', B 'AC','TH','PA','U ','NP','PU','AM','CM','BK','CF', C 'ES','FM','MD','NO','LR'/ PI=3.1415926536 1 format(a) write(6,*)'SHELX2TEP 1.0 (c) 1996 Mike Burnett (ORNL)' write(6,*)' (c) 1999-2004 Martin Kroeker (Freiburg)' write(6,*) +'please report bugs to martin@ruby.chemie.uni-freiburg.de' write(6,*) c *** get SHELX file name write (6,100) 100 format(' Enter SHELX file name or "exit": ',$) read (5,1) line if (line(1:4).eq.'exit' .or. line(1:4).eq.'EXIT') stop open(3,file=line,status='old',err=105) go to 110 105 write (6,108) line 108 format(' Cannot open ',a) stop c *** open output file 110 open(2,file='TEP.IN',status='unknown') c *** shelx93 and shelx97 differ in coordinate format oldshx=.false. write(6,111) 111 format(' Assume shelx93 (or earlier) file format ? (Y or N) ',$) read (5,1) char1 if (char1.eq.'y' .or. char1.eq.'Y') oldshx=.true. c *** ask about stereo stereo=.false. write (6,112) 112 format(' Stereo drawing? (Y or N) ',$) read (5,1) char1 if (char1.eq.'y' .or. char1.eq.'Y') stereo=.true. c *** title 115 read (3,1,end=999) line if (line(1:4) .ne. 'TITL') go to 115 write (2,1) line(6:) c *** cell parameters 118 read (3,1,end=999) line if (line(1:4) .ne. 'CELL') go to 118 open(4,status='scratch') write (4,1) line(5:iend(line)) backspace(4) read (4,*) f,A,B,C,ALPHA,BETA,GAMMA close(4) write (2,120) A,B,C,ALPHA,BETA,GAMMA 120 format('1',f8.4,5f9.4) c *** lattice latt=1 122 read (3,1,end=125) line if (line(1:4) .ne. 'LATT') go to 122 backspace(3) read (3,124) lab,latt 124 format(a,1x,i20) go to 128 125 rewind(3) c *** symmetry 128 numsymm=0 line='SYMM X,Y,Z' call addsymm(line) 130 read (3,1,end=132) line if (line(1:4) .ne. 'SYMM') go to 130 call addsymm(line) go to 130 132 rewind(3) c centrosymmetric if (latt.gt.0) then do 135 i=1,numsymm inew=i+numsymm do 135 j=1,3 ts(j,inew)=ts(j,i) do 135 k=1,3 fs(k,j,inew)=-fs(k,j,i) 135 continue numsymm=2*numsymm end if c lattice centering translations ilatt=iabs(latt) if (ilatt.gt.1) then isymm=numsymm if (ilatt.eq.2) call addlatt(isymm,.5,.5,.5) if (ilatt.eq.3) then call addlatt(isymm,2./3.,1./3.,1./3.) call addlatt(isymm,1./3.,2./3.,2./3.) end if if (ilatt.eq.4) then call addlatt(isymm,0.,.5,.5) call addlatt(isymm,.5,0.,.5) call addlatt(isymm,.5,.5,0.) end if if (ilatt.eq.5) call addlatt(isymm,0.,.5,.5) if (ilatt.eq.6) call addlatt(isymm,.5,0.,.5) if (ilatt.eq.7) call addlatt(isymm,.5,.5,0.) end if c print symmetry operators do 140 i=1,numsymm do 138 j=1,3 fsym(j,4)=ts(j,i) do 138 k=1,3 138 fsym(j,k)=fs(k,j,i) if (i.lt.numsymm) write (2,141) maksym(fsym) if (i.eq.numsymm) write (2,142) maksym(fsym) 140 continue 141 format(2X,a) 142 format('1 ',a) C *** ATOM TYPES NSFAC=0 143 READ (3,1,END=144) LINE C C INTERPRET SFAC CARD(S) FOR ORDINAL NUMBER 'FEATURE' C IF(LINE(1:4).NE.'SFAC') GOTO 143 NSFAC=NSFAC+1 ATSYM(NSFAC)(1:2)=' ' DO 2000 II=6,LEN(LINE) CDEBUG WRITE(*,*)'???',LINE(II:II),'???' IF (LINE(II:II).NE.' ')THEN IF(ATSYM(NSFAC)(1:1).EQ.' ') THEN ATSYM(NSFAC)(1:1)=LINE(II:II) ELSE ATSYM(NSFAC)=ATSYM(NSFAC)(1:1)//LINE(II:II) ENDIF ELSE IF (ATSYM(NSFAC)(1:1).NE.' ')NSFAC=NSFAC+1 ATSYM(NSFAC)(1:2)=' ' ENDIF 2000 CONTINUE IF (ATSYM(NSFAC)(1:1).EQ.' ') NSFAC=NSFAC-1 GOTO 143 144 REWIND(3) CDEBUG write (*,*) 'NSFAC=',NSFAC,':', (ATSYM(III),III=1,NSFAC) DO 2100 II=1,NSFAC CDEBUG WRITE(*,*) '###',ATSYM(II),'###' DO 2200 JJ=1,103 IF(ATSYM(II).EQ.SYMBLS(JJ)) THEN ORDNL(II)=JJ CDEBUG WRITE(*,*) '===> ', ORDNL(II) GOTO 2100 ENDIF 2200 CONTINUE 2100 CONTINUE c *** free variables 145 read (3,1,end=999) line if (line(1:4) .ne. 'FVAR') go to 145 open(4,status='scratch') write (4,1) line(5:iend(line))//' 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.' backspace(4) read (4,*) fvar close(4) c *** atoms numatoms = 0 150 read (3,1,end=999) lab if ((lab .eq. 'AFIX') .or. (lab .eq. 'PART') .or. & (lab .eq. 'MOLE') .or. (lab .eq. 'SAME')) go to 150 if ((lab .eq. 'FMAP') .or. (lab .eq. 'BOND') .or. & (lab .eq. 'PLAN') .or. (lab .eq. 'CONF')) go to 150 if (lab .eq. 'REM ') goto 150 if ((lab .eq. 'HKLF') .or. (lab .eq. 'END ') .or. & lab .eq. ' ') then write (2,152) 152 format('ORIGIN',10x,'98',7x,'98', & ' 0.0 0.0 0.0',/,' 0.01000',53x,'7',/, & 'CENTER',10x,'99',7x,'99', & ' 0.5 0.5 0.5',/,'1 0.01000',53x,'7') norg = numatoms + 1 go to 195 end if 154 backspace(3) if (oldshx) then read (3,155) lab, itype, x, y, z, dum, U11, U22, char1 else read (3,156,err=157) lab, itype, x, y, z, dum, U11, U22, char1 endif 155 format(a4,i5,6f10.0,1x,a1) 156 format(a4,i5,3f12.0,3f11.0,a1) call fvarparm(0.,x,y,z) call fvarparm(U11,U22,0.,0.) numatoms = numatoms + 1 if (char1 .eq. '=') then if (oldshx) then read (3,158) U33, U23, U13, U12 else read (3,159) U33, U23, U13, U12 endif 158 format(5x,4f10.0) 159 format(5x,4f11.0) call fvarparm(U33,U12,U13,U23) write (2,160) lab,itype,x,y,z,U11,U22,U33,U12,U13,U23 160 format(a4,5x,i9,8x,'1',3f9.5,/,6f9.5,8x,'8') call calcUeq(ueq) else write (2,162) lab, itype, x, y, z, U11*8.*PI*PI 162 format(a4,5x,i9,8x,'2',3f9.5,/,f9.5,53x,'6') end if go to 150 195 continue c *** ORTEP-III instructions c *** 201 write (2,200) 200 format('# Initiate plotting',/,6x,'201') c *** 301 write (2,205) 205 format('# Drawing parameters') if (stereo) then write (2,210) 210 format(6x,'301',6x,'2.4',6x,'8.0',6x,'15.',6x,'.25') else write (2,215) 215 format(6x,'301',6x,'8.0',6x,'8.0',6x,'15.',6x,'.50') end if c *** 400 ATOMS Array write (2,220) -numatoms 220 format('# Include all input atoms in drawing',/, & 6x,'401',3x,'155501',i4,'55501') c *** 406 may be needed to complete structure write (2,222) numatoms, numatoms 222 format('# Complete structure with symmetry equivalent atoms, ', & 'if necessary',/,'# 406',2(8x,'1',i9),' 2.00') c *** vertices for cell outline, if desired write (2,223) norg, -norg 223 format('# Vertices for cell outline (add ', & 'appropriate lines to 1001 and 812)',/, & '# 401',i4,'55501',i4,'66601') c *** 506 write (2,225) 225 format('# Default origin and coordinate system',/,6x,'506') c *** 502 (commented out template) write (2,230) 230 format('# Use following instruction to change orientation',/, & '# 502',8x,'3',7x,'20',8x,'2',7x,'-5') c *** 604 write (2,235) 235 format('# Automatic scale and position',/,6x,'604') c *** 503 (stereo only) if (stereo) write (2,240) 240 format('# Stereo rotation for left eye view',/,6x,'503',8x,'2', & 6x,'2.7',/,'# Start save sequence',/,5x,'1101') c *** 1001 write (2,245) numatoms, numatoms 245 format('# Calculate overlap',/,2x,'2',2x,'1001',/,9x, & ' 1',i3,' 1',i3,' 3 0.9 1.52 .03') write (2,246) norg, norg, norg, norg, A-.01, A+.01, norg, norg, & norg, norg, B-.01, B+.01, norg, norg, norg, norg, C-.01, C+.01 246 format('# Add following lines to take cell outline into account', & /,'# 2',6x,4i3,' 1',2f6.2,' .01', & /,'# 2',6x,4i3,' 1',2f6.2,' .01', & /,'# ',6x,4i3,' 1',2f6.2,' .01') c *** 700 write (2,250) 250 format('# Draw anisotropic atoms [Feature #2 (col. 27) is "1"] ', & /,'# with open-octant ellipsoid',/, & 2x,'1',3x,'716',/,35x,'1',8x,'1',8x,'2') write (2,251) 251 format('# Use this instead to add labels', &' (.2" high, .5" to the right, .3" above)',/, & '#',1x,'1',3x,'716',36x,'0.2',6x,'0.5',6x,'0.3',/, & '#',34x,'1',8x,'1',8x,'2') write (2,255) 255 format('# Draw isotropic atoms [Feature #2 is "2"] ', & 'with outline only',/, & 2x,'1',3x,'714',/,35x,'2',8x,'2',8x,'2') c *** 800 C write (2,260) numatoms, numatoms C 260 format('# Draw bonds',/,2x,'2',3x,'812',/,9x,' 1',i3,' 1', C & i3,' 3 0.9 1.52 .03') C C write (2,260) 260 format('# Draw bonds by atom type [Feature #1] ',/,2x,'2',3x, &'812',17x,'1') do 261 i=1,nsfac do 261 j=i,nsfac write (2,262) atsym(i),atsym(j) 262 format('# Bonds between ',A2,' and ',A2) if (i+j .ne. nsfac*2) then write (2,263) i,i,j,j else write (2,264) i,i,j,j endif 263 format(2x,'2',6x,i3,i3,i3,i3,' 3 0.9 2.5 .03') 264 format(9x,i3,i3,i3,i3,' 3 0.9 2.50 .03') 261 continue C write (2,1260) 1260 format('# Add this to draw unit cell outline ',/,'#',1x,'2',3x, &'812',17x,'1') write (2,1261) norg, norg, norg, norg, A-.01, A+.01, norg, norg, & norg, norg, B-.01, B+.01, norg, norg, norg, norg, C-.01, C+.01 1261 format('# 2',6x,4i3,' 1',2f6.2,' .01', & /,'# 2',6x,4i3,' 1',2f6.2,' .01', & /,'# ',6x,4i3,' 1',2f6.2,' .01') C C c *** 1102, 202, 503, 1103 (stereo only) if (stereo) write (2,265) 265 format('# End save sequence',/,5x,'1102',/, & '# Shift plot origin',/,6x,'202',4x,'2.375',/, & '# Stereo rotation for right eye view',/,6x,'503',8x,'2', & 5x,'-2.7',/,'# Repeat save sequence',/,5x,'1103') c *** 202 and -1 write (2,270) 270 format('# Terminate plotting',/,6x,'202',/, & '# End of instructions',/,7x,'-1') endfile 2 close(2) write (6,300) 300 format(' TEP.IN created. Check file thoroughly for accuracy.') 999 stop 157 write (6,*)'Error reading atom cards. Retrying...' if (oldshx) then oldshx=.false. else oldshx=.true. endif goto 154 end subroutine calcUeq(ueq) c *** Code courtesy of Dr. Beverly R. Vincent, Molecular Structure Corp. common /CELL/ A, B, C, ALPHA, BETA, GAMMA common /Uij/ U11, U22, U33, U12, U13, U23 real c4,c5,c6,be,ga,ax,bx,cx,fct(6) RAD = 0.0174533 ! Degrees to Radians c4=cos(ALPHA * RAD) c5=cos(BETA * RAD) c6=cos(GAMMA * RAD) be=(c4*c6-c5)/((sqrt(1.0-c4*c4))*(sqrt(1.0-c6*c6))) ga=(c4*c5-c6)/((sqrt(1.0-c4*c4))*(sqrt(1.0-c5*c5))) ax=1.0/(A*(sqrt(1.0-c5*c5))*(sqrt(1.0-ga-ga))) bx=1.0/(B*(sqrt(1.0-c4*c4))*(sqrt(1.0-ga-ga))) cx=1.0/(C*(sqrt(1.0-c4*c4))*(sqrt(1.0-be*be))) fct(1)= A * A /3.0 fct(2)= B * B /3.0 fct(3)= C * C /3.0 fct(4)= 2.0* A * B*c6/3.0 fct(5)= 2.0* A * C*c5/3.0 fct(6)= 2.0* B * C*c4/3.0 ueq = U11*fct(1)*ax*ax + U22*fct(2)*bx*bx + * U33*fct(3)*cx*cx + U12*fct(4)*ax*bx + * U13*fct(5)*ax*cx + U23*fct(6)*bx*cx end subroutine addsymm(card) common /SYMM/ fs(3,3,96), ts(3,96), numsymm character*80 card character*24 sympart(3) numsymm=numsymm+1 ipart=1 do 100 jk=1,3 do 100 kl=1,24 100 sympart(jk)(kl:kl)=' ' jk=5 200 if (card(jk:jk).eq.' ') go to 500 lm=1 do 300 kl=jk,80 CMK270800 if(card(kl:kl).eq.' ' .or. card(kl:kl).eq.',') then if( card(kl:kl) .eq. ' '.and. kl.lt.80) go to 290 if(card(kl:kl).eq.',' .or.kl.eq.80) then jk=kl go to 400 end if sympart(ipart)(lm:lm)=card(kl:kl) lm=lm+1 290 continue 300 continue 400 ipart=ipart+1 500 jk=jk+1 if (jk.lt.80) go to 200 do 600 isymp=1,3 call tepsym(sympart(isymp),numsymm,isymp) 600 continue return end function iend(string) c *** returns position of last non-space character in string character string*(*) do 100 i=len(string),1,-1 if (string(i:i) .ne. ' ') then iend = i return end if 100 continue iend = 1 return end character*(*) function maksym(gp) c *** returns character string representation of symmetry operator dimension gp(3,4) character*1 xyz(3) character*5 fract(23) character*12 part(3) data fract/'1/24','1/12','1/8','1/6','5/24','1/4','7/24','1/3', * '3/8','5/12','11/24','1/2','13/24','7/12','5/8','2/3', * '17/24','3/4','19/24','5/6','7/8','11/12','23/24'/ data xyz/'x','y','z'/ do 500 i=1,3 part(i) = ' ' iff = 0 do 300 j=1,3 if (ifix(gp(i,j)) .ne. 0) then if (ifix(gp(i,j)) .eq. -1) * part(i) = part(i)(1:iend(part(i))) // '-' // xyz(j) if (ifix(gp(i,j)) .eq. 1 .and. iff .eq. 0) * part(i) = part(i)(1:iend(part(i))) // ' ' // xyz(j) if (ifix(gp(i,j)) .eq. 1 .and. iff .eq. 1) * part(i) = part(i)(1:iend(part(i))) // '+' // xyz(j) iff = 1 end if 300 continue gpval = gp(i,4) if (gpval.eq.1.) goto 500 if(gpval.gt..01 .or. gpval.lt.-.01) then if (gpval.lt.0.) then part(i) = part(i)(1:iend(part(i))) // '-' else part(i) = part(i)(1:iend(part(i))) // '+' end if gpval=abs(gpval) tfour=1./24. do 400 mm=1,23 tf=float(mm)*tfour if (gpval.gt.(tf-.01) .and. gpval.lt.(tf+.01)) iw=mm 400 continue part(i) = part(i)(1:iend(part(i))) // fract(iw) end if 500 continue maksym = part(1) // part(2) // part(3) return end subroutine tepsym(txt,num,kk) c *** parses character string representation of symmetry operators c *** and stores information common /SYMM/ fs(3,3,96), ts(3,96), numsymm character*24 txt logical twodig c *** convert txt to upper case do 101 i=1,24 iascii=ichar(txt(i:i)) if (iascii.ge.97.and.iascii.le.122) txt(i:i)=char(iascii-32) 101 continue c *** default value of ts if no fraction specified 202 ts(kk,num)=0. c *** look for and interpret a/b style fraction n=index(txt,'/') if (n.gt.0) then c get denominator k=ichar(txt(n+1:n+1))-48 m=ichar(txt(n+2:n+2))-48 if (m.ge.0 .and. m.le.9) then iden=k * 10 + m else iden=k end if c get numerator twodig=.false. ksign=1 k=ichar(txt(n-1:n-1))-48 if (n-2.ge.1) then m=ichar(txt(n-2:n-2))-48 if (m.ge.0 .and. m.le.9) twodig=.true. if (txt(n-2:n-2).eq.'-') ksign=-1 if (n-3.ge.1) then if (txt(n-3:n-3).eq.'-') ksign=-1 end if end if if (twodig) then inum=ksign * (m * 10 + k) else inum=ksign * k end if ts(kk,num)=float(inum)/float(iden) end if c *** look for and interpret decimal style fraction n=index(txt,'.') if (n.gt.0) then c get post decimal point portion k=ichar(txt(n+1:n+1))-48 m=ichar(txt(n+2:n+2))-48 if (m.ge.0 .and. m.le.9) then ts(kk,num)=float(k) * .1 + float(m) * .01 else ts(kk,num)=float(k) * .1 end if c get sign ksign=1 if (n-1.ge.1) then if (txt(n-1:n-1).eq.'-') ksign=-1 end if if (n-2.ge.1) then if (txt(n-2:n-2).eq.'-') ksign=-1 end if ts(kk,num)=float(ksign) * ts(kk,num) end if c *** interpret xyz portion of symmetry operation do 303 i=1,24 if (txt(i:i).eq.'X') then fs(1,kk,num)=1. if (i.ge.2) then if(txt(i-1:i-1).eq.'-') fs(1,kk,num)=-1. end if end if if (txt(i:i).eq.'Y') then fs(2,kk,num)=1. if (i.ge.2) then if(txt(i-1:i-1).eq.'-') fs(2,kk,num)=-1. end if end if if (txt(i:i).eq.'Z') then fs(3,kk,num)=1. if (i.ge.2) then if(txt(i-1:i-1).eq.'-') fs(3,kk,num)=-1. end if end if 303 continue return end subroutine addlatt(isymm,d1,d2,d3) c *** creates new symmetry operators from old by adding lattice c *** centering translations dimension dlatt(3) common /SYMM/ fs(3,3,96), ts(3,96), numsymm if (numsymm.ge.96) then write (6,200) return end if 200 format(' Maximum number of symmetry operators (96) exceeded.') dlatt(1)=d1 dlatt(2)=d2 dlatt(3)=d3 do 222 i=1,isymm inew=i+numsymm do 222 j=1,3 ts(j,inew)=ts(j,i)+dlatt(j) if (ts(j,inew).gt.1.) ts(j,inew)=ts(j,inew)-1. if (ts(j,inew).lt.-1.) ts(j,inew)=ts(j,inew)+1. do 222 k=1,3 fs(k,j,inew)=fs(k,j,i) 222 continue numsymm=isymm+numsymm return end subroutine fvarparm(a1,a2,a3,a4) c *** set atom parameters from FVAR values c *** if a1=0, a2, a3, and a4 are xyz positional parameters common /FVRS/ fvar(10),ueq real t(4) t(1)=a1 t(2)=a2 t(3)=a3 t(4)=a4 do 200 i=1,4 abst=abs(t(i)) if (a1.eq.0.) go to 100 c skip following if doing positional parameters if (t(i).lt.0. .and. abst.gt..5 .and. abst.lt.5.) then t(i)=ueq*abst go to 200 end if 100 iwhich=t(i)/10. if (iwhich.ge.0) then fact=t(i)-iwhich*10. else fact=abst-abs(iwhich)*10. end if if (fact.lt.5) then if (iwhich.eq.0) go to 200 if (iwhich.eq.1) t(i)=t(i)-10. if (iwhich.gt.1) t(i)=fvar(iwhich)*fact if (iwhich.lt.-1) t(i)=(1.-fvar(-iwhich))*fact end if 200 continue if (a1.ne.0.) a1=t(1) if (a2.ne.0.) a2=t(2) if (a3.ne.0.) a3=t(3) if (a4.ne.0.) a4=t(4) return end