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 PDBSET C C--------------------------------------------------------------------------- C C Yet another program to do various useful manipulations on a PDB file C Phil Evans C MRC LMB, Cambridge, September 1992 C C Functions & options C In the description below, optional items are in [], alternatives are C separated by |, keywords are in uppercase, parameters (ie numbers) C are in lowercase. The input itself is case-insensitive for keywords C (but parameters eg chain IDs must of course be the correct case). C In the output file,the chain ID is always uppercase. C C 1) Divide residue ID into chain ID + residue number (if it begins C with a non-digit) (for output from O). This is ALWAYS done, so the C output file always has a valid numerical residue number C C 2) CELL a b c [alpha beta gamma] C Read cell dimensions & make CRYST1 & SCALE header records. These C will replace any CRYST1 & SCALE lines already present in file. C The CRYST1 line should have the spacegroup in it, so a SPACEGROUP C command is recommended. C C 3) ORTHOGONALIZATION (or NCODE) orthogonalization_code C Define code to generate orthogonalization matrix from input C cell. This is not normally required, and only has an effect if C a CELL command is also given C Code :- C = 1 axes along a, c* x a, c* (Brookhaven standard, default) C = 2 axes along b, a* x b, a* C = 3 axes along c, b* x c, b* C = 4 axes along a+b, c* x (a+b), c* C = 5 axes along a*, c x a*, c ( Rollett ) C = 6 axes along a, b*, a x b* C = 7 axes along a*, b, a* x b (TNT convention, C probably not very useful here C since TNT has its own converter C program) C C 4) SPACEGROUP spacegroup_name C read spacegroup name (not essential, but put into CRYST1 line on output) C C 5) SYMGEN Spacegroup_name|Spacegroup_number|Symmetry_operation|NCS C Generate chains with these symmetry operations applied. If the C operations are given explicitly, several SYMMETRY commands may be C given. The identity operation must be specified explicitly if required. C Use the CHAIN command to rename them. Note that, except for NCS, these C symmetry operations apply to fractional coordinates, so the C orthogonalization operation must be know to the program, either from C CRYST1 and/or CELL lines in the input coordinate file, or from a CELL C command. If the keyword NCS is given, then a series of TRANSFORM C commands should be given to define the non-crystallographic symmetry C operations to be used C C 6) RENUMBER [INCREMENT] start|increment [residue range] C [CHAIN old_chain [TO new_chain]] C renumber or add constant to residue numbers in given range. The residue C range is given as 1st_residue_number [TO] last_residue_number. C If the CHAIN keyword is present, the renumbering applies only to C this chain. The option "TO new_chain" clause causes the chain identifier C to be changed. Note that renumbering is done after chain renaming C specified by the CHAIN command, so the chain specified here ("old_chain") C is the chain ID after any renaming. NB there is NO check that different C RENUMBER commands are mutually exclusive. To avoid problems with C recursive renumbering, if more than one RENUMBER commands would apply C to a residue, only the first will be done. C (defaults all residues, all chains) C eg RENUMBER 35 ! renumber all residues, starting from 35 C RENUMBER INCREMENT -5 102 TO 110 CHAIN C ! subtract 5 from C ! residues 102 to 110 in chain C C RENUMBER 101 1 TO 78 CHAIN A TO B C ! renumber residues 1 to 78 in chain A from 101 (to 178), C ! changing the chain identifier to B C C 7) CHAIN [SYMMETRY Nsym] [old_chain] new_chain C change chain ID to given value. C If only one value given, change all chains to this value C If SYMMETRY keyword given, this applies to this symmetry operation only. C A series of CHAIN commands may be given C eg CHAIN Q ! change all chains to Q C CHAIN SYMMETRY 2 A B ! change chain generated from chain A C ! by symmetry operation 2 to B C C 8) BFACTOR [subkey] B_factor C set B-factor (default 20.0) C Subkeys: C ALWAYS (default) reset all B-factors C DEFAULT reset B-factors if = 0.0 C MINIMUM reset if less than given values C MAXIMUM reset if greater than given value C RANGE truncate to given range C C 9) OCCUPANCY [subkey] Occupancy C set occupancy (default 1.0) C Subkeys: C ALWAYS (default) reset all occupancies C DEFAULT reset occupancies if = 0.0 C MINIMUM reset if less than given values C C 10) SELECT [subkeys] C Subkeys: C CHAIN select only specified chain(s) C eg SELECT CHAIN C ! select only chain C C OCCUPANCY [] C select only atoms with occupancy .gt. minimum_occupancy C [ default = 0.0]. This can be used to strip out C dummy atoms with zero occupancy C BFACTOR [] C select only atoms with Bfactor less than C [default = 99.0] C C 11) ROTATE [INVERT] [MATRIX|EULER|POLAR] values C Define rotational transformation, either as MATRIX (this keyword may C be omitted) followed by 9 numbers (r11 r12 r13 r21 r22 r23 r31 r32 r33), C by keyword EULER followed by Eulerian angles alpha, beta gamma (as in C ALMN), or by keyword POLAR followed by polar angles omega, phi, kappa C (as in POLARRFN). This transformation will be applied to all atoms. C The SHIFT command may be used to define a translation in addition. C The transformation defined by ROTATE & SHIFT, or by TRANSFORM, is C applied after any SYMMETRY operation. Multiple definitions of ROTATE C or TRANSFORM, or of SHIFT will NOT be concatenated: only the last will C be effective. C The subkey INVERT causes the inverse transformation to be applied C Note that an INVERT instruction if present will apply to both C ROTATE & SHIFT C C 12) SHIFT [INVERT] [FRACTIONAL] tx ty tz C Define translation transformation (added AFTER rotation). If the C keyword FRACTIONAL is present, the translation is assumed to be C in fractional coordinates, otherwise orthogonal A. C The subkey INVERT causes the inverse transformation to be applied C Note that an INVERT instruction if present will apply to both C ROTATE & SHIFT C C 13) TRANSFORM [INVERT] [FRACTIONAL] C r11 r12 r13 r21 r22 r23 r31 r32 r33 tx ty tz C TRANSFORM [INVERT] ODB [O_database_filename] C Define transformation, equivalent to ROTATE MATRIX + SHIFT. If the C keyword FRACTIONAL is present, the translation is assumed to be C in fractional coordinates, otherwise orthogonal A. C The subkey ODB causes the transformation to be read from a file in the C format of an O datablock transformation. C The subkey INVERT causes the inverse transformation to be applied. C If a SYMGEN NCS command is given before TRANSFORM commands, these C are collected together to generate multiple NCS-symmetry related chains C C 14) REMARK anything C Just gets echoed to output coordinate file C C 15) XPLOR [subkeys] C The input file is assumed to come from Xplor: the following C operations are then done:- C (a) all hydrogens are removed, unless subkeyword HYDROGEN is present C (b) dummy atoms (X .gt. 9000) are removed C (c) the segment identifier (columns 73-76) is used as the CHAIN name C for any chain renaming (etc) commands: thus in this case references C to chains in other commands may have up to 4 characters and are C case-sensitive. Unless renamed, the first character of the segment C identifier is put in the chain ID and made uppercase. C (d) the residue number is read correctly for numbers .ge. 1000 C C 16) PICK atom1 atom2 . . . C Define atom names to be included: all other atoms will be omitted C eg PICK CA to choose C-alpha only C Note that the atomname is case-sensitive C C 17) SEQUENCE [PDB|SINGLE] [sequence file name] C Write out sequence to a file (default file name SEQUENCE). This can C be edited to give a sequence for Xplor or O, etc. If the keyword "PDB" C is present, the sequence is written in PDB SEQRES format, split by chains C If SINGLE is given, the sequence is written in single-letter code C C 18) OUTPUT [switches] C Set output options C a) XPLOR duplicate the chain ID as an Xplor segid, to make the file C suitable for direct input into Xplor C C 19) UTOB C Convert Us on input file to B (B = 8 pi**2 u**2) C C 20) ELEMENT . . . C Define list of 2-character element names to be left-justified C in atomnames, eg MG, FE, Zn. Note that the element name is case-sensitive C The PDB convention defines the first 2 characters of the atomname C as the element name, but Xplor & O put them in the wrong place C "CA" is NOT accepted, as this conflicts with Calpha: you will have C decide what to do with these yourself C C 21) REORTHOGONALIZE [[FROM] ] [TO] C Change orthogonalization convention for coordinates by converting C to fractional in the input convention (FROM) and reorthogonalizing C in the output convention (TO). If the FROM Ncode is omitted, the C orthogonalization will be taken from the input (PDB) file C as SCALEn lines, or the default of Ncode = 1 will be used. If the C cell is not present in the input file, a CELL command must be given C here. is compulsory. See above for Ncodes. C C 22) REPLACE RESIDUE BY C Globally replace residue type, eg REPLACE RESIDUE CYS BY CYH C Useful for renaming according to dictionary conventions of C different programs. The residue names will be right-justified C before use to allow for single character names C eg replace residue C by CYT C REPLACE ATOM BY [IN ] C Replace atom name by new one, optionally only in specified residue C name. Note that replace tests are done in the order given, so an C "IN " command must allow for previous REPLACE RESIDUE C commands . Note also that leading spaces must be given in atom names C eg REPLACE ATOM " O" BY " OW" IN HOH C C 23) EXCLUDE [subkeys] C Exclude some things, depending on subkey: C C HYDROGENS exclude hydrogen atoms (as for the Xplor option) C HEADERS exclude all lines except ATOM & HETATM lines C The default is to copy them from the input file C C 24) COM C Will calculate the centre of mass and maximum distance from it C of the coordinates output. This may be useful for determining the C rotation function integration radius. (Not done by default since C it requires an intermediate file.) C program pdbset C ============== C C implicit none C C things for parser ---- integer maxtok parameter (maxtok=100) character key*4,cvalue(maxtok)*4,line*256,atline*80 integer ibeg(maxtok),iend(maxtok),ityp(maxtok), . idec(maxtok),ntok real fvalue(maxtok) logical lend C real cell(6), ro(4,4), rf(4,4), $ trnsfm(4,4), amf(4,4), rmat(3,3) C Selection integer maxchs parameter (maxchs=20) integer nselch character selchn(maxchs)*4 real qmin, bfmax integer nqmin, nbfmax logical lqmin, lbfmax C C C Pick integer maxpat parameter (maxpat=20) integer nselat character*4 selatm(maxpat) C C C Residue renumbering integer maxren parameter (maxren=50) integer numren, incren(maxren), nrsren(2,maxren) character chnren(maxren)*4, allchn*4, chnrnw(maxren)*4 logical lincr(maxren) C C Symmetry integer maxsym parameter (maxsym=192) integer nsym, nsymp, numsgp, lsymgn real rsym(4,4,maxsym), rfsym(4,4,maxsym) character spgnam*10, pgname*10 C C Chains integer maxchn parameter (maxchn=50) integer nchain, mchain character oldchn(maxchn)*4, newchn(maxchn,maxsym)*4, $ chnold*4, chnnew*4, chains(maxchn)*4, segid*4 c c Element names integer maxelm parameter (maxelm=50) character*2 elmnts(maxelm) integer nelmnt C C Sequence integer maxseq parameter (maxseq=5000) character seqnce(maxseq)*4, seqfil*80, seqnc1(maxseq)*1 logical lseqnc integer nseqn, nsqchn(2,maxchn), kseqpd c c Reorthogonalization integer ncdin, ncdout real roin(4,4), rfin(4,4), roout(4,4), rfout(4,4) c c Replace integer maxrtp parameter (maxrtp=200) integer nreplc, lflrtp(maxrtp) character*4 oldrtp(maxrtp), newrtp(maxrtp), resatp(maxrtp), $ tname C integer iun,iunout,i,j,k,m,n,iser,ires, jres, kres, ierr, lflag, $ jtok, jchwld, isym, ncode, natom, lunit, istat, nbi, iunscr real x(3),xyzmin(3),xyzmax(3),angles(3),xx(3),xcom(3) real q,b,bfac,occ, det, xdummy, utob, bfacmx, maxdist character atnam*4,restyp*4,idch*1,resno*4, filnam*80, insert*1 character insold*1 character*15 codes(7) C integer jbfac, jqfac logical lcell, newres, ltrnsf, lfrtrn, lorthg, lspnam, lxplor, $ lnodum, lhydro, loutxp, lnewch, linvrt, lutob, $ lngnmi, lngnmo, lreort, lnewhd, com, lnohed c logical eqlstr integer lenstr external eqlstr, lenstr C C----------------------------------- C Data statements data codes/'a,c*xa,c*','b,a*xb,a*','c,b*xc,b*', . 'a+b,c*x(a+b),c*','a*,cxa*,c','a,b*,axb*','a*,b,a*xb'/ C marker for all chains data allchn /'*'/ C default transformation data trnsfm $ /1.0,0.0,0.0,0.0, $ 0.0,1.0,0.0,0.0, $ 0.0,0.0,1.0,0.0, $ 0.0,0.0,0.0,1.0/ C Unit number for reading O files data lunit/11/, iunscr /20/ data com /.false./, xcom /3*0.0/ C C----------------------------------- C C========================================================================== C Initializations C========================================================================== call ccpfyp CALL CCPRCS(6,'PDBSET','$Date: 1997/03/11 11:14:38 $') iun=1 iunout = 2 call ccpdpn(iun,'XYZIN','READONLY','F',0,0) call ccpdpn(iunout,'XYZOUT','NEW','F',0,0) C n=0 lseqnc = .false. kseqpd = 0 seqfil = 'SEQUENCE' nseqn = 0 nsqchn(1,1) = 1 loutxp = .false. lcell = .false. lnewhd = .false. lspnam = .false. lxplor = .false. lnodum = .false. lhydro = .false. lreort = .false. lnohed = .false. nreplc = 0 ncdin = -1 ncdout = -1 nchain = 0 mchain = 0 ltrnsf = .false. lfrtrn = .false. linvrt = .false. lutob = .false. utob = ((4.0*atan(1.0))**2)*8.0 jbfac = 0 jqfac = 0 qmin = 0.0 lqmin = .false. bfmax = 99.0 lbfmax = .false. do 1, i=1,maxren lincr(i) = .false. 1 continue do 2, i = 1,3 xyzmin(i) = +10000000. xyzmax(i) = -xyzmin(i) 2 continue xdummy = 1499. nselch = 0 nselat = 0 numren = 0 nelmnt = 0 spgnam = ' ' nsym = 0 lsymgn = 0 bfac = 20.0 bfacmx = 99.0 occ = 1.0 ncode = 1 ierr = 0 insold = ' ' C C========================================================================== C Read command input C========================================================================== 100 line=' ' ntok=maxtok call parser( . key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.) if(lend) go to 200 if(ntok.eq.0) go to 100 C C........................................................................ if (key .eq. 'END') then go to 200 C........................................................................ elseif (key .eq. 'REMA') then if (ntok .gt. 1) then k = min(iend(ntok), ibeg(2)+71) write (iunout, '(A,A)') 'REMARK ',line(ibeg(2):k) else write (iunout, '(A)') 'REMARK ' endif C........................................................................ elseif (key .eq. 'CELL') then C CELL read unit cell, & set flag to make headers if (ntok .lt. 4) then write (6,'(a)') ' ** Too few numbers on CELL command **' ierr = ierr+1 go to 100 endif cell(4) = 90.0 cell(5) = 90.0 cell(6) = 90.0 do 110, i=1,6 call gttrea(i+1,cell(i),lflag,ntok,ityp,fvalue) if (lflag .ne. 0 .and. i .le. 3) go to 199 110 continue lcell = .true. C........................................................................ elseif (key .eq. 'ORTH' .or. key .eq. 'NCOD') then C ORTHOGONALIZATION or NCODE, read orthogonalization code call gttint(2,ncode,lflag,ntok,ityp,fvalue) C........................................................................ elseif (key .eq. 'SPAC') then C SPACEGROUP, read spacegroup name for output spgnam = line(ibeg(2):iend(2)) call ccpupc(spgnam) lspnam = .true. C........................................................................ elseif (key .eq. 'SYMG') then C SYMGEN read symmetry operations to be applied call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'NCS') then c NCS keyword, set flag for non-crystallographic symmetry expansion lsymgn = -1 else jtok = 2 call rdsymm(jtok,line,ibeg,iend,ityp,fvalue,ntok, . spgnam,numsgp,pgname,nsym,nsymp,rsym) if (spgnam .ne. ' ') lspnam = .true. lsymgn = +1 endif C........................................................................ elseif (key .eq. 'SYMM') then write (6, '(a/a)') $ ' *** SYMMETRY command withdrawn ***', $ ' *** Use SYMGEN command: read write-up ***' C........................................................................ elseif (key .eq. 'CHAI') then C CHAIN set chain ID to override input n = 2 call ccpupc(cvalue(2)) k = 0 if (cvalue(2) .eq. 'SYMM') then call gttint(n+1,k,lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 if (k .gt. maxsym .or. k .le. 0) then write (6, '(a)') ' ** Illegal symmetry number **' ierr = ierr+1 go to 100 endif n = n+2 endif if (ntok .lt. n) then write (6, '(a)') ' ** Missing chain ID **' ierr = ierr+1 go to 100 elseif (ntok .eq. n) then chnnew = line(ibeg(n):iend(n)) chnold = '*' else chnold = line(ibeg(n):iend(n)) chnnew = line(ibeg(n+1):iend(n+1)) endif C C Index chain rename instructions on old chain name if (nchain .gt. 0) then do 120, j = 1, nchain if (chnold .eq. oldchn(j)) go to 121 120 continue endif C new entry nchain = nchain+1 if (nchain .gt. maxchn) then write (6, '(a)') ' ** Too many chain commands **' ierr = ierr+1 go to 100 endif C set defaults for new entry j = nchain do 122, i = 1, maxsym newchn(j, i) = chnold 122 continue C 121 oldchn(j) = chnold if ( k .gt. 0) then newchn(j, k) = chnnew else do 125, k = 1, maxsym newchn(j, k) = chnnew 125 continue endif C........................................................................ elseif (key .eq. 'BFAC') then C BFACTOR, set B-factor value C jbfac = 0 no change C = 1 reset always (ALWAYS) C = 2 reset if = 0.0 (DEFAULT) C = 3 reset if .lt. given value (MINIMUM) C = 4 reset if greater than given value (MAXIMUM) C = 5 truncate to given range (RANGE) jbfac = 1 n = 2 call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'ALWA') then jbfac = 1 n = 3 elseif (cvalue(2) .eq. 'DEFA') then jbfac = 2 n = 3 elseif (cvalue(2) .eq. 'MINI') then jbfac = 3 n = 3 elseif (cvalue(2) .eq. 'MAXI') then jbfac = 4 n = 3 elseif (cvalue(2) .eq. 'RANG') then jbfac = 5 n = 3 endif call gttrea(n,bfac,lflag,ntok,ityp,fvalue) if (lflag .gt. 0) go to 199 if (jbfac .eq. 5)then n = n+1 call gttrea(n,bfacmx,lflag,ntok,ityp,fvalue) if (lflag .gt. 0) go to 199 endif C........................................................................ elseif (key .eq. 'OCCU') then C OCCUPANCY, set occupancy value C jqfac = 0 no change C = 1 reset always (ALWAYS) C = 2 reset if = 0.0 (DEFAULT) C = 3 reset if .lt. given value (MINIMUM) jqfac = 1 n = 2 call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'ALWA') then jqfac = 1 n = 3 elseif (cvalue(2) .eq. 'DEFA') then jqfac = 2 n = 3 elseif (cvalue(2) .eq. 'MINI') then jqfac = 3 n = 3 endif call gttrea(n,occ,lflag,ntok,ityp,fvalue) if (lflag .gt. 0) go to 199 C........................................................................ elseif (key .eq. 'RENU') then C RENUMBER, read residue number increment numren = numren + 1 if (numren .gt. maxren) then write (6, '(a)') ' ** Too many ranges **' ierr = ierr+1 go to 100 endif n = 2 call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'INCR') then lincr(numren) = .true. n = 3 endif incren(numren) = 0 call gttint(n,incren(numren),lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 nrsren(1,numren) = -1000000 nrsren(2,numren) = +1000000 if (ntok .gt. n .and. ityp(n+1) .eq. 2) then C residue range call gttint(n+1,nrsren(1,numren),lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+2 call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'TO') n = n+1 call gttint(n,nrsren(2,numren),lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 endif 129 n = n+1 chnren(numren) = allchn chnrnw(numren) = chnren(numren) if (ntok .ge. n) then call ccpupc(cvalue(n)) if (cvalue(n) .ne. 'CHAI') then write (6, '(a)') ' ** Missing keyword CHAIN **' ierr = ierr+1 go to 100 else n = n+1 if (ntok .ge. n) then C chain chnren(numren) = cvalue(n)(1:1) chnrnw(numren) = chnren(numren) if (ntok .gt. n) then n = n+1 call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'TO') then n = n+1 if (ntok .ge. n) then chnrnw(numren) = cvalue(n)(1:1) else write (6, '(a)') ' ** Missing CHAIN value **' ierr = ierr+1 endif endif endif else write (6, '(a)') ' ** Missing CHAIN value **' ierr = ierr+1 endif endif endif C........................................................................ elseif (key .eq. 'SELE') then C SELECT n = 2 call ccpupc(cvalue(n)) 131 if (n .le. ntok) then if (cvalue(n) .eq. 'CHAI') then k = 0 do 130, i = n+1, ntok k = k+1 if (k .gt. maxchs) then write (6, '(a)') ' ** Too many SELECT chains**' ierr = ierr+1 go to 100 endif selchn(k) = line(ibeg(i):iend(i)) 130 continue nselch = k elseif (cvalue(n) .eq. 'OCCU') then lqmin = .true. if (n .lt. ntok) then call gttrea(n+1,qmin,lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 endif elseif (cvalue(n) .eq. 'BFAC') then lbfmax = .true. if (n .lt. ntok) then call gttrea(n+1,bfmax,lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 endif endif n = n+1 go to 131 endif C........................................................................ elseif (key .eq. 'PICK') then C PICK select only certain atom names if (ntok .gt. 1) then k = 0 do 135, i = 2, ntok k = k+1 if (k .gt. maxpat) then write (6, '(a)') ' ** Too many PICK atoms **' ierr = ierr+1 go to 100 endif selatm(k) = line(ibeg(i):iend(i)) 135 continue nselat = k endif C........................................................................ elseif (key .eq. 'ROTA') then C ROTATE define rotation matrix as MATRIX, EULER, or POLAR ltrnsf = .true. n = 1 k = 1 143 n = n+1 if (n .le. ntok) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'INVE') then C Set invert flag if inverted transformation is required linvrt = .true. else if (cvalue(n) .eq. 'MATR') then k = 1 n = n+1 elseif (cvalue(n) .eq. 'EULE') then k = 2 n = n+1 elseif (cvalue(n) .eq. 'POLA') then k = 3 n = n+1 endif if (k .eq. 1) then C Read 9 matrix elements do 140, j = 1,3 do 141, i = 1,3 call gttrea(n, trnsfm(j,i), lflag,ntok,ityp $ ,fvalue) if (lflag .ne. 0) go to 199 n = n+1 141 continue 140 continue else C Read 3 angles do 145, i = 1,3 call gttrea(n, angles(i), lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 145 continue C Make matrix if (k .eq. 2) then C Euler angles call eulmat(rmat,angles) elseif (k .eq. 3) then C Polar angles call polmat(rmat,angles) endif call strotn(rmat, trnsfm) endif n = ntok endif go to 143 endif C........................................................................ elseif (key .eq. 'SHIF') then C SHIFT define translation ltrnsf = .true. n = 1 153 n = n+1 if (n .le. ntok) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'INVE') then C Set invert flag if inverted transformation is required linvrt = .true. else if (cvalue(n) .eq. 'FRAC') then C If keyword FRACTIONAL present, set flag that translation is fractional lfrtrn = .true. n = n+1 else lfrtrn = .false. endif C Read 3 vector elements do 150, i = 1,3 call gttrea(n, trnsfm(i,4), lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 150 continue n = ntok endif go to 153 endif C........................................................................ elseif (key .eq. 'TRAN') then C TRANSFORM define 12 elements of transformation ltrnsf = .true. n = 1 163 n =n+1 if (n .le. ntok) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'INVE') then C Set invert flag if inverted transformation is required linvrt = .true. elseif (cvalue(n) .eq. 'ODB') then C Read transformation from O database file if (ntok .gt. n) then filnam = line(ibeg(n+1):) else filnam = 'ODB' endif call gtrtdb(lunit, filnam, trnsfm, istat) if (istat .ne. 0) then write (6, '(a)') $ ' **** Error in reading O datablock ****' go to 199 endif lfrtrn = .false. n = ntok else if (cvalue(n) .eq. 'FRAC') then C If keyword FRACTIONAL present, set flag that translation is c fractional lfrtrn = .true. n = n+1 else lfrtrn = .false. endif C Read 9 matrix elements do 160, j = 1,3 do 161, i = 1,3 call gttrea(n, trnsfm(j,i), lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 161 continue 160 continue C Read 3 vector elements do 162, i = 1,3 call gttrea(n, trnsfm(i,4), lflag,ntok,ityp,fvalue) if (lflag .ne. 0) go to 199 n = n+1 162 continue n = ntok endif go to 163 endif C Transformation now in trnsfm if (lsymgn .lt. 0) then c If SYMGEN NCS command given previously, increment number of symmetry c operations and store nsym = nsym+1 if (nsym .gt. maxsym) then call ccperr(1,'*** Too many symmetry operations ***') endif if (linvrt) then C Invert transformation matrix if required call rbrinv(trnsfm, amf) call ccpmvr(trnsfm, amf, 16) write (6,'(/a/)') $ ' Transformation has been inverted' endif do 164, j = 1,4 do 165, i = 1,4 rfsym(i,j,nsym) = trnsfm(i,j) 165 continue 164 continue endif C........................................................................ elseif (key .eq. 'XPLO') then C XPLOR option, set flags C Subkeys: HYDRogen keep hydrogens lxplor = .true. lnodum = .true. lhydro = .true. xdummy = 9000. n = 2 call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'HYDR') then lhydro = .false. endif C........................................................................ elseif (key .eq. 'SEQU') then C SEQUENCE option, write sequence to specified file C Option: PDB write sequence in PDB SEQRES format C SINGLE write sequence as single letter code lseqnc = .true. n = 2 168 if (n .le. ntok) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'PDB') then kseqpd = +1 n = n+1 go to 168 elseif (cvalue(n) .eq. 'SING') then kseqpd = +2 n = n+1 go to 168 else seqfil = line(ibeg(n):) go to 169 endif endif 169 continue C........................................................................ elseif (key .eq. 'OUTP') then C OUTPUT options C Subkeys: XPLOR set segid field n = 2 170 if (n .le. ntok) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'XPLO') then loutxp = .true. endif n = n+1 go to 170 endif C........................................................................ elseif (key .eq. 'UTOB') then C UTOB option, convert Us to B lutob = .true. C........................................................................ elseif (key .eq. 'ELEM') then C ELEMENT option, list 2-character element names to left-justify C according to PDB rules if (ntok .gt. 1) then do 180, n = 1, ntok-1 elmnts(n) = cvalue(n+1)(1:2) if (lenstr(elmnts(n)) .ne. 2) then ierr = ierr+1 write (6, '(/a,a/)') $ ' *** ELEMENT name must have 2-characters: ', $ elmnts(n) elseif (elmnts(n) .eq. 'CA') then ierr = ierr+1 write (6, '(/a/)') $ ' *** ELEMENT name cannot be CA,'// $ ' this clashes with Calpha ***' endif 180 continue nelmnt = n-1 endif C........................................................................ elseif (key .eq. 'REOR') then lreort = .true. n = 2 k = 2 190 if (n .le. ntok) then if (ityp(n) .eq. 2) then c Number if (k .eq. 1) then ncdin = nint(fvalue(n)) else ncdout = nint(fvalue(n)) endif elseif (ityp(n) .eq. 1) then call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'FROM') then k = 1 elseif (cvalue(n) .eq. 'TO') then k = 2 else write (6, '(/a/)') ' *** Illegal subkey ***' ierr = ierr+1 endif endif n = n+1 go to 190 endif C........................................................................ elseif (key .eq. 'REPL') then n = 2 k = 0 if (ntok .lt. 5) then write (6, '(/a/)') $ ' *** Too few keywords in REPLACE ***' go to 199 endif nreplc = nreplc+1 if (nreplc .gt. maxrtp) then write (6, '(/a/)') $ ' *** Too many REPLACE commands ***' go to 199 endif lflrtp(nreplc) = 0 resatp(nreplc) = ' ' 194 if (n .le. ntok) then call ccpupc(cvalue(n)) if (n .eq. 2) then if (cvalue(n) .eq. 'RESI') then lflrtp(nreplc) = 1 elseif (cvalue(n) .eq. 'ATOM') then lflrtp(nreplc) = 2 else go to 195 endif elseif (n .eq. 3) then c Note that residue types in this program begin with space, then up to c 3 characters, but atomnames are 4 characters tname = line(ibeg(n):iend(n)) if (lflrtp(nreplc) .eq. 1) then c Right-justify residue name call rjstfy(oldrtp(nreplc), tname) else oldrtp(nreplc) = tname endif elseif (n .eq. 4) then if (cvalue(n) .ne. 'BY') then go to 195 endif elseif (n .eq. 5) then tname = line(ibeg(n):iend(n)) if (lflrtp(nreplc) .eq. 1) then call rjstfy(newrtp(nreplc), tname) else newrtp(nreplc) = tname endif elseif (n .eq. 6) then if (cvalue(n) .ne. 'IN') then go to 195 endif elseif (n .eq. 7) then resatp(nreplc) = ' '//line(ibeg(n):iend(n)) endif n = n+1 go to 194 endif go to 100 195 write (6, '(/a/)') $ ' *** Illegal syntax in REPLACE ***' go to 199 else if (key .eq. 'COM') then com = .true. call qopen(iunscr, 'DISKIO1', 'SCRATCH') call qmode(1, 2, nbi) go to 100 C........................................................................ elseif (key .eq. 'EXCL') then C EXCLUDE C Options: C HYDROGENS n = 2 call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'HYDR') then lhydro = .true. elseif (cvalue(n) .eq. 'HEAD') then lnohed = .true. endif C........................................................................ else write (6, '(/'' Unknown keyword: '',A4/)') key endif go to 100 C........................................................................ 199 ierr = ierr+1 go to 100 C........................................................................ C C========================================================================== C Check & reflect commands C========================================================================== 200 if (ierr .gt. 0) then call ccperr(1,'** Input error **') endif C lngnmi = .false. lngnmo = .false. C C---- Get orthogonalization matrix from file if present lorthg = .false. if (.not. lcell) then call getorm(iun, rf, ro, cell, lflag) if (lflag .lt. 0) then C Error call ccperr(1,'** Error in reading PDB headers **') elseif (lflag .gt. 0) then C Set flag to say that orthogonalization matrix is available lorthg = .true. C---- If spacegroup name changed, and cell present, force new headers if (lspnam .and. mod(lflag, 2) .eq. 1) lnewhd = .true. if (lflag .ne. 2) lcell = .true. endif else call ortmat(cell, ncode, ro, rf) lorthg = .true. lnewhd = .true. endif C C---- Reorthogonalization option if (lreort) then if (.not.lcell) then call ccperr(1,'*** CELL must be given ***') endif c Input orthogonalization if (ncdin .le. 0) then c not specified, has it been read from file? if (lorthg) then c yes, copy matrices call ccpmvr(roin, ro, 16) call ccpmvr(rfin, rf, 16) ncdin = 1 else c No call ortmat(cell, ncode, roin, rfin) ncdin = ncode endif else call ortmat(cell, ncdin, roin, rfin) endif c Output orthogonalization if (ncdout .le. 0) then write (6, '(/a/)') $ '*** No output REORTHOGONALIZATION given ***' call ccperr(1,'*** REORTHOGONALIZATION TO ??? ***') endif call ortmat(cell, ncdout, roout, rfout) ncode = ncdout lnewhd = .true. c write (6, '(/a,a,a,a/)') $ ' Coordinates will be reorthogonalized from ', $ codes(ncdin),' to ', codes(ncdout) endif C C---- New orthogonalization for header if (lnewhd) then write (6,'(/a)') ' CRYST1 & SCALE header will be written' write (6,'(7x,a,i5/7x,a,a)') $ ' Orthogonalization code = ',ncode, $ ' Orthogonal axes x,y,z along ', codes(ncode) C make PDB header C write (iunout,6000) cell,spgnam 6000 format('CRYST1',3f9.3,3f7.2,1x,a11) write (iunout,6001) (i,(rf(i,j),j=1,3),i=1,3) 6001 format(2('SCALE',i1,4x,3f10.6,5x,' 0.00000'/), $ ('SCALE',i1,4x,3f10.6,5x,' 0.00000')) endif C C---- EXCLUDE HEADERS option if (lnohed) write (6,'(/a)') $ ' Header lines will not be copied from input file' C C---- B-factor & occupancy reset if (jbfac .eq. 1) then write (6, '(/a,f8.2,a)') ' B-factors will be reset to ',bfac else if (jbfac .eq. 2) then write (6, '(/a,f8.2,a)') ' B-factors will be reset to ',bfac, $ ', ONLY if zero' else if (jbfac .eq. 3) then write (6, '(/a,f8.2,a)') ' B-factors will be reset to ',bfac, $ ', if less than this value' else if (jbfac .eq. 4) then write (6, '(/a,f8.2,a)') ' B-factors will be reset to ',bfac, $ ', if greater than this value' else if (jbfac .eq. 5) then write (6, '(/a,f8.2,a,f8.2,a)') $ ' B-factors will be reset to between', $ bfac,' and', bfacmx, $ ', if outside this range' endif if (jqfac .eq. 1) then write (6, '(/a,f8.2,a)') ' Occupancy will be reset to ',occ else if (jqfac .eq. 2) then write (6, '(/a,f8.2,a)') ' Occupancy will be reset to ',occ, $ ', ONLY if zero' else if (jqfac .eq. 3) then write (6, '(/a,f8.2,a)') ' Occupancy will be reset to ',occ, $ ', if less than this value' endif if (lutob) then write (6,'(/a)') $ ' Input Us will be converted to Bs' endif C C-----Renumbering if (numren .gt. 0) then do 220, i=1,numren if (lincr(i)) then write (6,'(/a,i5,2(a,i8),a,a,a)') $ ' Increment residue numbers by ',incren(i), $ ' for residue range ', nrsren(1,i), ' to',nrsren(2,i), $ ' in chain ', chnren(i), '(* == all chains)' else write (6,'(/a,i5,2(a,i8),a,a,a)') $ ' Renumber residue numbers from ',incren(i), $ ' for residue range ', nrsren(1,i), ' to',nrsren(2,i), $ ' in chain ', chnren(i), '(* == all chains)' endif if (chnrnw(i) .ne. chnren(i)) then write (6, '(5x,a,a1)') 'Chain renamed to ',chnrnw(i) endif if (lenstr(chnren(i)) .gt. 1) lngnmi = .true. if (lenstr(chnrnw(i)) .gt. 1) lngnmo = .true. 220 continue endif C C---- Chain renaming jchwld = 0 if (nchain .gt. 0) then if (nsym .eq. 0) then write (6,'(/13X,A,5X,A)') 'Old chain','New chain' else write (6,'(/13X,A,5X,A)') 'Old chain','Symmetry chains' endif k = max(1, nsym) do 230, j = 1, nchain if (lenstr(oldchn(j)) .gt. 1) lngnmi = .true. do 231, i=1,k if (lenstr(newchn(j,i)) .gt. 1) lngnmo = .true. 231 continue write (6, '((3X,A,3X,A,96(2X,A)))') $ 'Rename chain',oldchn(j),' to ',(newchn(j,i),i=1,k) if (oldchn(j) .eq. '*') jchwld = j 230 continue endif C C---- Symmetry c lsymgn = 0 no symmetry expansion, +1 crystallographic, c -1 non-crystallographic c For NCS symmetry, matrices are already in array rfsym if (lsymgn .gt. 0) then C Must have orthogonalization matrices if (.not. lorthg) then call ccperr(1, $ '** For symmetry expansion, Cell must be defined **') endif do 240, n = 1,nsym C For each symmetry operation n, generate symmetry for C orthogonal coordinates C [rfsym] = [ro] [rsym] [rf] call matmln(4, amf, rsym(1,1,n), rf) call matmln(4, rfsym(1,1,n), ro, amf) 240 continue elseif (lsymgn .lt. 0) then c For NCS symmetry, matrices are already in array rfsym ltrnsf = .false. write (6, '(//a/)') $ ' NCS-Symmetry-related coordinates will be transformed'// $ ' as follows:' do 241, isym = 1, nsym write (6, 6241) isym, $ ((rfsym(i,j,isym), j = 1,4), i=1,3) 6241 format(' Symmetry operation ',i6/ $ 7x,'( ',3f10.6,' ) ( x ) ( ',f9.3,' )'/ $ 7x,'( ',3f10.6,' ) ( y ) + ( ',f9.3,' )'/ $ 7x,'( ',3f10.6,' ) ( z ) ( ',f9.3,' )'/) 241 continue endif C C---- Check transformation if (ltrnsf) then C C---- If translation was fractional, convert to orthogonal if (lfrtrn) then if (lorthg) then c Multiply translation vector by fractional to orthogonal matrix ro call trnfrm(trnsfm(1,4), ro) else write (6, '(/a,a/)') $ ' **** No orthogonalization matrix available', $ ' for fractional translation ****' call ccperr(1,'** Missing orthogonalization **') endif endif C if (linvrt) then C Invert transformation matrix if required call rbrinv(trnsfm, amf) call ccpmvr(trnsfm, amf, 16) endif C call tstrot(trnsfm, det, lflag) if (lflag .ne. 0) then write (6,'(A/A, F12.5)') $ ' **** Transformation matrix is not a rotation ****', $ ' Determinant = ',det if (lflag .lt. 0) then write (6,'(A)') ' **** Matrix reverses hand ****' endif call ccperr(1,'** Illegal transformation **') endif c if (linvrt) $ write (6,'(/a/)') $ ' Transformation has been inverted' write (6, 6201) ((trnsfm(i,j), j = 1,4), i=1,3) 6201 format(/' Coordinates will be transformed as follows:'// $ 7x,'( ',3f10.6,' ) ( x ) ( ',f9.3,' )'/ $ 7x,'( ',3f10.6,' ) ( y ) + ( ',f9.3,' )'/ $ 7x,'( ',3f10.6,' ) ( z ) ( ',f9.3,' )'/) endif C C C---- Selection if (nselch .gt. 0) then write (6,'(/A,50(2X,A))') $ ' Select only chains : ',(selchn(i), i=1,nselch) do 245, i=1,nselch if (lenstr(selchn(i)) .gt. 1) lngnmi = .true. 245 continue endif if (lqmin) then write (6,'(/A,F10.4)') $ ' Select only atoms with occupancy .gt. ',qmin endif if (lbfmax) then write (6,'(/A,F10.4)') $ ' Select only atoms with Bfactor .lt. ',bfmax endif C C---- Pick if (nselat .gt. 0) then write (6,'(/A,20(2X,A))') $ ' Select only atoms : ',(selatm(i), i=1,nselat) endif C C---- Xplor if (lxplor) then write (6,'(a)') $ ' Input file is Xplor type, remove dummy atoms' endif if (lhydro) then write (6,'(a)') $ ' Hydrogen atoms will be removed' else write (6,'(a)') $ ' Hydrogen atoms will be kept' endif if (loutxp) then write (6, '(a)') $ ' Chain identifier will be written as Xplor SegID' endif c---- Element if (nelmnt .gt. 0) then do 250, i = 1, nelmnt write (6, '(a,a)') $ ' Left-justify atoms of element ',elmnts(i) 250 continue endif C if (lngnmi .and. .not.lxplor) then call ccperr(1, $ '*** Input chain names can only be longer than'// $ ' 1 character with Xplor input ***') endif if (lngnmo .and. .not.loutxp) then call ccperr(1, $ '*** Output chain names can only be longer than'// $ ' 1 character with Xplor output ***') endif c if (nreplc .gt. 0) then write (6, '(/a)') $ ' Replace residue types or atom names' do 251, i = 1, nreplc if (lflrtp(i) .eq. 1) then write (6, '(5x,a,a,a,a)') $ 'Replace ', oldrtp(i),' by ',newrtp(i) elseif (lflrtp(i) .eq. 2) then write (6, '(5x,a,a,a,a,a,a)') $ 'Replace ', oldrtp(i),' by ',newrtp(i), $ ' in ', resatp(i) endif 251 continue endif C C========================================================================== C Do it C========================================================================== jres = -1000000 natom = 0 chnold = '?' nqmin = 0 C C------------------------------------------------------------- C Symmetry loop isym = 1 300 continue C------------------------------------------------------------- C Atom loop 400 continue C Don't use RBROOK to read atoms, because residue names from O screw it up read (iun,'(a)',end=500) atline if (atline(1:4) .eq. 'ATOM' .or. atline(1:4) .eq. 'HETA') then if (lxplor) then C Special read option for Xplor output, residue numbers .ge. 1000 are in C different column if (atline(27:27) .eq. ' ') then read (atline,fmt=6102) $ iser,atnam,restyp,resno,x,q,b,segid 6102 format(6x,i5,1x,a4,a4,2x,a4,1x,3x,3f8.3,2f6.2,6x,a) else read (atline,fmt=6103) $ iser,atnam,restyp,resno,x,q,b,segid 6103 format(6x,i5,1x,a4,a4,3x,a4,3x,3f8.3,2f6.2,6x,a) endif idch = segid(1:1) c Make chain ID uppercase call ccpupc(idch) insert = ' ' else c Normal input read (atline,fmt=6101) $ iser,atnam,restyp,idch,resno,insert,x,q,b 6101 format(6x,i5,1x,a4,a4,1x,a1,a4,a1,3x,3f8.3,2f6.2) segid = idch endif c c---- Omit hydrogens if required if (lhydro) then if (atnam(1:1).eq.'H' .or. atnam(2:2).eq.'H') go to 400 endif c c---- Omit dummys if required if (lnodum) then if (abs(x(1)) .gt. xdummy) then go to 400 endif endif c c---- Atom accepted C C---- Split residue number if required C If resno begins with non-numeric character, split into chain-ID & number do 420, j=1,4 if (resno(j:j) .ne. ' ') go to 421 420 continue C blank residue nnumber, stop call ccperr(1,' *** Residue number blank ***') C 421 if (lge(resno(j:j),'0').and.lle(resno(j:j),'9')) then C first non-blank character is digit, so leave label alone continue else C not digit, so extract as chain ID idch = resno(j:j) segid = idch resno(j:j) = ' ' endif C C---- List of chains in input file (after unpacking residue numbers) if (segid .ne. ' ' .and. isym .eq. 1) then if (mchain .gt. 0) then do 402, j = 1,mchain if (segid .eq. chains(j)) go to 403 402 continue endif mchain = mchain+1 if (mchain .le. maxchn) then chains(mchain) = segid endif 403 continue endif C C---- Select options if (nselch .gt. 0) then do 410, j = 1, nselch if (segid .eq. selchn(j)) go to 411 410 continue go to 400 endif 411 continue if (lqmin) then if (q .le. qmin) then nqmin = nqmin+1 go to 400 endif endif if (lbfmax) then if (b .gt. bfmax) then nbfmax = nbfmax+1 go to 400 endif endif C---- Pick if (nselat .gt. 0) then do 415, j = 1, nselat if (eqlstr(atnam, selatm(j))) go to 416 415 continue go to 400 endif 416 continue C C---- Option to fix up element component of atomname if (nelmnt .gt. 0) then call jfyelm(atnam, elmnts, nelmnt) endif C C---- Option to convert Us to Bs if (lutob) then b = utob*b endif C C---- Reset Bfactor & occupancy C Bfactor option, set B & occupancy if (jbfac .eq. 1) then b = bfac elseif (jbfac .eq. 2) then if (b .eq. 0.0) b = bfac elseif (jbfac .eq. 3) then b = max(b, bfac) elseif (jbfac .eq. 4) then b = min(b, bfac) elseif (jbfac .eq. 5) then b = min(max(b, bfac), bfacmx) endif if (jqfac .eq. 1) then q = occ elseif (jqfac .eq. 2) then if (q .lt. 0.000001) q = occ elseif (jqfac .eq. 3) then q = max(q,occ) endif C---- Is this a new chain? lnewch = .false. if (segid .ne. chnold) then if (chnold .ne. '?') then lnewch = .true. endif chnold = segid endif C C---- Rename chains C Override chain-ID if CHAIN command given if (nchain .gt. 0) then do 430, j = 1, nchain if (segid .eq. oldchn(j)) go to 431 430 continue C this chain not explicitly in rename list, is there a wild-card rename? j = 0 if (jchwld .gt. 0) j = jchwld 431 if (j .gt. 0) then C Yes, there is a rename command active for this chain, so do it, C unless this would rename to '*' if (newchn(j, isym) .ne. allchn) then segid = newchn(j, isym) endif endif endif C C---- Pick up residue number C residue number is in resno read (resno, '(i4)') ires if (jres .ne. ires .or. insert.ne.insold) then newres = .true. jres = ires insold = insert else newres = .false. endif C C---- Renumber option if (numren .gt. 0) then do 440, j=1,numren C Renumber option if (ires .ge. nrsren(1,j) .and. $ ires .le. nrsren(2,j)) then if (chnren(j) .eq. allchn .or. $ chnren(j) .eq. segid) then C Yes, renumber if (lincr(j)) then ires = ires + incren(j) else if (newres) then kres = incren(j) incren(j) = incren(j) + 1 endif ires = kres insert = ' ' endif write (resno, '(i4)') ires if (chnren(j) .ne. allchn) segid = chnrnw(j) go to 445 endif endif 440 continue endif 445 continue c c---- Reorthogonalization option if (lreort) then c coordinate -> fractional by input orthogonalization call trnfrm(x, rfin) c orthogonalize with output orthogonalization call trnfrm(x, roout) endif C C---- Coordinate transformations if (nsym .gt. 0) then C Symmetry generation call trnfrm(x, rfsym(1,1,isym)) endif if (ltrnsf) then call trnfrm(x, trnsfm) endif C do 450, i = 1,3 xyzmin(i) = min(xyzmin(i), x(i)) xyzmax(i) = max(xyzmax(i), x(i)) 450 continue C c---- Residue replace if (nreplc .gt. 0) then do 455, i = 1, nreplc if (lflrtp(i) .eq. 1) then if (restyp .eq. oldrtp(i)) then restyp = newrtp(i) endif elseif (lflrtp(i) .eq. 2) then if (resatp(i) .eq. ' ' $ .or. (resatp(i) .eq. restyp)) then if (atnam .eq. oldrtp(i)) then atnam = newrtp(i) endif endif endif 455 continue endif C---- Sequence if (lseqnc) then if (newres) then nseqn = nseqn + 1 if (nseqn .gt. maxseq) then write (6, '(a,a,i8)') $ ' *** Too many residues in sequence, ', $ 'increase maximum MAXSEQ =',maxseq call ccperr(1,'** Too many residues **') endif if (lnewch) then C New chain started, so put a marker residue in the sequence seqnce(nseqn) = ' END' C note mchain is new current chain C record pointer in sequence store to end of last chain & beginning of new one nsqchn(2, mchain-1) = nseqn-1 nseqn = nseqn + 1 nsqchn(1, mchain) = nseqn if (nseqn .gt. maxseq) then write (6, '(a,a,i8)') $ ' *** Too many residues in sequence, ', $ 'increase maximum MAXSEQ =',maxseq call ccperr(1,'** Too many residues **') endif endif seqnce(nseqn) = restyp endif endif C---- Chain ID must be uppercase idch = segid(1:1) call ccpupc(idch) C---- If not Xplor output, clear segid if (.not.loutxp) segid = ' ' C---- Write out natom = natom+1 write (iunout,fmt=6111) $ natom,atnam,restyp,idch,resno,insert,x,q,b,segid 6111 format('ATOM ',i5,1x,a4,a4,1x,a1,a4,a1,3x,3f8.3,2f6.2,6X,A) C write coordinates to scratch file for another pass if (com) then call qwritr(1,x,3) xcom (1) = xcom (1) + x (1) xcom (2) = xcom (2) + x (2) xcom (3) = xcom (3) + x (3) end if else if (atline(1:3) .eq. 'TER') then C---- Always copy TER line write (iunout,'(a)') atline elseif (isym .eq. 1) then if (.not. lnohed) then C---- Copy non-ATOM lines from input to output, unless EXCLUDE HEADERS set C except for CRYST1 & SCALE lines if they have been replaced if (lnewhd) then if (atline(1:5) .ne. 'SCALE' .and. $ atline(1:5) .ne. 'CRYST') then write (iunout,'(a)') atline endif else write (iunout,'(a)') atline endif endif endif endif go to 400 C---- End atom loop 500 continue if (isym .lt. nsym) then C---- Symmetry loop isym = isym+1 rewind (unit=iun) go to 300 endif C C---- End everything write (6,'(/1x,i8,a)') natom,' atoms copied' if (lqmin) then write (6,'(1x,i8,a,f8.3,a)') $ nqmin, ' atoms with occupancy less than ',qmin, $ ' omitted' endif if (lbfmax) then write (6,'(1x,i8,a,f8.3,a)') $ nbfmax, ' atoms with Bfactor greater than ',bfmax, $ ' omitted' endif do 510, i = 1,3 x(i) = 0.5 * (xyzmin(i) + xyzmax(i)) 510 continue if (natom .gt. 0) then write (6,6500) $ (xyzmin(i),xyzmax(i),x(i),xyzmax(i)-xyzmin(i), i=1,3) 6500 format(//' Coordinate limits in output file:'/ $ 7x,' Minimum Maximum Centre Range'/ $ 7x,'On X : ',4f10.2/ $ 7x,'On Y : ',4f10.2/ $ 7x,'On Z : ',4f10.2/) if (com) then call qseek (1, 1, 1, 1) maxdist = 0.0 xcom(1) = xcom(1)/natom xcom(2) = xcom(2)/natom xcom(3) = xcom(3)/natom do 511 i = 1, natom call qreadr(iunscr, xx, 3, ierr) maxdist = max (maxdist, + sqrt ((xcom(1)-xx(1))**2 + (xcom(2)-xx(2))**2 + + (xcom(3)-xx(3))**2)) 511 continue if (ierr .ne. 0) call ccperr(1, + ' *** read error in scratch file ***') write (6, 6501) xcom, maxdist C Format per amore: 6501 format (/, + ' Center of Mass:',3F9.2 / + ' Maximal distance from Center of Mass:',F9.2 /) end if if (mchain .gt. 0) then k = min(maxchn, mchain) write (6,6510) mchain, (chains(i), i = 1,k) 6510 format(/' Number of chains in input file = ',i6// $ ' Original chain IDs: ',50(2x,a)) endif write (6,'(/)') endif C---- Sequence if (lseqnc) then nsqchn(2, mchain) = nseqn write (6, '(a,a)') ' Sequence will be written to file ', $ seqfil(1:lenstr(seqfil)) call ccpdpn(9, seqfil, 'NEW', 'F', 0, 0) if (kseqpd .eq. +1) then C Write sequence as PDB SEQRES, each chain separately do 520, n = 1, mchain write (9, 6600) chains(n) 6600 format ('REMARK Chain ',A) k = nsqchn(1,n) j = 1 m = nsqchn(2,n) - nsqchn(1,n) + 1 521 if (k .le. nsqchn(2,n)) then write (9,6601) j, chains(n), m, $ (seqnce(i), i = k, min(k+12, nsqchn(2,n))) 6601 format ('SEQRES',i4,1x,a1,1x,i4,1x,13a4) j = j+1 k = k+13 go to 521 endif 520 continue elseif (kseqpd .eq. +2) then C Write sequence as single-letter code, each chain separately c Convert to single letter code do 530, i = 1, nseqn seqnc1(i) = ' ' call res3to1(seqnce(i)(2:4), seqnc1(i)) 530 continue do 535, n = 1, mchain write (9, 6610) chains(n) 6610 format ('>chain_',A) k = nsqchn(1,n) j = 1 m = nsqchn(2,n) - nsqchn(1,n) + 1 536 if (k .le. nsqchn(2,n)) then write (9,6611) $ (seqnc1(i), i = k, min(k+59, nsqchn(2,n))) 6611 format (6(1x,10a1)) j = j+1 k = k+60 go to 536 endif 535 continue else c or approximately O format write (9, '(1000(10(A4,2X)/))') (seqnce(i), i = 1,nseqn) endif endif C close(unit=iun,status='keep') close(unit=iunout,status='keep') if (lseqnc) close(unit=9,status='keep') call ccperr(0,' === Normal completion PDBSET ===') C end C C subroutine polmat(rotn,angles) C ============================== C C************************************************************** C Sets rotation matrix ROTN expressing a rotation through an C angle KAPPA right-handedly about an axis C with polar angles OMEGA, PHI (OMEGA with Z-axis, projection C giving angle PHI with X-axis). C These angles give direction cosines DIRCOS(I). C C Input: C ANGLES angles omega, phi, kappa in degrees C C Output: C ROTN (3,3) rotation matrix C C*************************************************************** C C real rotn(3,3),dircos(3),angles(3) real omega,phi,kappa C real conv,snom,csom,snph,csph,snka,cska,epsijk integer i,j,k,k1 C conv=atan(1.)/45. omega=angles(1) phi=angles(2) kappa=angles(3) snom=sin(omega*conv) csom=cos(omega*conv) snph=sin(phi*conv) csph=cos(phi*conv) snka=sin(kappa*conv) cska=cos(kappa*conv) dircos(1)=snom*csph dircos(2)=snom*snph dircos(3)=csom do 40 i=1,3 do 30 j=1,3 k1=6-i-j k=k1 if ((k1.lt.1).or.(k1.gt.3)) k=3 epsijk=((i-j)*(j-k)*(k-i))/2 rotn(i,j)=dircos(i)*dircos(j)* $ (1.0-cska)-epsijk*dircos(k)*snka if (i.eq.j) rotn(i,j)=rotn(i,j)+cska 30 continue 40 continue return end C C subroutine eulmat(rotn,angles) C ============================== C C Sets rotation matrix ROTN expressing a rotation through an C angles alpha, beta, gamma (degrees) = ANGLES C C Input: C ANGLES alphas, beta, gamma in degrees C C Output: C ROTN (3,3) rotation matrix C C C*************************************************************** C real rotn(3,3),angles(3) real alpha,beta,gamma real conv,sina,cosa,sinb,cosb,sing,cosg C alpha=angles(1) beta=angles(2) gamma=angles(3) conv=atan(1.0)/45. sina=sin(alpha*conv) cosa=cos(alpha*conv) sinb=sin(beta*conv) cosb=cos(beta*conv) sing=sin(gamma*conv) cosg=cos(gamma*conv) rotn(1,1)=cosg*cosb*cosa-sing*sina rotn(2,1)=cosg*cosb*sina+sing*cosa rotn(3,1)=-cosg*sinb rotn(1,2)=-sing*cosb*cosa-cosg*sina rotn(2,2)=-sing*cosb*sina+cosg*cosa rotn(3,2)=sing*sinb rotn(1,3)=sinb*cosa rotn(2,3)=sinb*sina rotn(3,3)=cosb return end C C C subroutine strotn(a, b) C ====================== C C STore ROTatioN C Transfer rotational elements of 3x3 matrix A to 4x4 matrix B C Other elements of B left unchanged C real a(3,3), b(4,4) C integer i,j C do 1, j = 1,3 do 2, i = 1,3 b(i,j) = a(i,j) 2 continue 1 continue return end C C C subroutine getorm(iun, rf, ro, cell, lflag) C =========================================== C C Read Brookhaven file open on stream iun, search for CRYST1 & SCALE lines C Return when first ATOM line found C C Input: C iun unit number for reading C C Output: C rf(4,4) orthogonal to fractional matrix C ro(4,4) fractional to orthogonal matrix C cell(6) unit cell (only if present on file) C lflag = 0 if neither cell nor matrix found C .lt. 0 if error C = 1 cell only C = 2 Scale only (no cell returned) C = 3 both found, matrices taken from SCALE lines C C C integer iun, lflag real rf(4,4), ro(4,4), cell(6) C character line*80 integer nsc, i, j real vol, rall(3,3,6) C character*1 c123(3) data c123 /'1','2','3'/ C nsc = 0 lflag = 0 C 10 read (iun, '(a)') line if (line(1:4) .eq. 'ATOM' .or. line(1:6) .eq. 'HETATM') then go to 100 elseif (line(1:6) .eq. 'CRYST1') then C CRYST1 line, set cell read (line, 6001) cell 6001 format(6x,3f9.3,3f7.2) lflag = lflag + 1 elseif (line(1:5) .eq. 'SCALE') then C SCALE line, set rf matrix nsc = nsc+1 if (line(6:6) .ne. c123(nsc)) then lflag = -1 return endif read (line, 6002) (rf(nsc,j), j=1,4) 6002 format(6x,4x,3f10.5,5x,f10.5) if (nsc .eq. 3) then do 20, j = 1,3 rf(4,j) = 0.0 20 continue rf(4,4) = 1.0 lflag = lflag+2 endif endif go to 10 C C------- C All header lines read, process matrices C If cell was read, make matrix from cell 100 if (lflag .eq. 1) then call rbfro1(cell,vol,rall) do 210, j=1,4 do 211, i=1,4 ro(i,j)=0.0 if (i.lt.4.and. j.lt.4) ro(i,j)=rall(i,j,1) 211 continue 210 continue ro(4,4)=1.0 C Invert call rbrinv(ro,rf) C elseif (lflag .gt. 0) then call rbrinv(rf, ro) endif C rewind (unit=iun) return end C C C subroutine matmln(n,a,b,c) C ========================= C C Multiply two nxn matrices C a = b . c C integer n real a(n,n),b(n,n),c(n,n) integer i,j,k C do 1 i=1,n do 2 j=1,n a(j,i)=0. do 3 k=1,n a(j,i)= b(j,k)*c(k,i)+a(j,i) 3 continue 2 continue 1 continue return end C C C subroutine trnfrm(x, rtmat) C =========================== C C Transform vector x by 4x4 matrix rtmat (assumed to be essentially 3x4) C real x(3), rtmat(4,4) C real y(3) integer i,j C C do 1, j = 1,3 y(j) = rtmat(j,4) do 2, i=1,3 y(j) = y(j) + rtmat(j,i) * x(i) 2 continue 1 continue C do 3, i=1,3 x(i) = y(i) 3 continue return end C C C subroutine tstrot(rmat, det, lflag) C ================================== C C Test if 3x3 part of 4x4 matrix is actually a rotation C C Input: C rmat(4,4) matrix to be tested C C Output: C det determinant C lflag = 0 OK C = 1 or 3 C C real rmat(4,4), det integer lflag C real amat(3,3), errlim, test integer i,j,k logical fail C data errlim/0.0001/ C lflag = 0 C det = rmat(1,1)*rmat(2,2)*rmat(3,3) $ + rmat(1,2)*rmat(2,3)*rmat(3,1) $ + rmat(1,3)*rmat(2,1)*rmat(3,2) $ - rmat(1,3)*rmat(2,2)*rmat(3,1) $ - rmat(1,2)*rmat(2,1)*rmat(3,3) $ - rmat(1,1)*rmat(2,3)*rmat(3,2) C if (abs(det-1.0) .gt. errlim) then lflag = 1 endif C C Multiply by transpose do 1, i=1,3 do 2, j=1,3 amat(i,j) = 0.0 do 3, k=1,3 amat(i,j) = amat(i,j) +rmat(i,k)*rmat(j,k) 3 continue 2 continue 1 continue C C Test for identity fail = .false. do 11, i=1,3 do 12, j=1,3 if (i .eq. j) then test = 1.0 else test = 0.0 endif if (abs(amat(i,j) - test) .gt. errlim) fail = .true. 12 continue 11 continue C if (fail) lflag = lflag + 2 C if (det .lt. errlim) then C Negative determinant lflag = - lflag endif C return end c c c SUBROUTINE RJSTFY(A,B) C ====================== C C Right justify string b into a C CHARACTER*(*) A,B INTEGER I,J,LA,LB,MA,LENSTR EXTERNAL LENSTR C LA=LEN(A) LB=LENSTR(B) A = ' ' C MA = MAX(1,LA-LB+1) J = LB DO 1, I=LA,MA,-1 A(I:I) = B(J:J) J = J-1 1 CONTINUE RETURN END c c c logical function eqlstr(a,b) c ============================ c c Returns .true. if strings a & b are equal, after justification c character*(*) a, b c character*80 c, d c call rjstfy(c, a) call rjstfy(d, b) if (c .eq. d) then eqlstr = .true. else eqlstr = .false. endif return end subroutine gtrtdb(lunit, file, trnsfm, istat) c ============================================= c c Read an O datablock file to get transformation c c trnsfm is a 4x4 matrix which premultiplies (x y z 1)T c c ( x' ) ( 1,1 1,2 1,3 1,4 ) ( x ) c ( y' ) = ( 2,1 2,2 2,3 2,4 ) ( y ) c ( z' ) ( 3,1 3,2 3,3 3,4 ) ( z ) c ( 1 ) ( 4,1 4,2 4,3 4,4 ) ( 1 ) c c ( R1 R4 R7 R10:tx ) ( x ) c = ( R2 R5 R8 R11:ty ) ( y ) c ( R3 R6 R9 R12:tz ) ( z ) c ( 0 0 0 1 ) ( 1 ) c c c On entry: c lunit unit number to use for reading c file file name c c On exit: c trnsfm(4,4) 4x4 transformation matrix c istat = 0 OK c = -1 end of file before all found c = +1 blank file c = +2 illegal data block (wrong number of items, etc) c c integer lunit, istat real trnsfm(4,4) character file*(*) c integer nitems, i, j character blknam*80, type*1, format*80 c call opnodb(lunit, file, blknam, type, nitems, format, istat) c if ((type .ne. 'R' .and. type .ne. 'r' ) .or. nitems .ne. 12) then istat = +2 endif if (istat .ne. 0) return c read (lunit, format) ((trnsfm(i,j), i=1,3), j=1,3), $ (trnsfm(i,4),i=1,3) trnsfm(4,1) = 0.0 trnsfm(4,2) = 0.0 trnsfm(4,3) = 0.0 trnsfm(4,4) = 1.0 close (unit=lunit) return c end subroutine opnodb( $ lunit, file, blknam, type, nitems, format, istat) c =================================================================== c c Open an O data block file & read header line c c On entry: c lunit unit number to open file on c file file name c c On exit: c blknam data block name c type block type c nitems number of items c format format c istat = 0 OK c = -1 end of file c c integer lunit, nitems, istat character*(*) file, blknam, type, format c C things for parser ---- integer maxtok parameter (maxtok=100) character cvalue(maxtok)*4,line*256 integer ibeg(maxtok),iend(maxtok),ityp(maxtok), . idec(maxtok),ntok real fvalue(maxtok) c c istat = 0 call ccpdpn(lunit, file, 'READONLY', 'F', 0,0) c 10 read (lunit, '(a)',end=910) line ntok=-maxtok call parse( . line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok) if(ntok.le.0) then go to 10 endif c blknam = line(ibeg(1):iend(1)) type = line(ibeg(2):iend(2)) nitems = nint(fvalue(3)) format = line(ibeg(4):iend(ntok)) c 900 return 910 istat = -1 go to 900 end c c c subroutine jfyelm(atnam, elmnts, nelmnt) c ======================================== c c Justify element name in atom name c If first 2 non-blank characters in atnam match one of the elements c in the list elmnts (2-character atomnames only), then c leftjustify atomname c c On entry: c atnam*4 atom name c elmnts*2 element names to justify c nelmnt number of element names c c On exit: c atnam modified element names c implicit none integer nelmnt character atnam*4, elmnts(nelmnt)*2 c integer i,j character*4 name c do 2, j = 1, 4 if (atnam(j:j) .ne. ' ') go to 3 2 continue return c 3 do 10, i = 1, nelmnt if (atnam(j:j+1) .eq. elmnts(i)) then c Found name = atnam(j:) atnam = name return endif 10 continue c return end c c c subroutine ortmat(cell, ncode, ro, rf) c ====================================== c c Return orthogonalization matrix c c On entry: c cell(6) cell dimensions c ncode orthogonalization code C = 1 axes along a, c* x a, c* (Brookhaven standard, default) C = 2 axes along b, a* x b, a* C = 3 axes along c, b* x c, b* C = 4 axes along a+b, c* x (a+b), c* C = 5 axes along a*, c x a*, c ( Rollett ) C = 6 axes along a, b*, a x b* C = 7 axes along a*, b, a* x b (TNT) c c On exit: c ro(4,4) orthogonalization matrix c rf(4,4) inverse of ro, fractionalization matrix c implicit none c real cell(6), ro(4,4), rf(4,4) integer ncode c real rall(3,3,6), vol integer i, j, k call rbfro1(cell,vol,rall) do 10, j=1,4 do 11, i=1,4 ro(i,j)=0.0 11 continue 10 continue c if (ncode .le. 6) then do 210, j=1,3 do 211, i=1,3 ro(i,j)=rall(i,j,ncode) 211 continue 210 continue elseif (ncode .eq. 7) then c ncode 7 for TNT, permute ncode 2, rows of matrix are 3,1,2 do 310, j=1,3 do 311, i=1,3 k = i-1 if (k .le. 0) k = k+3 ro(i,j)=rall(k,j,2) 311 continue 310 continue else write (6, '(/a,i6,a/)') $ ' *** Illegal NCODE: ',ncode,' ***' call ccperr(1,'*** Illegal NCODE ***') endif c ro(4,4)=1.0 c C Invert call rbrinv(ro,rf) c return end