[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ccp4bb]: A very clumsy solution to symmetry
*** For details on how to be removed from this list visit the ***
*** CCP4 home page http://www.dl.ac.uk/CCP/CCP4/main.html ***
Can you look at this abortion of a msymlb3?
At least it is returning SG I21 instead of I1211 for this f2mtz
I will attach the whole f2mtz.f
At least it works for this script..
#!/bin/csh -f
#
.../f2mtz HKLIN /y/people/reynolds/rfl/gregmyo3.rfl <<eof
title Convert Myoglobin from rfl to mtz
symmetry I21
cell 124.92 42.63 92.2 90.0 92.16 90.0
format '(f4.0,f4.0,f4.0,f10.2,x,f10.2)'
skipline 0
labout H K L FP PP
ctypout H H H F P
PNAME CCP4I_proj
DNAME
end
eof
.../f2mtz HKLIN /y/people/reynolds/rfl/gregmyo3.rfl <<eof
title Convert Myoglobin from rfl to mtz
symmetry 'I 1 21 1'
cell 124.92 42.63 92.2 90.0 92.16 90.0
format '(f4.0,f4.0,f4.0,f10.2,x,f10.2)'
skipline 0
labout H K L FP PP
ctypout H H H F P
PNAME CCP4I_proj
DNAME
end
eof
C*************************************************************
C
C This code is distributed under the terms and conditions of the
C CCP4 licence agreement as `Part ii)' software. See the conditions
C in the CCP4 manual for a copyright statement.
C
C**************************************************************
C
program f2mtz
* Formatted_2_MTZ
* Convert a formatted reflection file of some sort to mtz format.
* Written at the CCP4/ESF workshop in Cambridge, Oct-1991.
* Morten Kjeldgaard.
C
integer maxlab, maxtok, maxfrm
* Max number of labels
parameter (maxlab = 80)
* Max number of tokens for 'parser'
parameter (maxtok = 2*maxlab+1)
* Max length of user-specified format string
parameter (maxfrm = 200)
integer i, k, ct, lun, dummy, ntok, nskip
integer ibeg(maxtok), iend(maxtok), ityp(maxtok), idec(maxtok)
integer nsymx, nsympx, nspgrx, nzero
real fvalue(maxtok), cell(6), adata(maxtok-1),SCAL(MAXLAB)
character str*500,date*8,date6*6
character*64 PROJECTNAME,DATASETNAME
character*64 PNAME(maxlab),DNAME(maxlab)
character*4 key, cvalue(maxtok)
character spgrnx*10, pgnamx*10
character lsprgo(maxlab)*30, ctprgo(maxlab)
integer nlprgo
character frmat*200
real rsymx(4,4,192)
real nan
logical qisnan
external qnan,qisnan
logical lsymm, lfile, lcell, labout, lctyp, xpfreer(maxlab), lend
data lsymm, lfile, lcell, labout, lctyp /5*.false./
data lun /1/
data frmat /'*'/
data xpfreer /maxlab*.false./
ntok = maxtok
nsymx = 0
nskip = 0
DO 1 I = 1,MAXLAB
1 SCAL(I) = 1.0
call ccpfyp
call mtzini
CALL CCPRCS(6,'f2mtz','$Date: 1999/12/03 14:40:36 $')
write (*,*) 'Maximum number of columns in this version:',maxlab
C Set default PROJECTNAME and DATASETNAME
CALL CCPDAT(DATE)
WRITE(DATE6,'(A2,A2,A2)') DATE(1:2),DATE(4:5),DATE(7:8)
5 I = INDEX(DATE6,' ')
IF (I.GT.0) THEN
WRITE(DATE6(I:I),'(A1)') '0'
GOTO 5
ENDIF
PROJECTNAME = 'unknown'
DATASETNAME = 'unknown'//DATE6
* Open mtz file
call lwopen(1,'HKLOUT')
* Parse user input
10 continue
str = ' '
ntok = maxtok
call parser (key, str, ibeg, iend, ityp, fvalue, cvalue,
$ idec, ntok, lend, .true.)
if (lend) go to 500
* TITLE keyword
if (key .eq. 'TITL') then
call lwtitl (1, str(ibeg(2):iend(ntok)), 0)
goto 10
* CELL keyword - read in the cell parameters
else if (key .eq. 'CELL') then
call rdcell (2, ityp, fvalue, ntok, cell)
call lwcell(1, cell)
lcell = .true.
go to 10
* LABOUT keyword
else if (key .eq. 'LABO') then
k = 0
if (ntok.gt.maxlab) call ccperr (1, 'Too many labels')
do 20 i = 2, ntok
k = k+1
lsprgo(k) = str(ibeg(i):iend(i))
20 continue
nlprgo = (ntok-1)
write (*,*) 'Number of columns to be read in:', nlprgo
labout = .true.
goto 10
* CTYPOUT keyword
else if (key .eq. 'CTYP') then
k = 0
if (ntok.gt.maxlab) call ccperr (1, 'Too many columns')
do 22 i = 2, ntok
k = k+1
ctprgo(k) = str(ibeg(i):iend(i))
call ccpupc (ctprgo(k))
if (ctprgo(k).eq.'X') then
* fudge for X-PLOR free R input in opposite sense, indicated
C by column type X, which is converted to W. the column
C will have its value fixed if xpfreer is true
ctprgo(k) = 'I'
xpfreer(k) = .true.
endif
22 continue
lctyp = .true.
goto 10
C
* SCALE keyword
else if (key .eq. 'SCAL') then
k = 0
do 21 i = 2, ntok
k = k+1
SCAL(k) = fvalue(i)
21 continue
goto 10
* SYMM keyword
else if (key .eq. 'SYMM') then
call rdsymm (2, str, ibeg, iend, ityp, fvalue, ntok,
$ spgrnx, nspgrx, pgnamx, nsymx, nsympx, rsymx)
lsymm = .true.
write(6,'(a,2x,a,/,a,2x,a,/,3I5)')
+ ' spgrnx',spgrnx,'pgnamx',pgnamx, nspgrx, nsymx, nsympx
goto 10
* FILE keyword
else if (key .eq. 'FILE') then
str = str(ibeg(2):iend(2))
call ccpdpn(lun, str, 'READONLY', 'F', dummy, 0)
lfile = .true.
goto 10
* SKIP keyword
else if (key .eq. 'SKIP') then
if (ntok .le. 1) goto 10
if (ityp(2) .eq. 2) then
call keynum(1, 2, str, ibeg, iend, ityp, ntok)
nskip = int(fvalue(2))
write (*,*)'First', nskip,' lines of file will be skipped.'
endif
goto 10
* FORMAT keyword
else if (key .eq. 'FORM') then
if (ntok .ne. 2)
+ call fatal(lun,'FORMAT must have one string argument')
if (iend(2)-ibeg(2).gt.maxfrm) then
call fatal(lun,'FORMAT format string too long')
else
frmat = str(ibeg(2):iend(2))
endif
if (index(frmat,'i').ne.0 .or. index(frmat,'I').ne.0) then
write (*,FMT='(A,/,A)') ' *** Format contains an "I". ***',
+ ' Beware you cannot use I(or A) format. Replace with F?.0'
call ccperr (1, 'Bad keyword: '//key)
end if
if (index(frmat,'a').ne.0 .or. index(frmat,'A').ne.0) then
write (*,FMT='(A,/,A)') ' *** Format contains an "A". ***',
+ ' Beware you cannot use A(or I) format. Replace with F'
call ccperr (1, 'Bad keyword: '//key)
end if
goto 10
C
C ==========================================
C---- PNAME >> PROJECTNAME [pname] <string>
C ==========================================
C
C if given with DATASET then harvest will o/p a file
C this project_name should be always used for the
C one structure determination
C no default
C
ELSEIF (KEY.EQ.'PNAM') THEN
if (ntok .ge. 2) PROJECTNAME = str(ibeg(2):iend(2))
goto 10
C
C ===========================================
C---- DNAME >> DATASETNAME [dname] <string>
C ===========================================
C
C if given with PROJECT then harvest will o/p a file
C this dataset_name is the name of one of the diffraction
C data sets used in a particular project
C no default
C
ELSEIF (KEY.EQ.'DNAM') THEN
if (ntok .ge. 2) DATASETNAME = str(ibeg(2):iend(2))
goto 10
* END keyword
else if (key .eq. 'END') then
goto 500
* else error
else
call ccperr (1, 'Bad keyword: '//key)
endif
*** End of input parsing.
500 continue
if (.not. lcell) call fatal(lun,'CELL must be specified')
if (.not. lsymm) call fatal(lun,'SYMMETRY must be specified')
if (.not. labout) call fatal(lun,'LABOUT must be specified')
if (.not. lctyp) call fatal(lun,'CTYPOUT must be specified')
if (.not. lfile) call ccpdpn(lun, 'HKLIN', 'READONLY',
$ 'F', dummy, 0)
* Store the project name and dataset name in the mtz header:
if (PROJECTNAME.eq.'unknown')
+ call ccperr(2,' PROJECTNAME not assigned')
call lwid(1,PROJECTNAME,DATASETNAME)
* Store the column labels in the mtz header:
call lwclab(1, lsprgo, nlprgo, ctprgo, 0)
* Associate columns with datasets
do i=1,nlprgo
pname(i) = PROJECTNAME
dname(i) = DATASETNAME
enddo
call lwidas(1, nlprgo, pname, dname , 0)
* Store the symmetry in the mtz header:
call lwsymm (1, nsymx, nsympx, rsymx, spgrnx(1:1), nspgrx,
$ spgrnx, pgnamx)
* Store a history header:
call lwhstl (1, ' ')
*** Read all reflections from input file...
if (nskip .ne. 0) then
do 25 i = 1,nskip
read (lun,*)
25 continue
endif
ct = 0
call qnan(nan)
30 continue
if (frmat.eq.'*') then
do 32 i=1,nlprgo
C You can default data to `missing' with free format.
adata(i) = nan
32 continue
read (lun, *, end=50, err=1000) (adata(i),i=1,nlprgo)
else
read (lun, frmat, end=50, err=1000) (adata(i),i=1,nlprgo)
endif
* scale if required
nzero = 0
do 31 i = 1,nlprgo
if (qisnan(adata(i))) goto 31
if (adata(i) .eq. 0.0) nzero = nzero + 1
adata(i)= scal(i)*adata(i)
if (xpfreer(i)) then
* fiddle X-PLOR or SHELX free R to opposite sense
if (adata(i) .ge. 0.000001) then
* exclude
adata (i) = 0.0
else
* include
adata (i) = 1.0
endif
endif
31 continue
if (nzero .eq. nlprgo) then
write(*,*) ' ****Reflection rejected: record contains only ',
+ ' zeros (blank line?)****.'
else
ct = ct+1
call lwrefl (1, adata)
endif
goto 30
* Close mtz file, print full header:
50 call lwclos (1, 3)
close(unit=lun,status='keep')
write (*,'(i7, '' input reflections processed'')') ct
call ccperr(0,'Normal termination')
1000 write (*, FMT='(A,/,A)') '*** Read error',
+ 'Check FORMAT -- especially the need for REALs'
write (str, '(''problems reading reflection '',i7)') ct
call fatal (lun,str)
end
c=======================================================================
subroutine fatal (lun,str)
* Fatal -- issue fatal error message and die.
integer lun
character str*(*)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
close(unit=lun,status='keep')
call ccperr(1,str)
end
C =========================================================
SUBROUTINE MSYMLB3(IST,LSPGRP,NAMSPG_CIF,NAMSPG_CIFS,
+ NAMPG,NSYMP,NSYM,RlSymmMatrx)
C =========================================================
C
C---- Get symmetry operations for spacegroup LSPGRP from library file
C on stream IST, logical name SYMOP.
C
C In the library file, the header for each entry may contain
C Order not guaranteed, but must start with:
C LSPGRP NLINS and contain either NAMSPG or NAMSPG_CIF
C
C LSPGRP NLINS NLINP NAMSPG NAMPG CRYSTAL NAMSPG_CIF
C
C where LSPGRP spacegroup number
C NLINS total number of lines of symmetry operators.
C NLINP number of LINES of primitive symmetry operators
C Not used now..
C NAMSPG_CIF spacegroup name
C NAMPG name of corresponding pointgroup
C Not used now..
C
C On entry:
C IST stream number to read file
C LSPGRP spacegroup number
C NAMSPG or NAMSPG_CIF
C any acceptable spacegroup name: this will be used to
C identify the spacegroup if possible
C
C Returns
C LSPGRP spacegroup number
C NAMSPG_CIF full spacegroup name
C NAMSPG_CIFS name without any spaces
C NAMPG pointgroup name ( obtained from pgdefn - not 100% reliable!)
C NSYMP number of primitive symmetry operations obtained from pgdefn-
C only different from NSYM in non-primitive spacegroups
C NSYM total number of symmetry operations
C RlSymmMatrx(4,4,NSYM) Symmetry Rotation/translation matrices
C
C .. Parameters ..
INTEGER NPARSE
PARAMETER (NPARSE=200)
C ..
C .. Scalar Arguments ..
INTEGER IST,LSPGRP,NSYM,NSYMP
CHARACTER NAMPG* (*),NAMSPG_CIF* (*),NAMSPG_CIFS* (*),
+ NAMSAV*20,NAMSAV1*20
C ..
C .. Array Arguments ..
REAL RlSymmMatrx(4,4,*),ROTCHK(4,4)
C ..
C .. Local Scalars ..
INTEGER I,IFAIL,ISG,NLIN,NTOK
CHARACTER LINE*400,LINERR*400
C INTEGER NLINS
C ..
C .. Local Arrays ..
REAL FVALUE(NPARSE)
INTEGER IBEG(NPARSE),IDEC(NPARSE),IEND(NPARSE),ITYP(NPARSE),
+ IRCHK(192)
CHARACTER CVALUE(NPARSE)*4
LOGICAL NAMFIT
C ..
C .. External Subroutines ..
EXTERNAL CCPDPN,CCPUPC,PARSE,SYMFR2,LERROR
C ..
C .. Intrinsic Functions ..
INTRINSIC NINT
C ..
IFAIL = 0
CALL CCPDPN(IST,'SYMOP','READONLY','F',0,IFAIL)
C
NTOK = 0
NSYM = 0
C Remove all ' 1' then all remaining ' ' from long SG name
C
ILEN = LENSTR(NAMSPG_CIF)
NAMSPG_CIFS = ' '
C
NAMSAV = NAMSPG_CIF(1:1)
NAMSAV1= NAMSPG_CIF
write(6,'(A,/,3x,A8,/,3x,A8,/,3x,A8)')
+ ' 1 NAMSPG_CIF,NAMSPG_CIFS,NAMSAV',NAMSPG_CIF,NAMSPG_CIFS,NAMSAV
IF(ILEN.GE.2) THEN
J = 1
DO I = 2,ILEN
NAMSAV = NAMSAV(1:J)//NAMSAV1(I:I)
J = J + 1
IF( NAMSAV1(I:I+1).EQ.' 1') THEN
NAMSAV1(I+1:I+1) = ' '
END IF
write(6,'(A,i3,/,3x,A8,/,3x,A8)')
+ ' 2 NAMSAV1,NAMSAV',i,NAMSAV1,NAMSAV(1:J)
END DO
C
C Now search and destroy any remaining spaces
NAMSPG_CIFS = NAMSAV
J = 1
DO I = 2,ILEN
IF( NAMSAV(I:I).NE.' ') THEN
NAMSPG_CIFS = NAMSPG_CIFS(1:J)//NAMSAV(I:I)
J = J + 1
END IF
write(6,'(A,i3,/,3x,A8,/,3x,A8)')
+ ' 3 NAMSAV,NAMSPG_CIFS',i,NAMSAV,NAMSPG_CIFS(1:J)
END DO
END IF
C
write(6,'(A,/A,3x,A,3x,A)') '4 names',
+ NAMSPG_CIF, NAMSPG_CIFS,NAMSAV
10 CONTINUE
C
C---- Find correct space-group in file.
C Each space-group has header line of space-group number,
C number of line of symmetry operations for non-primitive
C and primitive cells.
C
READ (IST,FMT='(A)',ERR=30,END=30) LINE
CALL CCPUPC(LINE)
NTOK = -NPARSE
CALL PARSE(LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK)
C
C---- Compulsory Fields are space group number,
C number of lines,
C spacegroup name
C
IF (ITYP(1).NE.2 .OR. ITYP(2).NE.2 )
+ CALL LERROR(2,-1,'MSYMLB3: Error in format of SYMOP file: '
+ // LINE)
ISG = NINT(FVALUE(1))
NLIN = NINT(FVALUE(2))
c NLINS = NINT(FVALUE(3))
C
C---- Check for spacegroup name given - may be anywhere on line..
C---- Record the longest name on the line;
C---- This will be returned as NAMSPG_CIF
C
NAMLGTH = 1
NAMSAV = ' '
NAMFIT = .false.
C Set NAMSAV to longest name ana maybe find a match to NAMSPG_CIF..
DO 15 ITOK = 3,NTOK
C Spacegroup name must begin P A B C F I "H " R
IF (LINE(IBEG(ITOK):IBEG(ITOK)).NE.'P' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'A' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'B' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'C' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'F' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'I' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)+1).NE.'H ' .AND.
+ LINE(IBEG(ITOK):IBEG(ITOK)).NE.'R' ) GO to 15
C Ooo - get rid of CUBIC and any PG
IF(LINE(IBEG(ITOK):IBEG(ITOK)+1).EQ.'CU' .OR.
+ LINE(IBEG(ITOK):IBEG(ITOK)+1).EQ.'PG' ) GO to 15
C
LGTHCHK = MAX(NAMLGTH ,(IEND(ITOK)-IBEG(ITOK)+1))
IF(LGTHCHK .GT. NAMLGTH ) NAMSAV = LINE(IBEG(ITOK) :IEND(ITOK))
NAMLGTH = LGTHCHK
IF (NAMSPG_CIF.NE.LINE(IBEG(ITOK):IEND(ITOK)) .AND.
+ NAMSPG_CIFS.NE.LINE(IBEG(ITOK):IEND(ITOK))) GO TO 15
C Found a suitable match to space group name
NAMFIT = .true.
15 CONTINUE
C Move this - we need to check ALL tokens, not just till we find a match
IF(NAMFIT) GO TO 40
C
C
C---- No name match; check for spacegroup number if given
C
IF (LSPGRP.GT.0 .AND. LSPGRP.EQ.ISG) GO TO 40
C
C---- Not this one, skip NLIN lines
C
DO 20 I = 1,NLIN
READ (IST,FMT=*)
20 CONTINUE
C try again
GO TO 10
C
40 CONTINUE
C----- Reset space group name to longest on offer
LSPGRP = ISG
NAMSPG_CIF = NAMSAV
C
C Repeat the removal of all spaces from SG name
C
ILEN = LENSTR(NAMSPG_CIF)
NAMSPG_CIFS = ' '
C
NAMSAV = NAMSPG_CIF(1:1)
NAMSAV1= NAMSPG_CIF
write(6,'(A,/,3x,A8,/,3x,A8,/,3x,A8)')
+ ' 5 NAMSPG_CIF,NAMSPG_CIFS,NAMSAV',NAMSPG_CIF,NAMSPG_CIFS,NAMSAV
IF(ILEN.GE.2) THEN
J = 1
DO I = 2,ILEN
NAMSAV = NAMSAV(1:J)//NAMSAV1(I:I)
J = J + 1
IF( NAMSAV1(I:I+1).EQ.' 1') THEN
NAMSAV1(I+1:I+1) = ' '
END IF
write(6,'(A,i3,/,3x,A8,/,3x,A8)')
+ ' 6 NAMSAV1,NAMSAV',i,NAMSAV1,NAMSAV(1:J)
END DO
C
C Now search and destroy any remaining spaces
NAMSPG_CIFS = NAMSAV
J = 1
DO I = 2,ILEN
IF( NAMSAV(I:I).NE.' ') THEN
NAMSPG_CIFS = NAMSPG_CIFS(1:J)//NAMSAV(I:I)
J = J + 1
END IF
write(6,'(A,i3,/,3x,A8,/,3x,A8)')
+ ' 7 NAMSAV,NAMSPG_CIFS',i,NAMSAV,NAMSPG_CIFS(1:J)
END DO
END IF
C
write(6,'(A,/A,3x,A,3x,A)') '8 names',
+ NAMSPG_CIF, NAMSPG_CIFS,NAMSAV
C
C
C---- Space-group found, convert NLIN lines of
C symmetry operators to matrices
C
C read all sym ops at once; PGDEFN will sort out primitive and non primitive
DO 50 I = 1,NLIN
READ (IST,FMT='(A)') LINE
C Convert line to matrices
NSYM = NSYM + 1
CALL SYMFR2(LINE,1,NSYM,RlSymmMatrx)
50 CONTINUE
C
C-----Endeavor to test these sym ops form a closed group
C
do n = 1,nsym
irchk(n) = 0
end do
C
DO ISYM = 1,NSYM
C Determinant should be +1 or -1
CALL DETERM(DET,RlSymmMatrx(1,1,ISYM))
if(abs(det).lt.0.5) GO TO 25
DO JSYM = 1,NSYM
CALL MATMULNM(4,4,ROTCHK,RlSymmMatrx(1,1,ISYM),
+ RlSymmMatrx(1,1,JSYM))
C Check ROTCHK is also a symop
C---- Check This RSM Matrx for rotation and translation
C
IGOOD = 0
DO 90 N = 1,NSym
DO 95 I = 1,3
DO 98 J = 1,4
DCHK = ABS(ROTCHK(I,J) - RlSymmMatrx(I,J,N))
C
C---- This may be needed for translation components; no harm for others.
C
IF(J.EQ.4)
+ DCHK = ABS(MOD(ROTCHK(I,J) - RlSymmMatrx(I,J,N)
+ +99.5,1.0)-0.5)
IF (DCHK.LT.0.01) GO TO 98
C
C--- This MTZ symm op no good - off to check the next..
C
GO TO 90
98 CONTINUE
95 CONTINUE
C
C---- Found a good match - now check next ISM
C
IGOOD = 1
irchk(n) = irchk(n) + 1
GO to 80
90 Continue
C
C
C---- If this symmetry operator is missing no point going on. Try next SG
C
IF(IGOOD.EQ.0) THEN
Write(6,'(A,I4,A,I4,A)') ' Symm operators ', ISYM,' and',
+ JSYM,' have a problem.'
GO TO 35
END IF
80 CONTINUE
C
C
END DO
END DO
c
C
C Check all symmetry operators are equally likely to be generated
DO 100 N = 1,NSym
if(irchk(n) .ne.nsym ) go to 35
100 continue
C
call PGDEFN(NAMPG,NSYMP,NSYM,RlSymmMatrx,.FALSE.)
C
CLOSE (IST)
RETURN
C
25 CONTINUE
WRITE (LINERR,FMT='(A,A,I5,A)')
+ 'MSYLB3: Problem with sym op - determinant ne -+1',
+ ' space group number',LSPGRP,' in SYMOP file'
CALL LERROR(2,-1,LINERR)
C
30 CONTINUE
WRITE (LINERR,FMT='(A,A,I5,A)')
+ 'MSYLB3: No symmetry information for space group ',
+ ' number',LSPGRP,' in SYMOP file'
CALL LERROR(2,-1,LINERR)
C
35 CONTINUE
WRITE (LINERR,FMT='(A,A,I5,A)')
+ 'MSYLB3: Symmetry operators are not a closed group',
+ ' Something wrong for space group number',
+ LSPGRP,' in SYMOP file'
CALL LERROR(2,-1,LINERR)
END
C
C
------------------------------------------------------------------
Eleanor J.Dodson, Chemistry Department, University of York, U.K.
Tel: Home +44 (1904) 42 44 49, work: +44 (1904) 43 25 65
------------------------------------------------------------------