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 C VECTORS C C Phil Evans C MRC LMB Cambridge C July 1990 C C Generate all Patterson vectors from a list of input atoms, and C produce a list of all vectors which fall within the volume of the C Patterson calculated (optionally reads Patterson map header to find this out) C suitable for plotting into the Patterson map C C------------------------------------------------------------------------- C C Input: C C 1) Patterson map (MAPIN): the header from this file is used to define C the volume of space into which the vectors will be put. The program C assumes that the symmetry stored in the map header is the Patterson C symmetry for the true Patterson space-group. This is used to move C vectors into the volume of the map, but may be overridden by the C PSYMMETRY command C If the command XYZLIM is present, a map will not be read, and C the volume, cell & symmetry will be taken from the command input: C in this case, XYZLIM, CELL & PSYMMETRY commands must be present C C 2) Control input (DATA) C C TITLE title for the run, written to the output vector file C C SYMMETRY space-group name OR number OR symmetry operation C Space-group symmetry for the atoms (NOT the Patterson) C This may be given in 3 ways:- C (a) space-group name C (b) space-group number C (c) symmetry operations as in International Tables, separated C by '*', on a series of SYMMETRY lines if necessary C For options (a)& (b), symmetry operations are read from the C library file SYMOP C C PSYMMETRY space-group name OR number OR symmetry operation C Space-group symmetry for the Patterson: by default this is taken C from the Patterson header, but may overridden by this command C This may be given in 3 ways, as for SYMMETRY C PSYMMETRY P1 will give the unique set of vectors, provided they C lie within the cell volume given C C NCODE orthogonalization code C Vectors are written out to a PDB file in an orthogonal frame C defined by this code. The default is NCODE = 1, which is the C usual frame for Brookhaven files, so normally this command C can be omitted C C XZYLIM xmin xmax ymin ymax zmin zmax C Limits of Patterson volume in fractions of the unit cell. C If this command is present, a Patterson map file will not C be read, and CELL, PSYMMETRY commands must be present; a C GRID command is optional C C CELL a b c alpha beta gamma C Unit cell, by default taken from map header. This command C is compulsory if no map is given (ie if XYZLIM is present) C C GRID nx ny nz C Map grid used to convert vectors to grid coordinates. By default, C taken from map header. Optional if no map is given (else set to C 100 100 100) C C ATOM atomname x y z C Input an atom to be used to predict vectors. "atomname" is a C 1-character unique identifier for this site, x,y,z are the C FRACTIONAL coordinates C C END C End of input, also end-of-file will do C C------------------------------------------------------------------------- C C Output: C C The vectors are written to a standard PDB file (XYZOUT) as orthogonal C coordinates with atom name "VAB" where A & B are the atomnames of the 2 atoms C involved. This file may be plotted on the Patterson using Pluto, with C only the identity symmetry operation defined there. Duplicate vectors are C removed, but all equivalent vectors in the map volume should be present C C------------------------------------------------------------------------- C C VAX example: C C $ vectors mapin x16patt6a xyzout x16vecs.pdb C $ go C SYMMETRY 214 C ATOM B 0.0446 0.0344 -0.0317 C ATOM C -0.0051 0.0508 -0.0194 C END C $! C $! Contour Patterson with vectors C $! C $plot: pluto mapin X16patt6a xyzin X16vecs.pdb plot X16patt6a C $ go C TITLE Vectors for sites B & C C SYMM X,Y,Z C MAP SCALE 2 C CONTRS 10 20 30 40 50 60 80 100 200 400 600 800 C MODE BELOW 15 DASHED C SECTNS 0 18 C GRID 14.81 14.81 C INPUT BROOK C SOLID C RADII ATOMS VB -0.3 VC -0.5 C PLOT C $! C C------------------------------------------------------------------------- C PROGRAM VECTORS C =============== C C Map things CHARACTER*80 TITLE INTEGER NSEC,IUVW(3),MXYZ(3),NW1,NU1,NU2,NV1,NV2, . LSPGRP,LMODE REAL CELL(6),RHMIN,RHMAX,RHMEAN,RHRMS REAL FLMMAP(2,3),FBOX(3) INTEGER NCT(3) LOGICAL MAP C C Atoms INTEGER MAXATM PARAMETER (MAXATM=200) CHARACTER*1 ATNAM(MAXATM) REAL X(3,MAXATM) INTEGER NATOMS, ISYM(MAXATM) C C Vectors INTEGER MAXVEC PARAMETER (MAXVEC=10000) CHARACTER*4 VECNAM(MAXVEC) REAL V(3,MAXVEC) INTEGER NVEC, IJSYM(2,MAXVEC) CHARACTER*72 VTITLE C C Orthogonalization INTEGER NCODE REAL RO(3,3), ROR(3,3,6), VOLUME C C Symmetry C Symmetry from map file INTEGER NPASYM,LPASGP REAL PATROT(4,4,192) CHARACTER*10 NAMPSG,PGNAMP C Space-group symmetry for atoms INTEGER NATSYM,LATSGP REAL ATMROT(4,4,192) CHARACTER*10 NAMSGP,PGNAME C C Things for PARSER ---- INTEGER MAXTOK PARAMETER (MAXTOK=100) CHARACTER KEY*4,CVALUE(MAXTOK)*4,LINE*80 INTEGER IBEG(MAXTOK),IEND(MAXTOK),ITYP(MAXTOK),IDEC(MAXTOK),NTOK REAL FVALUE(MAXTOK) LOGICAL LEND C C Main input keywords INTEGER NKEYS PARAMETER (NKEYS=9) CHARACTER*9 KEYS(NKEYS) C C General INTEGER I,IZ,IZZ,J,K,L,IX,IY,IFLAG,IRES,NRESET,NREJCT REAL Y(3),Z(3),U(3),VTEST(3),DTEST,FR,FD,Q,BISO,B(6),VEC_LEN LOGICAL CELMUL,PRINT,LPASYM CHARACTER*4 VNAM,RESNAM,ID,CHNNAM*1,RESNO,INSCOD*1,ALTCOD*1, . SEGID C C Externals EXTERNAL CCPERR,CCPFYP,CCPMVI,CCPRCS,CCPUPC,CHKKEY,GTPINT, . GTPREA,MATVC4,MATVEC,MRDHDR,MSYMOP,PARSER,RBFRO1, . RDSYMM,VDIF,VECTST,XYZADVANCE,XYZATOM, . XYZCLOSE,XYZCOORD,XYZINIT,XYZOPEN C C Data Statements DATA KEYS/'TITLE','SYMMETRY','PSYMMETRY','NCODE','ATOM', . 'XYZLIM','CELL','GRID','END'/ DATA VTITLE/'Patterson vectors'/ C C******************************************************************** C CALL CCPFYP CALL CCPRCS(6,'VECTORS','$Date: 2002/08/12 08:55:18 $') CALL XYZINIT KDUMMY = 80 KFAIL = 0 C C Initializations NATOMS = 0 MAP = .TRUE. NATSYM = 0 NPASYM = 0 LPASYM = .TRUE. SymHist = 0 NCODE = 1 PRINT = .TRUE. MXYZ(1) = 100 MXYZ(2) = 100 MXYZ(3) = 100 CELL(1) = -100. LPASGP = 0 ALTCOD = ' ' INSCOD = ' ' SEGID = ' ' RESNO = ' ' IZZ = 0 BISO = 1.0 C DTEST is limit for vector coincidence in A DTEST=0.2 C 10 NTOK=MAXTOK LINE=' ' CALL PARSER( . KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,LEND,.TRUE.) C IF(LEND) GO TO 1000 C Ignore blank line IF(NTOK.EQ.0) GO TO 10 C C Make uppercase CALL CCPUPC(KEY) CALL CHKKEY(KEY,KEYS,NKEYS,K) IF(K.EQ.0) THEN WRITE(6,FMT='(A)') ' ** Unrecognized keyword **' GO TO 10 ELSE IF(K.LT.0) THEN WRITE(6,FMT='(A)') ' ** Ambiguous keyword **' GO TO 10 ENDIF C C 19 GO TO . (20,30,35,40,50,60,70,80,1000),K C--------------------------------------------------------------------- C C TITLE 20 VTITLE = LINE(IBEG(2):) GO TO 10 C--------------------------------------------------------------------- C C SYMMETRY read atom symmetry 30 CONTINUE J=2 CALL RDSYMM(J,LINE,IBEG,IEND,ITYP,FVALUE,NTOK, . NAMSGP,LATSGP,PGNAME,NATSYM,NSYMP,ATMROT) GO TO 10 C C--------------------------------------------------------------------- C C PSYMMETRY read Patterson symmetry 35 CONTINUE IF (LPASYM) THEN C If LPASYM .true., then this is the first PSYMMETRY command, so clear C symmetry count NPASYM = 0 LPASYM = .FALSE. ENDIF J=2 CALL RDSYMM(J,LINE,IBEG,IEND,ITYP,FVALUE,NTOK, . NAMPSG,LPASGP,PGNAMP,NPASYM,NPASYP,PATROT) GO TO 10 C C--------------------------------------------------------------------- C C NCODE read orthogonalization code C 40 CALL GTPINT(2,NCODE,NTOK,ITYP,FVALUE) GO TO 10 C--------------------------------------------------------------------- C C ATOM read atom name (1-character) x,y,z C 50 CONTINUE IF (NTOK .LT. 5) THEN WRITE(6,*) ' *** Too few fields for ATOM ***' GO TO 10 ENDIF C IF (NATOMS .GE. MAXATM) THEN WRITE(6,*) ' *** More than ',MAXATM,' atoms given ***' GO TO 1000 ENDIF C NATOMS = NATOMS + 1 ATNAM(NATOMS) = CVALUE(2)(1:1) ISYM(NATOMS) = 1 DO 51, I=1,3 CALL GTPREA(I+2,X(I,NATOMS),NTOK,ITYP,FVALUE) 51 CONTINUE GO TO 10 C--------------------------------------------------------------------- C XYZLIM set fractional limits for Patterson volume, as alternative C to reading Patterson header 60 K=1 DO 61, J=1,3 DO 62, I=1,2 K=K+1 CALL GTPREA(K,FLMMAP(I,J),NTOK,ITYP,FVALUE) 62 CONTINUE 61 CONTINUE MAP=.FALSE. GO TO 10 C C--------------------------------------------------------------------- C CELL read unit cell 70 DO 71, I=1,6 CALL GTPREA(I+1,CELL(I),NTOK,ITYP,FVALUE) 71 CONTINUE GO TO 10 C C--------------------------------------------------------------------- C GRID read sampling grid 80 DO 81, I=1,3 CALL GTPINT(I+1,MXYZ(I),NTOK,ITYP,FVALUE) 81 CONTINUE GO TO 10 C C--------------------------------------------------------------------- C C GO, calculate vectors C 1000 WRITE(6,'(1X,I8,A,/)') NATOMS,' atoms read' IF (MAP) THEN WRITE(6,*) 'Reading Patterson map header' C C Open Patterson map file & read header C CALL MRDHDR(3,'MAPIN',TITLE,NSEC,IUVW,MXYZ,NW1,NU1,NU2,NV1,NV2, . CELL,LSPGRP,LMODE,RHMIN,RHMAX,RHMEAN,RHRMS) C C Read symmetry of Patterson IF (LPASYM) THEN CALL MSYMOP(3,NPASYM,PATROT) LPASYM = .FALSE. ENDIF C C Store limits of Patterson map in fractional coordinates FLMMAP(1,IUVW(1)) = NU1 FLMMAP(2,IUVW(1)) = NU2 FLMMAP(1,IUVW(2)) = NV1 FLMMAP(2,IUVW(2)) = NV2 FLMMAP(1,IUVW(3)) = NW1 FLMMAP(2,IUVW(3)) = NW1+NSEC-1 ENDIF C C Check that CELL & Patterson symmetry have been given IF (CELL(1) .LT. 0.) THEN WRITE(6,*) '****** No unit cell given ******' CALL CCPERR(1,' Program Terminated ') ENDIF IF (LPASYM) THEN WRITE(6,*) '****** No Patterson symmetry given ******' CALL CCPERR(1,' Program Terminated ') ENDIF C C CELMUL is .true. if the map volume is more than 1 unit cell in any direction CELMUL = .FALSE. DO 1001, J=1,3 DO 1002, I=1,2 C Convert to fractional if Map read, else already fractional IF (MAP) FLMMAP(I,J) = FLMMAP(I,J)/REAL(MXYZ(J)) 1002 CONTINUE C FBOX is length of box in each direction FBOX(J) = FLMMAP(2,J) - FLMMAP(1,J) C NCT is the number of unit cells in each direction, C = 0 if smaller than cell NCT(J) = INT(FBOX(J)) IF (NCT(J) .GE. 1) CELMUL = .TRUE. 1001 CONTINUE C DO 1012, I=1,3 1012 VTEST(I) = DTEST/CELL(I) C C WRITE(6,6090) . 'Space group for atoms : ',NAMSGP,', ',NATSYM, . ' symmetry operations' IF (LPASGP .EQ. 0) LPASGP = LSPGRP WRITE(6,6010) . 'Space group no. for Patterson: ',LPASGP,', ',NPASYM, . ' symmetry operations' 6090 FORMAT(/,1X,A,A,A,I4,A,/) 6010 FORMAT(/,1X,A,I6,A,I4,A,/) IF (.NOT.MAP) WRITE(6,6091) FLMMAP 6091 FORMAT(' Limits of Patterson volume (fractional): ',3(2X,2F8.3)) C C Set up matrices needed C C Orthogonalization C ROR contains all fractional to orthogonal matrices C Put the chosen one (NCODE) in RO CALL RBFRO1(CELL,VOLUME,ROR) CALL CCPMVI(RO,ROR(1,1,NCODE),9) WRITE (6,6100) NCODE, ((RO(I,J),J=1,3),I=1,3) 6100 FORMAT(/, . ' Fractional to Orthogonal matrix for orthogonalization code ', . I3,//,3(11X,3F10.5,/)) C C C Generate all symmetry related atoms K = NATOMS IF (PRINT) THEN WRITE(6,6040) 6040 FORMAT( . ' All symmetry related atoms:',//, . ' Atom Sym Satom x y z',/) ENDIF IF (PRINT) THEN L = 1 DO 1020, J=1,NATOMS WRITE(6,6050) J,L,J,(X(I,J),I=1,3) 6050 FORMAT(1X,3I5,4X,3F10.5) 1020 CONTINUE ENDIF C DO 1030, L=2,NATSYM DO 1030, J=1,NATOMS IF (K .GE. MAXATM) THEN WRITE(6,*) . ' *** More than ',MAXATM,' symmetry-generated atoms ***' GO TO 1100 ELSE K = K+1 CALL MATVC4(X(1,K),ATMROT(1,1,L),X(1,J)) ATNAM(K) = ATNAM(J) ISYM(K) = L IF (PRINT) THEN WRITE(6,6050) J,L,K,(X(I,K),I=1,3) ENDIF ENDIF 1030 CONTINUE C 1100 NATOMS = K C C Generate all vectors, testing them for coincidence C NVEC = 0 NRESET = 0 NREJCT = 0 C DO 1200, J=1, NATOMS-1 DO 1200, I=J+1,NATOMS C Y is basic vector CALL VDIF(Y, X(1,J), X(1,I)) C Make vector name IF (LLE(ATNAM(I),ATNAM(J))) THEN VNAM = 'V'//ATNAM(I)//ATNAM(J) ELSE VNAM = 'V'//ATNAM(J)//ATNAM(I) ENDIF C Now try all Patterson symmetry operations to try to put it back into the C volume of the map, including translations DO 1210, L=1,NPASYM CALL MATVC4(Z, PATROT(1,1,L), Y) C See if this vector can be put back into the map volume by cell translation DO 1220, K=1,3 FR = Z(K) - FLMMAP(1,K) FD = AMOD(FR, 1.0) IF (FD.LT.0.0) FD = FD + 1.0 C FD is coordinate relative to start of box C Test against box size, jump out as soon as it fails IF (FD .GT. FBOX(K)+VTEST(K)) GO TO 1210 C If vector is just outside box, put it back on edge IF (FD .GT. FBOX(K)) THEN FD = FBOX(K) NRESET = NRESET+1 ENDIF C Required translation is FD - FR Z(K) = Z(K) + FD - FR 1220 CONTINUE C We have a vector Z in the map box IF (CELMUL) THEN C At least one direction has more than 1 unit cell, so we must loop C multiple cell translations DO 1230, IX=0,NCT(1) DO 1230, IY=0,NCT(2) DO 1230, IZ=0,NCT(3) U(1) = Z(1) + REAL(IX) U(2) = Z(2) + REAL(IY) U(3) = Z(3) + REAL(IZ) C Test against top of box DO 1240, K=1,3 IF (U(K) .GT. FLMMAP(2,K)) GO TO 1230 1240 CONTINUE C Test this vector against all previous vectors, add to list if OK CALL VECTST(U,V,NVEC,VNAM,VECNAM,MAXVEC,VTEST,IFLAG) C IFLAG = -1 if list full IF (IFLAG .LT. 0) GO TO 2000 1230 CONTINUE C ELSE C Map volume is all .lt. 1 cell, no need to loop translations C Test this vector against all previous vectors, add to list if OK CALL VECTST(Z,V,NVEC,VNAM,VECNAM,MAXVEC,VTEST,IFLAG) C IFLAG = -1 if list full IF (IFLAG .LT. 0) GO TO 2000 ENDIF IF (IFLAG .EQ. 0) THEN C A vector has been added to the list, store symmetry numbers of atoms IJSYM(1,NVEC) = ISYM(I) IJSYM(2,NVEC) = ISYM(J) ELSEIF (IFLAG .GT. 0) THEN NREJCT = NREJCT+1 ENDIF C 1210 CONTINUE 1200 CONTINUE C 2000 CONTINUE C C****************************************************************** C C Now write out all vectors to Brookhaven format file (orthogonal) C and (optionally) to printer C KFAIL = 0 IXYZOUT = 0 CALL XYZOPEN('XYZOUT','OUTPUT','PDB',IXYZOUT,KFAIL) C C WHAT SHALL WE DO WHEN CCIF IMPLIMENTED ? C LINE = 'REMARK '//VTITLE CALL WREMARK(IXYZOUT,LINE) CALL WBCELL(IXYZOUT,CELL,NCODE) WRITE(6,6205) NVEC, NREJCT, NRESET 6205 FORMAT(/,' Number of accepted vectors : ',I6,/ . ' Number of rejected duplicates : ',I6,/ . ' Number of vectors reset to edge of box : ',I6,/) IF (PRINT) WRITE(6,6210) 6210 FORMAT(//,' List of vectors',//, . ' Number Fractional coordinate ', . ' Orthogonal coordinate Grid coordinate ', . ' Vector Symmetry Vector_length',/) C Q=1.0 RESNAM = 'VEC' CHNNAM=' ' IRES = 1 ID = ' ' B(1) = BISO DO 2010 I=2,6 B(I) = 0.0 2010 CONTINUE C DO 2100, L=1,NVEC C Orthogonalize CALL MATVEC(Z,RO,V(1,L)) C CALL XYZATOM(IXYZOUT,L,VECNAM(L),RESNAM,CHNNAM,IRES, + RESNO,INSCOD,ALTCOD,SEGID,IZZ,ID) CALL XYZCOORD(IXYZOUT,'O','U',Z(1),Z(2),Z(3),Q,BISO,B) CALL XYZADVANCE(IXYZOUT,0,0,*2100,*2100) C IF (PRINT) THEN C Calculate grid units DO 2110, I=1,3 2110 NCT(I) = NINT(V(I,L) * MXYZ(I)) C Calculate vector lengths vec_len = z(1)*z(1) +z(2)*z(2)+z(3)*z(3) vec_len = sqrt(vec_len) WRITE (6,6300) L,(V(I,L), I=1,3), Z, NCT, VECNAM(L), + (IJSYM(I, L), I=1,2),vec_len 6300 FORMAT(1X,I6,4X,3F10.5,5X,3F8.2,5X,3I6,5X,A,2X,2I3,2X,F10.2) ENDIF 2100 CONTINUE C CALL XYZCLOSE(IXYZOUT) WRITE(6,*) ' ' CALL CCPERR(0,'Normal termination') END C C SUBROUTINE VECTST(Z,V,NVEC,VNAM,VECNAM,MAXVEC,VTEST,IFLAG) C ========================================================== C C Test new vector Z against all vectors so far, add to list if C not coincident with previous ones C C Z(3) vector to test C V(3,MAXVEC) list of vectors C NVEC number of vectors, incremented on output if this vector C is accepted C VNAM vector name, coincident vectors are rejected only if C they have different names C VECNAM(MAXVEC) list of vector names C MAXVEC maximum number of vectors C VTEST(3) limits for coincidence C IFLAG on return = 0 if OK (ie vector added to list) C = 1 if vector already present in list C =-1 if vector list full C INTEGER NVEC,MAXVEC,IFLAG REAL Z(3),V(3,MAXVEC),VTEST(3) CHARACTER*(*) VNAM,VECNAM(MAXVEC) C INTEGER I,J C INTRINSIC ABS C IF (NVEC .GT. 0) THEN DO 10, J=1,NVEC IF (VNAM .EQ. VECNAM(J)) THEN DO 20, I=1,3 IF (ABS(Z(I)-V(I,J)) .GT. VTEST(I)) GO TO 10 20 CONTINUE C We have 2 coincident vectors with the same name, so reject this one IFLAG = 1 RETURN C ====== ENDIF C 10 CONTINUE C ENDIF C This vector is accepted, so add it to the list IF (NVEC .GE. MAXVEC) THEN WRITE(6,*) ' *** Vector list full, maximum ',MAXVEC IFLAG = -1 RETURN C ====== ENDIF C NVEC = NVEC + 1 CALL VSET(V(1,NVEC), Z) VECNAM(NVEC) = VNAM IFLAG = 0 RETURN END C C