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 Reindex C C Reindex and/or change symmetry & reduce to asymmetric unit again C C The following operations are carried out on each reflection C 1) If unmerged data (multirecord file): reconstruct original C indices C 2) reindex if requested C 3) reject reflections with non-integral indices C 4) reduce reflection to asymmetric unit, using new symmetry (if C different). Symmetry cannot be changed for merged data. C C Currently, phase columns are not allowed. C Reflections with fractional indices after reindexing will be omitted C C Commands: C SYMMETRY new_spacegroup C REINDEX reindexing command C LEFTHAND allow hand inversion in reindexing C NOREDUCE don't reduce reindexed reflections to asymmetric unit C C C Note this program contains private copies of ASUSET & ASUPUT C (ASUSET1 & ASUPUT1) from symlib, so that they can have different C symmetry in common /RECSYM/ C C implicit none C integer mbatch, mcols parameter (mbatch=5000, mcols=500) INTEGER MSETS PARAMETER (MSETS=MCOLS) C real rsym(4,4,192) character spgnam*10, pgname*10, lattyp*1 integer numsgp, nsym, nsymp, nlaue real rsym1(4,4,192) character spgnam1*10, pgname1*10, lattyp1*1 integer numsgp1, nsym1, nsymp1, nlaue1 external ccprcs, lwhstl, matmul C C----Harvesting stuff C INTEGER NDATASETS,ISETS(MSETS),CSETID(MCOLS) REAL DCELL(6,MSETS),DWAVE(MSETS) CHARACTER*64 PNAME(MSETS),DNAME(MSETS) C C C----Column matching stuff C integer ncolgf(mcols),csetidgf(mcols),nset(mcols), + ncolql(mcols),csetidql(mcols), + ncolki(mcols),csetidki(mcols) integer nanogf,nanoql,nanoki logical lanogf,lanoql,lanoki C Things for PARSER ---- integer maxtok parameter (maxtok=100) character key*4,cvalue(maxtok)*4,line*200,linemat*200,linemat2*200 integer ibeg(maxtok),iend(maxtok),ityp(maxtok), . idec(maxtok),ntok real fvalue(maxtok) logical lend C C Reindex real hklmat(3,3), hklinv(3,3), hkltrn(3), + tmpmat(4,4), tmpmat1(3,3),tmpmat2(3,3),tmpmat3(3,3) integer ibeglhs(3),iendlhs(3),ibegrhs(3),iendrhs(3) logical lrindx, lhkltr C C Batch stuff integer nbatch, jbatch(mbatch) C C Locals integer mflin, mflout, iprint, ifail, ihkl(3), jhkl(3), $ isym, m, ierr, i, j, nrfout, nrjfrc, isort(5), negadf, ic real s, data(500), umat(3,3), umat2(3,3), $ cell(6), cell2(6), cell3(6), det, hkl(3), v(3) logical eof, lhand, lresym, lphase, lreduc, lphcof, lanodf C C integer maxcol parameter (maxcol=4) character lsprgi(maxcol)*30, ctprgi(maxcol)*1, $ ctyps(mcols)*1, clabs(mcols)*30 integer lookup(maxcol), nlprgi, ncmsym, ncol logical logmss (mcols) C C M/ISYM column is optional data nlprgi/maxcol/ data lsprgi/'H','K','L','M/ISYM'/ data ctprgi/'H','H','H','Y'/ data lookup/-1, -1, -1, 1/ C C call ccpfyp call ccprcs (6, 'REINDEX', '$Date: 2002/08/06 13:05:28 $') call mtzini C nsym1 = 0 lrindx = .false. lhkltr = .false. lhand = .false. lresym = .false. lreduc = .true. mflin = 1 mflout = mflin iprint = 1 ifail = 0 C C do icol = 1,mcols ncolgf(icol) = 0 ncolql(icol) = 0 ncolki(icol) = 0 nset(icol) = 0 C If these aren't got from LRCLID, all columns will be in dataset 0 csetid(icol) = 0 c Actually these are all unused at present.. csetidgf(icol) = 0 csetidql(icol) = 0 csetidki(icol) = 0 end do C ierr = 0 C 10 line=' ' ntok=maxtok call parser( . key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.) if(lend) go to 1000 if(ntok.eq.0) goto 10 call ccpupc(key) C----------------------------------------------------------------------- C SYMMETRY command; if (key .eq. 'SYMM') then call rdsymm(2,line,ibeg,iend,ityp,fvalue,ntok, . spgnam1,numsgp1,pgname1,nsym1,nsymp1,rsym1) lresym = .true. C----------------------------------------------------------------------- C REINDEX command Followed by HKL or AXIS elseif (key .eq. 'REIN') then C call ccpupc(line) IF(LINE(ibeg(2):ibeg(2)) .NE.'A') then C Check if the keyword HKL is given itok = 2 IF(LINE(ibeg(2):iend(2)) .EQ.'HKL') itok = 3 C symfr2a returns a reindexing matrix [H] which postmultiplies row vector C ie [h' k' l'] = [h ] [H] C [k ] C [l ] C C Set hklmat = HTransp such that [h'] = [HTransp] [h] C [k'] [k] C [l'] [l] C lrindx = .true. j = 1 call symfr2(line,ibeg(itok),j,tmpmat) if (j .ne. 1) then write (6, '(a)') ' *** Illegal REINDEX command ***' ierr = ierr + 1 else C copy to 3x3 matrix & translations lhkltr = .false. do 101, j = 1,3 hkltrn(j) = tmpmat(j,4) if (abs(hkltrn(j)) .gt. 0.05) lhkltr = .true. do 102, i = 1,3 hklmat(i,j) = tmpmat(i,j) 102 continue 101 continue endif END IF C----------------------------------------------------------------------- C REINDEX command IF(LINE(ibeg(2):iend(2)) .EQ.'AXIS') then C We need to strip out "*" symbols if present - symfr2 cant cope with them.. linemat = line(ibeg(1):iend(ntok)) lenmat = lenstr(linemat) linemat2 = line(ibeg(1):iend(2)) lenmat2 = lenstr(linemat2) +1 linemat2 = linemat2//' ' l2 = lenmat2 do l=iend(2)+1,lenmat if(linemat(l:l) .ne.'*' ) then linemat2 = linemat2(1:l2)//linemat(l:l) l2 = l2 + 1 end if end do call parser( . key,linemat2,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.) itok = 3 C C symfr2 returns a conversion matrix [M] C [a*new] = [M] [a*old] C [b*new] [b*old] C [c*new] [c*old] C C Since [hn kn ln] [a*n] = [ho ko lo] [a*o] C [b*n] [b*o] C [c*n] [c*o] C C [hn kn ln] [M] [a*o] = [ho ko lo] [a*o] C [b*o] [b*o] C [c*o] [c*o] C C and [hn kn ln] = [ho ko lo] [M]**-1 C lrindx = .true. j = 1 call symfr2(linemat2,ibeg(itok),j,tmpmat) if (j .ne. 1) then write (6, '(a)') ' *** Illegal REINDEX command ***' ierr = ierr + 1 else C copy to 3x3 matrix & translations lhkltr = .false. C Actually need the transpose of tmpmat**-1 - se above do j = 1,3 hkltrn(j) = tmpmat(j,4) if (abs(hkltrn(j)) .gt. 0.05) lhkltr = .true. do i = 1,3 tmpmat1(j,i) = tmpmat(i,j) end do end do call minv3(hklmat,tmpmat1,det) endif END IF elseif (key .eq. 'MATC') then C C Give 3 matches - from this the program can work out the C correct reindexing matrix. C Example: 2h' = h -k +l * 2k' = -h +k + l * 1l' = 1h' C This means that the reindexing matrix [H] satisfies these equations: C C [ 2 0 0 ] [ 1 1 -1] [ H ] C [ 0 2 0 ] = [-1 1 1] C [ 0 0 1 ] [ 1 0 0] C C C so [H] = [ 1 1 -1] **(-1) [ 2 0 0] C [-1 1 1] [ 0 2 0] C [ 1 0 0] [ 0 0 1] C C = [ 0.0 0.0 1.0] [ 2 0 0] C [ 0.5 0.5 0.0] [ 0 2 0] C [-0.5 0.5 1.0] [ 0 0 1] C C = [ 0 0 1] C [ 1 1 0] C [-1 1 1] C C C C what was this dooing here!!!? c lmatch = .true. j = 1 C Check there are just 3 matches specified. C Count "=" signs (3) and "*" signs (2) linemat = line(ibeg(2):iend(ntok)) lenmat = lenstr(linemat) neqsign = 0 nstar = 0 ibeglhs(1) =1 iendrhs(3) =lenmat do l=1,lenmat C Get rid of primes if(linemat(l:l) .eq.'''') linemat(l:l)=' ' if(linemat(l:l) .eq.'=') then neqsign = neqsign+1 iendlhs(neqsign) =l-1 ibegrhs(neqsign) =l+1 end if if(linemat(l:l) .eq.'*' .or. linemat(l:l) .eq.';') then nstar = nstar+1 iendrhs(nstar) =l-1 ibeglhs(nstar+1) =l+1 end if end do if( neqsign.ne.3 .or. nstar.ne.2 ) then write (6, '(a)') ' *** Illegal MATCH command ***' ierr = ierr + 1 end if linemat2 = linemat(ibeglhs(1):iendlhs(1))//','// + linemat(ibeglhs(2):iendlhs(2))//','// + linemat(ibeglhs(3):iendlhs(3)) call symfr2(linemat2,1,j,tmpmat) if (j .ne. 1) then write (6, '(a)') ' *** Illegal MATCH command ***' ierr = ierr + 1 else C copy to 3x3 matrix do 103, j = 1,3 hkltrn(j) = tmpmat(j,4) if (abs(hkltrn(j)) .gt. 0.05) lhkltr = .true. do 104, i = 1,3 tmpmat1(i,j) = tmpmat(i,j) 104 continue 103 continue endif write(6,'(//,a ,3(/,3f18.5))') + ' Matrix LHS', + ((tmpmat(i,j),j=1,3),i=1,3) linemat2 = linemat(ibegrhs(1):iendrhs(1))//','// + linemat(ibegrhs(2):iendrhs(2))//','// + linemat(ibegrhs(3):iendrhs(3)) j = 1 call symfr2(linemat2,1,j,tmpmat) if (j .ne. 1) then write (6, '(a)') ' *** Illegal MATCH command ***' ierr = ierr + 1 else C copy to 3x3 matrix do 105, j = 1,3 hkltrn(j) = tmpmat(j,4) if (abs(hkltrn(j)) .gt. 0.05) lhkltr = .true. do 106, i = 1,3 tmpmat2(i,j) = tmpmat(i,j) 106 continue 105 continue endif write(6,'(//,a ,3(/,3f18.5))') + ' Matrix RHS', + ((tmpmat2(i,j),j=1,3),i=1,3) call minv3(tmpmat3,tmpmat2,det) call matmul(tmpmat2,tmpmat3,tmpmat1) call transp(hklmat,tmpmat2) lrindx = .true. C----------------------------------------------------------------------- C LEFTHAND command, allow lefthanded matrix (why?) elseif (key .eq. 'LEFT') then lhand = .true. C----------------------------------------------------------------------- C NOREDUCE don't reduce reflections to an asymmetric unit elseif (key .eq. 'NORE') then lreduc = .false. C----------------------------------------------------------------------- elseif (key .eq. 'END') then go to 1000 C----------------------------------------------------------------------- else write (6 ,'(/a/)') ' **** Unrecognized command ****' ierr = ierr+1 endif go to 10 C----------------------------------------------------------------------- C----------------------------------------------------------------------- 1000 if (ierr .gt. 0) then call ccperr(1, '*** Error in control input ***') endif C if (lrindx) then C Invert matrix call minv3(hklinv,hklmat,det) C if (abs(det) .lt. 0.00001) then write (6, '(/a/)') $ ' **** Reindexing matrix has zero determinant ****' call ccperr(1,'** Illegal reindexing **') elseif (det .lt. 0.0) then call ccperr(2, $ ' !!!! Reindexing matrix INVERTS hand !!!!') if (.not.lhand) then call ccperr(1, $ ' !!!! You are NOT allowed to do this - Changing all signs $ in reindexing matrix') do j = 1,3 do i = 1,3 hklmat(i,j) = -hklmat(i,j) end do end do else write (6,'(/A)') $ ' !!!! Do you REALLY mean to do this ? !!!!' endif endif C write (6,'(/A)') $ ' Reflections will be reindexed, and unit cell recalculated' C write (6, 6101) ((hklmat(i,j), i=1,3), j=1,3) 6101 format(/' Reindexing transformation:'/ $ 5x,' h'' = ','( h k l ) (',3f9.5,' )'/ $ 2(27x,'(',3f9.5,' )'/)/) C write (6, 6111) ((hklmat(i,j), i=1,3), j=1,3) 6111 format(/' Real axes transformed by same matrix:'/ $ 5x,' a'' = ','( a b c ) (',3f9.5,' )'/ $ 2(27x,'(',3f9.5,' )'/)/) C C Invert to get transform for reciprocal axes and coordinates. call minv3(tmpmat1,hklmat,det) C write (6, 6121) ((tmpmat1(i,j), i=1,3), j=1,3) 6121 format(/' Reciprocal axes transformed by inverse matrix:'/, $ 5x,' a*'' = ',' (',3f9.5,' ) ( a*)', /, $ 27x,'(',3f9.5,' ) ( b*)', /, $ 27x,'(',3f9.5,' ) ( c*)'/) C C write (6, 6131) ((tmpmat1(i,j), i=1,3), j=1,3) 6131 format(/' Real coordinates transformed by same matrix:'/, $ 5x,' x'' = ',' (',3f9.5,' ) ( x)', /, $ 27x,'(',3f9.5,' ) ( y)', /, $ 27x,'(',3f9.5,' ) ( z)'/) C C if (lhkltr) write (6,'(/A/A/A,3F8.1/)') $ ' !!!! WARNING - you are translating the indices !!!!', $ ' !!!! orientation data may be in error !!!!', $ ' hkl translation vector : ',hkltrn C endif C call lropen(mflin, 'HKLIN', iprint, ifail) if (ifail .ne. 0) then go to 999 endif C C---- Open input file call lrassn(mflin, lsprgi, nlprgi, lookup, ctprgi) C Column for M/ISYM, if present ncmsym = lookup(4) call lrsymi(mflin, nsymp, lattyp, numsgp, spgnam, pgname) call lrsymm(mflin, nsym, rsym) call lrcell(mflin, cell) C Read batch information, if present call lrbats(mflin ,nbatch,jbatch) C If SYMMETRY changed; must also turn off reduction to asymmetric unit if(nbatch.eq.0 .and.lresym) lreduc = .false. C C---- Read column types (& labels) call lrclab(mflin, clabs, ctyps, ncol) C C Read project and dat set names CALL LRIDC(MFLIN,PNAME,DNAME,ISETS,DCELL,DWAVE,NDATASETS) C---- Get dataset ID for column and match to dataset header info. IF (NDATASETS.GT.0) + CALL LRCLID(MFLIN,CSETID,NCOL) C C---- Check for phase & anomalous difference columns lphase = .false. lphcof = .false. lanodf = .false. C---- Check for F+ F-, or I+ I- pairs nanogf = 0 lanogf = .false. nanoql = 0 lanoql = .false. nanoki = 0 lanoki = .false. C do 50, i = 1,ncol if (ctyps(i) .eq. 'P') then C Phase lphase = .true. elseif (ctyps(i) .eq. 'A') then C Phase coefficients lphcof = .true. elseif (ctyps(i) .eq. 'D') then C Anomalous difference lanodf = .true. elseif (ctyps(i) .eq. 'G') then C F+ F- pair; store column number ( labelled G) lanogf = .true. nanogf = nanogf + 1 ncolgf(nanogf) = i csetidgf(nanogf) = csetid(i) elseif (ctyps(i) .eq. 'L'.or.ctyps(i) .eq. 'M') then C Q+ Q- pair; store column number ( labelled L (SIGF) or M(SIGI)) lanoql = .true. nanoql = nanoql + 1 ncolql(nanoql) = i csetidql(nanoql) = csetid(i) elseif (ctyps(i) .eq. 'K') then C I+ I- pair; store column number lanoki = .true. nanoki = nanoki + 1 ncolki(nanoki) = i csetidki(nanoki) = csetid(i) endif 50 continue C C--- Check there are pairs of ncolgf and ncolki for the same dataset ichk = mod(nanogf,2) if(ichk .ne. 0) write (6, '(/a/)') $ ' **** Something WRONG with F+ F- pairing ****' if(ichk .ne. 0) call ccperr(1,' **** Check data set ****') C ichk = mod(nanoql,2) if(ichk .ne. 0) write (6, '(/a/)') $ ' **** Something WRONG with Q+ Q- pairing ****' if(ichk .ne. 0) call ccperr(1,' **** Check data set ****') C ichk = mod(nanoki,2) if(ichk .ne. 0) write (6, '(/a/)') $ ' **** Something WRONG with I+ I- pairing ****' if(ichk .ne. 0) call ccperr(1,' **** Check data set ****') C C--- Match partners do icset = 1,nanogf if(nset(ncolgf(icset)) .eq. 0) then do jcset = icset+1,nanogf if(csetid(ncolgf(icset)) .eq. csetid(ncolgf(jcset)) ) then nset(ncolgf(icset)) = ncolgf(jcset) nset(ncolgf(jcset)) = ncolgf(icset) go to 52 end if end do end if 52 continue C Has this worked? if(nset(ncolgf(icset)) .eq. 0) then write (6, '(/a,i3,a/)') + ' **** No match for column ncolgf ',icset,' ****' call ccperr(1,' **** Check data set ****') else write(6, '(/a,i3,3x,i3/)') + ' Column pairs- F+ F-:',ncolgf(icset), nset(ncolgf(icset)) end if end do do icset = 1,nanoql if(nset(ncolql(icset)) .eq. 0) then do jcset = icset+1,nanoql if(csetid(ncolql(icset)) .eq. csetid(ncolql(jcset)) .and. + ctyps(ncolql(icset)) .eq. ctyps(ncolql(icset))) then nset(ncolql(icset)) = ncolql(jcset) nset(ncolql(jcset)) = ncolql(icset) go to 54 end if end do end if 54 continue C Has this worked? if(nset(ncolql(icset)) .eq. 0) then write (6, '(/a,i3,a/)') + ' **** No match for column ncolql ',icset,' ****' call ccperr(1,' **** Check data set ****') else write(6, '(/a,i3,3x,i3/)') + ' Column pairs - SIGs:',ncolql(icset), nset(ncolql(icset)) end if end do do icset = 1,nanoki if(nset(ncolki(icset)) .eq. 0) then do jcset = icset+1,nanoki if(csetid(ncolki(icset)) .eq. csetid(ncolki(jcset)) ) then nset(ncolki(icset)) = ncolki(jcset) nset(ncolki(jcset)) = ncolki(icset) go to 56 end if end do end if 56 continue C Has this worked? if(nset(ncolki(icset)) .eq. 0) then write (6, '(/a,i3,a/)') + ' **** No match for column ncolki ',icset,' ****' call ccperr(1,' **** Check data set ****') else write(6, '(/a,i3,3x,i3/)') + ' Column pairs - I+ I-:',ncolki(icset), nset(ncolki(icset)) end if end do C if (lphase .or. lphcof) then write (6, '(//1x,18a/)') ('****',i=1,18) write (6, '(/a,a/)') $ ' **** This program cannot handle properly phase or ', $ 'phase coefficient columns - You will need to ', $ 'recalculate them ****' if(.not. lresym) write (6, '(/a/)') $ ' **** These columns will probably be WRONG, so beware' write (6, '(//1x,18a/)') ('****',i=1,18) if(lresym) write (6, '(/a/)') $ ' **** These columns will CERTAINLY BE WRONG, so beware' write (6, '(//1x,18a/)') ('****',i=1,18) endif C if (lanodf .and. lreduc) then write (6, '(/a,a/)') $ ' Anomalous difference values will be negated if hand', $ ' of indices hkl is changed in reducing to asymmetric unit' endif C if (lanogf .and. lreduc) then write (6, '(/a,a/)') $ ' Anomalous pairs will be swapped if hand', $ ' of indices hkl is changed in reducing to asymmetric unit' endif C C---- Process symmetry C Set up ASUGET with input symmetry write (6, '(/a/)') ' Input symmetry:' call asuset(spgnam, numsgp, pgname, nsym, rsym, nsymp, nlaue, $ .true.) call centric(nsym,rsym,0) if (lresym) then C Set up ASUPUT1 with new symmetry write (6, '(/a/)') ' Output symmetry:' call asuset1( $ spgnam1, numsgp1, pgname1, nsym1, rsym1, nsymp1, nlaue1, $ .true.) call centric(nsym1,rsym1,0) endif C if (lreduc) then write (6,'(/a/)') $ ' Reflection indices will be reduced to the asymmetric unit' else C--- Noreduce option if (nbatch .eq. 0) then write (6,'(/a/)') $ ' Reflection indices will NOT be reduced to the asymmetric unit' else write (6,'(/a/)') $ ' **** Unmerged data MUST be reduced to the asymmetric unit ****' call ccperr(1,'** Illegal NOREDUCE option **') endif endif C call lwopen(mflout, 'HKLOUT') call lwhstl (mflout,' ') if (lresym) then if (nbatch .eq. 0) then C Warn if change of symmetry with merged data write (6, '(///a,a/)') $ ' !!!! You are changing the symmetry of', $ ' merged data are you SURE you know what you are doing!!!!' call ccperr(2,'** Symmetry change of merged data **') endif lattyp1 = spgnam1(1:1) call lwsymm(mflout, $ nsym1, nsymp1, rsym1, lattyp1, numsgp1, spgnam1, pgname1) endif C C Flag file as unsorted do 80, i=1,5 isort(i) = 0 80 continue call lwsort(mflout, isort) C if (lrindx) then call unitmt(umat,3) C Update all dataset cells do i = 1,NDATASETS do j = 1,6 cell3(j) = dcell(j,i) enddo call newcel(hklmat, hklinv, cell3, umat, cell2, umat2, ierr) if (ierr .ne. 0) call ccperr(1, '*** ERROR in new cell ***') call LWIDC(mflout,PNAME(i),DNAME(i),cell2,DWAVE(i)) enddo C Update all global cell call newcel(hklmat, hklinv, cell, umat, cell2, umat2, ierr) if (ierr .ne. 0) then call ccperr(1, '*** ERROR in new cell ***') endif call lwcell(mflout, cell2) write (6,'(/A,6F8.2/)') $ ' New unit cell determined from reindexing: ',cell2 if (nbatch .gt. 0) then C Batch headers present, change cell etc in all of them call modbhd(mflin, mflout, nbatch, jbatch, hklmat, hklinv) endif endif C C************************************************************************** C************************************************************************** nrjfrc = 0 nrfout = 0 negadf = 0 C 100 call lrrefl(mflin, s, data, eof) if (eof) go to 200 call lrrefm (mflin, logmss) ihkl(1) = nint(data(lookup(1))) ihkl(2) = nint(data(lookup(2))) ihkl(3) = nint(data(lookup(3))) C--- Pick up "original" indices if (ncmsym .gt. 0) then C Assume this datum isn't missing if the column is assigned isym = nint(data(ncmsym)) m = isym/256 isym = isym - m * 256 call asuget(ihkl, jhkl, isym) else jhkl(1) = ihkl(1) jhkl(2) = ihkl(2) jhkl(3) = ihkl(3) endif C C--- Reindex if (lrindx) then C determine new indices do 130, i = 1,3 hkl(i) = jhkl(i) 130 continue call matvec(v,hklmat,hkl) call vsum(v,v,hkltrn) C We may now have fractional indices: such reflections are rejected do 140, i = 1,3 j = nint(v(i)) if (abs(v(i)-float(j)) .gt. 0.05) then C Reject nrjfrc = nrjfrc + 1 go to 100 else jhkl(i) = j endif 140 continue C endif C C---- Put reflection back into asymmetric unit if (lreduc) then if (lresym) then call asuput1(jhkl, ihkl, isym) else call asuput(jhkl, ihkl, isym) endif C Change sign of any anomalous difference columns, acentric reflections only if (lanodf) then if (mod(isym,2) .eq. 0) then call centr(ihkl, ic) if (ic .eq. 0) then negadf = negadf + 1 if (negadf .lt. 20) then write (6,'(a,3i5)') + ' Negate DelAno and/or swap F,I,Q columns',ihkl elseif (negadf .eq. 20) then write (6,'(a/)') ' . . . more negated . . .' endif do 150, i = 1,ncol if (ctyps(i) .eq. 'D' .and. .not.logmss(i)) then data(i) = -data(i) endif c if (ctyps(i) .eq. 'G' .and. .not.logmss(i)) then C + - pairs to swap if (ctyps(i) .eq. 'G' .or.ctyps(i) .eq. 'L' .or. + ctyps(i) .eq. 'K' .or.ctyps(i) .eq. 'M') then if(nset(i) .gt. i) then dsave = data(i) data(i) = data(nset(i)) data(nset(i)) = dsave end if c if(.not.logmss(i) .and. .not.logmss(nset(i))) c + write(6,*)data(i),data(nset(i)) endif 150 continue endif endif endif else ihkl(1) = jhkl(1) ihkl(2) = jhkl(2) ihkl(3) = jhkl(3) endif data(lookup(1)) = ihkl(1) data(lookup(2)) = ihkl(2) data(lookup(3)) = ihkl(3) if (ncmsym .gt. 0) then data(ncmsym) = m * 256 + isym endif C nrfout = nrfout+1 call lwrefl(mflout, data) go to 100 C 200 call lwclos(mflout, iprint) write (6, '(/a,i10/)') $ ' Number of reflections written to output file:',nrfout if (nrjfrc .gt. 0) then write (6, '(a,i10/)') $ ' Rejected reflections with fractional indices:',nrjfrc endif if (negadf .gt. 0) then write (6, '(a,i10/)') $ ' Reflections with negated anomalous difference:',negadf endif C call ccperr(0, ' ** Normal termination') C 999 call ccperr(1, ' ** Error in input **') end C C C ASUSET1, ASUPUT1 are private versions of routines which share C coomn block /RECSYM1/ instead of /RECSYM/ C subroutine asuset1( . spgnam,numsgp,pgname,msym,rrsym,msymp,mlaue,lprint) C ======================================================== C C Set up & store symmetry for later use in ASUPUT or ASUGET C C On input: C SPGNAM space-group name (not used) C NUMSGP space-group number (not used) C PGNAME point-group name (if returned from SYMOP.LIB) C MSYM total number of symmetry operations C RRSYM(4,4,MSYM) symmetry matrices (real-space) C C On output: C PGNAME point-group name C MSYMP number of primitive symmetry operations C MLAUE Laue group number C C C Arguments: integer numsgp, msym, msymp, mlaue real rrsym(4,4,192) CcMSYM) character*(*) spgnam, pgname logical lprint C C Common blocks C C Common block /RECSYM1/ C RSYM real-space symmetry operators C RSYMIV their inverse C NSYM number of symmetry operations C NSYMP number of primitive symmetry operations C NLAUE number of Laue group integer maxsym parameter (maxsym=192) common /recsym1/rsym(4,4,maxsym),rsymiv(4,4,maxsym), . nsym,nsymp,nlaue integer nsym,nsymp,nlaue real rsym,rsymiv save /recsym1/ C C Functions integer lenstr C Locals integer i, j, l character*100 strout, name*8, launam*8 C C Get point group and primitive-only operations call pgdefn(pgname,msymp,msym,rrsym,.false.) C C Get Laue group number call pgnlau(pgname,mlaue,launam) C C Store in common block nsym = msym nsymp = msymp nlaue = mlaue do 10, i=1,nsym call ccpmvi(rsym(1,1,i), rrsym(1,1,i), 16) C Invert all symmetry operators call invsym(rsym(1,1,i), rsymiv(1,1,i)) 10 continue C C Print if (lprint) then call blank('CURWIN',1) strout = ' Reciprocal space symmetry' call putlin(strout,'CURWIN') name = pgname l = lenstr(name) if (name(1:2) .eq. 'PG') name = name(3:l) l = lenstr(name) i = lenstr(spgnam) j = lenstr(launam) write (strout, 6001) . spgnam(1:i),numsgp,name(1:l),launam(1:j) 6001 format(' Space group: ',a,' (',i3,')',5x, . 'Point group: ',a,5x,'Laue group: ',a) call putlin(strout,'CURWIN') C if (nlaue .eq. 3) then strout = '[-1] hkl:l>=0 hk0:h>=0 0k0:k>=0' elseif (nlaue .eq. 4) then strout = '[2/m] hkl:k>=0, l>=0 hk0:h>=0' elseif (nlaue .eq. 6) then strout = '[mmm] hkl:h>=0, k>=0, l>=0' elseif (nlaue .eq. 7) then strout = '[4/m] hkl:h>=0, l>=0 with k>=0 if h=0'// . ' and k>0 if h>0' elseif (nlaue .eq. 8) then strout = '[4/mmm] hkl:h>=0, k>=0, l>=0 and h>=k' elseif (nlaue .eq. 9) then strout = '[-3] hkl:h>=0, k>0 00l:l>0' elseif (nlaue .eq. 10) then strout = '[312] hkl:h>=0, k>=0 with k<=h '// . 'for all l. If h = 0 l>=0' elseif (nlaue .eq. 11) then strout = '[321] hkl:h>=0, k>=0 with k<=h '// . 'for all l. If h = k l>=0' elseif (nlaue .eq. 12) then strout = '[6/m] hkl:h>=0, k>=0, l>=0 with k>=0'// . ' if h=0, and k> 0 if h>0' elseif (nlaue .eq. 13) then strout = '[6/mmm] hkl:h>=0, k>=0, l>=0 with h>=k' elseif (nlaue .eq. 14) then strout = '[m3] hkl:h>=0, k>=0, l>=0 with l>=h,'// . ' k>=h if l=h, k> h if l>h' elseif (nlaue .eq. 15) then strout = . '[m3m] hkl:h>=0, k>=0, l>=0 with k>=l, and l>=h' else write(strout, 6020) ' **** Illegal Laue group : ',nlaue 6020 format(1x,a,i6) call putlin(strout,'CURWIN') call ccperr(1,' STOP in SYMLIB.for 7711') endif C strout = 'Asymmetric unit: '//strout call putlin(strout,'CURWIN') C C Print symmetry operations call prtrsm(pgname, nsymp, rsymiv) C endif C return end C C C subroutine asuput1(ihkl,jhkl,isym) C ================================= C C Put reflection into asymmetric unit defined by call to ASUSET C C On input: C IHKL(3) input indices hkl C C On output: C JHKL(3) output indices hkl C ISYM symmetry number for output C odd numbers are for I+ C even numbers are for I- C real-space symmetry operation number L = (ISYM-1)/2 + 1 C C The real-space symmetry matrices are applied by premultiplying them C by a row vector hkl, ie (h'k'l') = (hkl)R C C C Arguments integer ihkl(3), jhkl(3), isym C C Common block /RECSYM1/ C RSYM real-space symmetry operators C RSYMIV their inverse C NSYM number of symmetry operations C NSYMP number of primitive symmetry operations C NLAUE number of Laue group integer maxsym parameter (maxsym=192) common /recsym1/rsym(4,4,maxsym),rsymiv(4,4,maxsym), . nsym,nsymp,nlaue integer nsym,nsymp,nlaue real rsym,rsymiv save /recsym1/ C C Routines integer inasu external inasu C Locals integer i,l,isgn character*100 linerr C C do 10, l=1,nsymp C h' = h R ie row vector h premultiplies symmetry matrix do 20, i=1,3 jhkl(i) = ihkl(1)*rsym(1,i,l) + ihkl(2)*rsym(2,i,l) . + ihkl(3)*rsym(3,i,l) 20 continue C Test this index against the asymmetric unit C Function INASU returns 0 if outside au, else +-1 depending on sign C isgn = inasu(jhkl, nlaue) if (isgn .ne. 0) go to 100 C Failed, try next symmetry 10 continue C C Shouldn't get here, can't reduce reflection to asymmetric unit write (linerr,'(A,3I4,A)') . 'ASUPUT: can''t put reflection ',ihkl, . ' into asymmetric unit' C Fatal error, stop in routine call lerror(2,-1,linerr) C C Succesful, multiply by sign 100 do 101, i=1,3 jhkl(i) = jhkl(i)*isgn 101 continue C if (isgn .gt. 0) then C I+, odd ISYM isym = l*2 - 1 else C I-, even ISYM isym = l*2 endif C return end C C C subroutine newcel( $ hklmat, hklinv, cell1, umat1, cell2, umat2, ierr) C ====================================================== C C Determine new cell produced by reindexing operation C C On entry: C hklmat(3,3) reindexing matrix HT, h' = HT h C (transpose of matrix H which postmultiplies C row vector hT, h'T = hT H) C hklinv(3,3) inverse of HT (hklmat**-1) C cell1(6) original cell C umat1(3,3) original orientation matrix U C C On exit: C cell2(6) new cell C umat2(3,3) new orientation matrix U' C ierr error flag, = 0 OK, .ne. 0 failure C = +1 left-handed umat1 C = -1 null umat1 C = -2 null cell C C implicit none C real hklmat(3,3), hklinv(3,3), cell1(6), umat1(3,3), $ cell2(6), umat2(3,3) integer ierr C C Locals real lambda, b(3,3), ub1(3,3), ub2(3,3), t1(3,3), t2(3,3), $ det, bm1(3,3), dtorad integer i, iaxes(3) external matmul C C lambda is a dummy wavelength which cancels itself out data lambda /1.0/ C ixes are dummy axis permutations data iaxes/1,2,3/ C C C------------------------------------------------------------ C ierr = 0 dtorad = atan(1.0) / 45.0 C call mkbmat(cell1, b, lambda, iaxes) c c Check matrices are valid call minv3(t1, umat1, det) if (det .lt. 0.0) then write (6, '(/a,g10.5/)') $ ' *** WARNING: Umat is left-handed - ',det ierr = +1 return endif if (abs(det) .lt. 0.9) then ierr = -1 return endif call minv3(t1, b, det) if (abs(det) .lt. 1.0e-20) then ierr = -2 return endif C Matrix UB call matmul(ub1, umat1, b) C C New indices h'T = (H h)T (where H = hklmat, T indicates transpose) C ie h' = HT h C we want (UB)' h' = (UB) h C so (UB)' = (UB) HT**-1 (HT**-1 = hklinv, (UB)' = ub2) C call matmul(ub2, ub1, hklinv) C Decompose new UB matrix into U & B to get new cell call transp(t1, ub2) C Reciprocal metric tensor G**-1 in T2 = (UB)T . (UB) call matmul(t2, t1, ub2 ) C Invert to get real metric tensor in T1 call minv3(t1,t2,det) C Get cell dimensions from real metric tensor C (CLCALC returns a,b,c/lambda, cos alpha, beta, gamma) call clcalc(cell2,t1) C Convert cell do 15, i = 1,3 cell2(i) = cell2(i) * lambda cell2(i+3) = acos(cell2(i+3))/dtorad 15 continue C New B matrix call mkbmat(cell2,b,lambda,iaxes) C BM1 = B**-1 call minv3(bm1,b,det) C U = (UB) . B**-1 call matmul(umat2, ub2, bm1) C return end C C C subroutine clcalc(cell,a) C ========================= C implicit none C C C Calculate cell (real or reciprocal) from metric tensor A C C C C .. Array Arguments .. real a(3,3),cell(6) C .. C .. Local Scalars .. integer i C .. C .. Intrinsic Functions .. intrinsic sqrt C .. C C do 10 i = 1,3 cell(i) = sqrt(a(i,i)) 10 continue C C---- Cosines of angles C cell(4) = a(2,3)/ (cell(2)*cell(3)) cell(5) = a(1,3)/ (cell(1)*cell(3)) cell(6) = a(1,2)/ (cell(1)*cell(2)) C C end C C subroutine mkbmat(cell, bmat, lamda, laxv) C ========================================== C C Makes reciprocal cell orthogonalization matrix BMAT from C direct cell CELL and wavelength LAMDA, using permutation C flags LAXV C C Based on subroutines CONDTR & MATFST C A. Messerschmidt C Max-Planck-Institut fuer Biochemie C D8033 Martinsried C West Germany C C FORTRAN 77: No exceptions. C C Phil Evans April 1988 C C Dummy variables real cell(6), bmat(3,3), lamda integer laxv(3) C C Local variables integer laxs(3) C real parm(9), star(9), winkel(3) C real aa,bb,cc,ca,cb,cg,sa,sb,sg real as,bs,cs,cas,cbs,cgs,sas,sbs,sgs equivalence (star(1), as), (star(2), bs), (star(3), cs), * (star(4),cas), (star(5),cbs), (star(6),cgs), * (star(7),sas), (star(8),sbs), (star(9),sgs) C real sine, cosang, pv, pvstar, dg2rd integer i, lax, lup, lin C C Multiply by this to convert degrees to radians: C parameter (dg2rd = 3.141592654 / 180.) C C **** Inline function definitions **** C sine(cosang) = sqrt( 1. - cosang * cosang) C aa = cell(1) bb = cell(2) cc = cell(3) ca = cos ( cell(4) * dg2rd ) cb = cos ( cell(5) * dg2rd ) cg = cos ( cell(6) * dg2rd ) sa = sin ( cell(4) * dg2rd ) sb = sin ( cell(5) * dg2rd ) sg = sin ( cell(6) * dg2rd ) C C Look at those equivalences above! C pv = sqrt(1. - ca * ca - cb * cb - cg * cg + 2. * ca * cb * cg) C C * AS,BS,CS,CAS,CBS,CGS * RECIPROCAL LATTICE (A*,B*,C*,COS(ALPHA*),...) C as = sa / (aa * pv) bs = sb / (bb * pv) cs = sg / (cc * pv) cas = (cb * cg - ca) / (sb * sg) cbs = (cg * ca - cb) / (sg * sa) cgs = (ca * cb - cg) / (sa * sb) do 200 i = 4, 6 C star(i+3) = sine(star(i)) 200 continue C do 300 i = 1, 3 winkel(i) = acos(star(i+3)) / dg2rd C C * Normalize Ewald-sphere, radius = 1 C star(i) = star(i) * lamda parm(i+3) = star(i) parm(i+6) = winkel(i) * dg2rd C 300 continue C pvstar = sqrt( 1. -cas*cas - cbs*cbs - cgs*cgs+ 2. *cas*cbs*cgs) C lax = iabs(laxv(1)) lup = iabs(laxv(2)) lin = iabs(laxv(3)) do 305 i = 1, 3 laxs(i)=isign(1,laxv(i)) 305 continue C C Form orientation matrix of reciprocal lattice C bmat(2, lax) = 0. bmat(3, lax) = 0. bmat(3, lup) = 0. bmat(1, lax) = star(lax) * laxs(1) bmat(1, lup) = star(lup) * star(3 + lin) * laxs(1) bmat(1, lin) = star(lin) * star(3 + lup) * laxs(1) bmat(2, lup) = star(lup) * star(6 + lin) * laxs(2) bmat(2, lin) =-star(lin) * laxs(2) / star(6 + lin) * * ( star(3 + lup) * star(3 + lin) - star(3 + lax) ) bmat(3, lin) = star(lin) / star(6 + lin) * laxs(3) * * pvstar return C ====== C end C C C subroutine unitmt(u,n) C ====================== C C Set nxn matrix U to unit matrix C integer n real u(n,n) C integer i,j C do 10, j = 1,n do 20, i = 1,n u(i,j) = 0.0 20 continue u(j,j) = 1.0 10 continue return end C C C subroutine modbhd(mflin, mflout, nbatch, jbatch, hklmat, hklinv) C ================================================================ C C Update all batch headers by reindexing transformation C C On entry: C mflin, mflout "streams" for hkl file C nbatch number of batches C jbatch(nbatch) list of batch numbers C hklmat(3,3) reindexing matrix HT C hklinv(3,3) inverse of HT C C implicit none C integer mflin, mflout, nbatch, jbatch(nbatch) real hklmat(3,3), hklinv(3,3) C C+++ orient.h +++ C C C C Orientation block data C C This contains slots for all information that seems to be essential C at present. Each group of parameters is padded at the end for future C expansion. C C Data in the orientation block are referred to the "Cambridge" C laboratory axis frame: x along the (idealized) X-ray beam, z along C usual rotation axis E1 (omega on 3-axis system). The matrix Q converts C a vector in the Madnes frame to the Cambridge frame. Note that the C laboratory frame is essentially defined by the vectors e1,e2,e3 & C source. It doesn't really seem necessary to carry through a whole lot C of crystal and beam tensors, particularly as we have integrated C intensities at this stage, but maybe someone will want to, using the C allocated padding C C The general orientation equation is C C x = R M U B h C C where x position in laboratory frame C R goniostat matrix C M missetting angle matrix (if relevant, see MISFLG) C PhiZ PhiY PhiX (PHIXYZ) C U crystal orientation matrix UMAT C B cell orthogonalization matrix, derived from cell dimensions C h reflection indices C C Note that the description below is NOT is the same order as in the C common block, in which all the integers come before all the reals C (flagged as I or R in the description below) C CI NWORDS number of words in orientation block CI NINTGR number of integers (first part of block, C includes these counts) CI NREALS number of reals CI IORTYP type of orientation block (for possible future use, now = 0) CI INTPAD(8) padding for future use (integers) C C--- Information for this crystal C CR CELL(6) cell dimensions (A & degrees) CI LBCELL(6) refinement flags for cell dimensions CR UMAT(3,3) orientation matrix U. If MISFLG .gt. 0, U is the C "standard" setting when PhiXYZ ==0 CI MISFLG status of "missetting" angles PHIXYZ C = 0 PHIXYZ not used, all orientation in UMAT C = 1 1 set of missetting angles (PHIXYZ(I,1)) C = 2 2 sets PHIXYZ(I,J), J=1,2 CR PHIXYZ(3,2) missetting angles at beginning & end of rotation CI JUMPAX reciprocal axis closest to principle goniostat axis E1 C (only used for printing) CI NCRYST crystal number: a crystal may contain several batches CI LBSETID dataset number (this indexes a list of datasets in the C file header) CI LCRFLG type of crystal mosaicity information C (=0 for isotropic, =1 anisotropic) C *** CRYDAT(12) equivalenced to following *** CR ETAD reflection width (full width) (degrees) (if LCRFLG=0) C or CR ETADH,ETADV horizontal & vertical reflection width (if LCRFLG=1) CR rest of CRYDAT: padding for crystal information (eg more complicated C mosaicity model) C *** C C--- Information for this batch C CI LDTYPE type of data C = 1 oscillation data (2D spots) C = 2 area detector data (3D spots) C = 3 Laue data CR DATUM(3) datum values of goniostat axes, from which Phi is measured C (degrees) CR PHISTT,PHIEND start & stop values of Phi (degrees) relative to datum CR PHIRANGE range of Phi values: typically this will be PHIEND-PHISTT, C but storing this explicitly allows a distinction C eg between a rotation of +160 degrees from a rotation C of -200 degrees CI JSCAXS goniostat scan axis number (=1,2,3, or =0 for C multiple axis scan CR SCANAX(3) rotation axis in laboratory frame (not yet implemented: C only relevant if JSCAXS=0) CR TIME1, TIME2 start & stop times in minutes CI NBSCAL number of batch scales & Bfactors plus SD's C (4 at present, BSCALE, BBFAC & sd's) C set = 0 if batch scales unset CR BSCALE batch scale CR BBFAC batch temperature factor C corresponding scale is exp(-2 B (sin theta/lambda)**2) CR SDBSCL sd (Bscale) CR SDBFAC sd (BBfac) CR BATPAD(11) padding for batch information C C--- Crystal goniostat information C CI NGONAX number of goniostat axes (normally 1 or 3) CI E1(3),E2(3),E3(3) vectors (in "Cambridge" laboratory frame, see below) C defining the NGONAX goniostat axes CC GONLAB(3) names of the three goniostat axes CR GONPAD(12) padding for goniostat information C C--- Beam information C CR SOURCE(3) Idealized (ie excluding tilts) source vector C (antiparallel to beam), in "Cambridge" laboratory frame CR S0(3) Source vector (antiparallel ! to beam), in C "Cambridge" laboratory frame, including tilts CI LBMFLG flag for type of beam information following C = 0 for ALAMBD, DELAMB only (laboratory source) C = 1 ALAMBD,DELAMB,DELCOR,DIVHD,DIVVD (synchrotron) C (other options could include white beam) C *** BEMDAT(25) equivalenced to following *** CR ALAMBD Wavelength in Angstroms CR DELAMB dispersion Deltalambda / lambda. CR DELCOR Correlated component of wavelength dispersion. CR DIVHD Horizontal beam divergence in degrees. CR DIVVD Vertical beam divergence (may be 0.0 for isotropic beam C divergence. CR rest of BEMDAT: padding for beam information C (*** How much here for Laue? ***) C *** C C--- Detector information C CI NDET number of detectors (current maximum 2) C -- for each detector CR DXn crystal to detector distance (mm) CR THETAn detector tilt angle (=Madnes:tau2) (degrees) CR DETLMn(2,2) minimum & maximum values of detector coordinates (pixels) C (i,j): i = 1 minimum, = 2 maximum C j = 1 Xdet, = 2 Ydet C The exact detector frame is not important, but Ydet C should be the axis ~ parallel to the pricipal C rotation axis CR DETPAD(34) padding for detector information C C C Lengths: MBLENG is total length of block, MBLINT, MBLREA are numbers C of integers & reals integer mbleng,mblint,mblrea parameter (mbleng=185,mblint=29,mblrea=156) integer nwords,nintgr,nreals,iortyp,lcell,misflg, . jumpax,ncryst,lcrflg,ldtype,jscaxs,nbscal,ngonax,lbmflg, $ ndet,LBSETID,intpad real cell,umatr,phxyz,crydat,etad,datum, $ phistt,phiend,PHIRANGE,scanax,time1,time2,gonpad, $ bscale,bbfac,sdscal,sdbb,batpad,e1,e2,e3, $ source,s0,bemdat,alambd,delamb,delcor,divhd,divvd, $ dx1,theta1,detlm1,dx2,theta2,detlm2,detpad character btitle*70, gonlab*8 C common /corien/btitle,gonlab(3) C common /orient/ nwords,nintgr,nreals,iortyp,lcell(6),misflg, . jumpax,ncryst,lcrflg,ldtype,jscaxs,nbscal,ngonax,lbmflg, $ ndet,LBSETID,intpad(8), + cell(6),umatr(3,3),phxyz(3,2),crydat(12),datum(3), $ phistt,phiend,scanax(3),time1,time2, $ bscale,bbfac,sdscal,sdbb,PHIRANGE,batpad(11), $ e1(3),e2(3),e3(3), $ gonpad(12),source(3),s0(3),bemdat(25), $ dx1,theta1,detlm1(2,2),dx2,theta2,detlm2(2,2),detpad(33) REAL RDUMMY (MBLENG) EQUIVALENCE (RDUMMY,NWORDS) save /corien/, /orient/ C C *** Note well the following equivalences: these are to allow for C alternative definitions of crystal mosaicity, beam parameters, etc equivalence (etad,crydat(1)) equivalence (bemdat(1),alambd) equivalence (bemdat(2),delamb) equivalence (bemdat(3),delcor) equivalence (bemdat(4),divhd) equivalence (bemdat(5),divvd) C *** C--- orient.h --- C Locals integer k, ibatch, ierr C C--- Loop batches do 10, k = 1, nbatch C C Read next batch header call lrbat(mflin, ibatch, rdummy, btitle, 0) if (ldtype .gt. 0 .and. ldtype .le. 3) then C C convert cell & orientation matrix Umat call newcel( $ hklmat, hklinv, cell, umatr, cell, umatr, ierr) if (ierr .ne. 0) then call ccperr(1, '*** ERROR in new cell in modbhd ***') endif C C Write modified batch header call lwbat(mflout, ibatch, rdummy, btitle) endif 10 continue C return end