C C PDBSET C Copyright (C) 1992 Phil Evans C C This code is distributed under the terms and conditions of the C CCP4 Program Suite Licence Agreement as a CCP4 Application. C A copy of the CCP4 licence can be obtained by writing to the C CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK. C C 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 new 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 two cells are diferent the transformation is: C X_out = [RO_out] [ symm] [ RF_in] X_in. C 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 RESET set q = 0.0 if .lt. Occ, else set q = 1.0 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 TRANSFORM [INVERT] FILE 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 Subkey FILE read from a formatted file containing a 3x3 matrix c followed by a translation vector 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 THIS IS NOW DONE BY DEFAULT 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 C 25) NOISE [maximum shift] C maximum shift (Angs) C defauls to 0.2 Angs, fails if greater than 0.5 Angs C Subkeys: C CHAIN select only specified chain(s) C eg NOISE 0.1 CHAIN C ! select only chain C C BFACTOR [] C select only atoms with B-factor greater than C PICK define atom names to be included: all other atoms will C be omitted C eg NOISE 0.1 PICK CA to choose C-alpha only C Note that the atomname is case-sensitive c========================================================================== c attempts to adapt to new rwbrook library (3.3 release) c Phil Evans 2/10/97 c c Problems:- c c 1) if cell changed, the library adds CRYST1 & SCALEn linesw, but copies c the old ones as well c c 2) the library at present doesn't write the spacegroup symbol into c the file c c 3) if we have both chain_ID and SegId, then we need to define what c the chain rename commands here should do. At present, chain_id will always c be the first character of segid c Phil Evans 28/6/2000 a fix of sorts c (a) for selections, use ChId unless blank .and. Segid set c (b) If segid blank, always set = ChId c c c========================================================================== 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 integer ibeg(maxtok),iend(maxtok),ityp(maxtok), . idec(maxtok),ntok real fvalue(maxtok) logical lend,anisflg C real cellpdb(6), cellnew(6), $ trnsfm(4,4), amf(4,4), rmat(3,3), $ tmp(4,4),tmpT(4,4),uf(4,4),rf_tru(4,4),ro_trU(4,4) C Selection integer maxchs parameter (maxchs=20) C number of side chain and main chain atoms integer matsc,matmc parameter (matsc=100,matmc=4) 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 ATRE, atom renumbering variables logical latren integer iserout C NOISE, random shift variables logical lnoise, lbfmin integer nnlchn, nnlatm, iseed, iseed2, nnoise real rnoise,bfmin,pxval character*4 nolatm(maxpat),nolchn(maxchs) external rmarin 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,lspgrp, + lspgrpPDB,nsymppdbs,nsympdbs real rsympdbs(4,4,maxsym) real rsym(4,4,maxsym), rfsym(4,4,maxsym) integer lrfsym(maxsym) character spgnam*10, pgname*10,SPGRPpdb*15,SPGNAM_CIFS*20, + NAMSPG_CIFS*20,nampg*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, chnnam*4 c c Element names integer maxelm,maxaa parameter (maxelm=50,maxaa=21) character*2 elmnts(maxelm) integer nelmnt C C Sequence integer maxseq parameter (maxseq=5000) character seqnce(maxseq)*4, seqfil*80, seqnc1(maxseq)*1, + seqid1(maxaa)*1,seqid3(maxaa)*4 logical lseqnc integer nseqn, nsqchn(2,maxchn), kseqpd,mwseqid(maxaa), + mwchn(maxchn),mwtot 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 character restp3*3 C integer iun,iunout,i,j,k,m,n,iser,ires, jres, kres, ierr, lflag, $ jtok, jchwld, isym, ncode, natom, lunit, istat, nbi, iunscr, $ ifail, iz, iter, icopy,iserchk,ireschk,natsc,natmc,iat,jter, $ iressc(matsc),iresmc(matmc),isermc(matmc),isersc(matsc), $ izsc(matsc),izmc(matmc) real x(3),xyzmin(3),xyzmax(3),xf(3),xyzminf(3),xyzmaxf(3), $ angles(3),xx(3),xcom(3),xsc(3,matsc),xmc(3,matmc), $ qsc(matsc),qmc(matmc) real q,biso,bfac,occ, det, xdummy,utob, bfacmx, maxdist, anisou(6) real anischk,fac,batmc,batsc character atnam*4,restyp*4,idch*1,resno*4, filnam*80, insert*1, $ atnamsc(matsc)*4,restypsc(matsc)*4,insertsc(matsc)*1, $ atnammc(matmc)*4,restypmc(matmc)*4,insertmc(matmc)*1, $ idchsc(matsc)*1,idchmc(matmc)*1,segidmc(matmc)*1, $ resnomc(matmc)*4,resnosc(matsc)*4,segidsc(matsc)*1, $ atmtypmc(matmc)*8,atmtypsc(matsc)*8, $ altcodmc(matmc)*1,altcodsc(matsc)*1 character insold*1, atmtyp*8, altcod*1, FILTYP*3 character*15 codes(7) C integer jbfac, jqfac logical lcellnew, lcellpdb,newres,ltrnsf,lfrtrn,lorthg, lspnam, $ lxplor, lnodum, lhydro, loutxp, lnewch, linvrt, lutob, lnowat, $ lngnmi, lngnmo, lreort, lnewhd, com, lnohed, lnosid, lfrtrf 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 /'*'/ data anisou /6*0.0/ 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 Give molecular weight of each type of amino acid DATA seqid3/'ALA ','ARG ','ASN ','ASP ','CYS ','GLN ','GLU ', + 'GLY ','HIS ','ILE ','LEU ','LYS ','MET ','ORN ', + 'PHE ','PRO ','SER ','THR ','TRP ','TYR ','VAL '/ DATA seqid1/'A','R','N','D','C','Q','E', 1 'G','H','I','L','K','M','X', 1 'F','P','S','T','W','Y','V'/ DATA MWSEQID/ 71 , 157 , 114 , 114 , 102 , 128 , 128, 1 57 , 138 , 113 , 113 , 129 , 131 , 114, 1 147 , 97 , 87 , 101 , 186 , 163 , 99 / C C----------------------------------- C C========================================================================== C Initializations C========================================================================== call ccpfyp call ccp4h_init() CALL CCPRCS(6,'PDBSET','$Date: 2008/09/17 13:46:49 $') c c Initialize coordinate input & output c init library common blocks in rwbrook.f call xyzinit() iun = 0 iunout = 0 ifail = 0 c open files in rwbrook.f. Read appropriate header info (PDB or (CIF)) call xyzopen( 'XYZIN', 'INPUT',' ', iun, ifail) call getorm(iun, rfin, roin, cellpdb, lflag) C lflag = 0 - No CRYST1 or SCALEi header if (lflag .eq. 0) then lcellpdb = .false. lorthg = .false. end if C lflag > 0 - CRYST1 read; if no SCALEi ncode defaults to 1 if (lflag .gt. 0) then lcellpdb = .true. lorthg = .true. end if c--- C n=0 lseqnc = .false. kseqpd = 0 seqfil = 'SEQUENCE' nseqn = 0 nsqchn(1,1) = 1 loutxp = .false. lcellnew = .false. lnewhd = .false. lspnam = .false. lxplor = .false. lnodum = .false. lhydro = .false. lreort = .false. lnohed = .false. lnosid = .false. lnowat = .false. nreplc = 0 ncdin = -1 ncdout = -1 nchain = 0 mchain = 0 ltrnsf = .false. lfrtrn = .false. linvrt = .false. latren = .false. lnoise = .false. rnoise = 0.2 nnoise = 0 nnlchn = 0 nnlatm = 0 pxval = 0.0 lutob = .false. utob = ((4.0*atan(1.0))**2)*8.0 jbfac = 0 jbfac = 0 iressc(1) = -999 iresmc(1) = -999 restypmc(1) = '?' restypsc(1) = '?' insertmc(1) = '?' insertsc(1) = '?' idchmc(1) = '?' idchsc(1) = '?' jqfac = 0 qmin = 0.0 lqmin = .false. bfmax = 99.0 lbfmax = .false. lbfmin = .false. mwtot = 0.0 do 1, i=1,maxren lincr(i) = .false. 1 continue do 2, i = 1,3 xyzmin(i) = +10000000. xyzmax(i) = -xyzmin(i) xyzminf(i) = +1000. xyzmaxf(i) = -xyzminf(i) 2 continue xdummy = 1499. nselch = 0 nselat = 0 numren = 0 nelmnt = 0 FILTYP = ' ' SPGRPpdb = ' ' NAMSPG_CIFS = ' ' NAMPG = ' ' spgnam = ' ' lspgrp = 0 lspgrpPDB = 0 nsym = 0 lsymgn = 0 bfac = 20.0 bfacmx = 99.0 occ = 1.0 ncode = 1 ierr = 0 insold = ' ' do i = 1,4 do j = 1,4 do k = 1,96 rsym(i,j,k) = 0.0 rsympdbs(i,j,k) = 0.0 enddo tmp(i,j) = 0.0 uf(i,j) = 0.0 tmpT(i,j) = 0.0 trnsfm(i,j) = 0.0 enddo rsym(i,i,1) = 1.0 rsympdbs(i,i,1) = 1.0 tmp(i,i) = 1.0 trnsfm(i,i) = 1.0 enddo CALL RBSPGRP(SPGRPpdb) IF(SPGRPpdb.NE.' ') + call MSYMLB3(1,LSPGRPpdb,SPGRPpdb,NAMSPG_CIFS, + NAMPG,NSYMPpdbs,NSYMpdbs,RSYMpdbs) c Pass through TER lines iter = 0 C Set flag to indicate end_of_file jter = 0 C C========================================================================== C Read command input C========================================================================== 100 line=' ' ntok=maxtok C calls parser in parser.f, to parse command line(s) call parser( . key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.) C check for end-of-file lend ==.true. if(lend) go to 200 C check for blank line, if so try again if(ntok.eq.0) go to 100 C C........................................................................ C end of input tocken if (key .eq. 'END') then go to 200 C........................................................................ C optional REMARK. echo to output file. elseif (key .eq. 'REMA') then if (ntok .gt. 1) then k = min(iend(ntok), ibeg(2)+71) line = 'REMARK '//line(ibeg(2):k) C associate remark with INPUT unit number, so that mmdb will then copy C it from the input channel to the output channel. Sending it directly C to the output channel would cause it to be overwritten (NDS). call wremark(iun, line(1:lenstr(line))) else line = 'REMARK ' call wremark(iun, line(1:lenstr(line))) endif C........................................................................ C CELL option. generate CRYSTl and SCALEn from input 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 cellnew(4) = 90.0 cellnew(5) = 90.0 cellnew(6) = 90.0 do 110, i=1,6 call gttrea(i+1,cellnew(i),lflag,ntok,ityp,fvalue) if (lflag .ne. 0 .and. i .le. 3) go to 199 110 continue lcellnew = .true. C........................................................................ C ORTHOGONALISE option. define orthogonalisation code. C only effective if new CELL defined. 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........................................................................ C SPACEGROUP option. define spacegroup. should use with CELL elseif (key .eq. 'SPAC') then C SPACEGROUP, read spacegroup name for output spgnam = line(ibeg(2):iend(2)) call ccpupc(spgnam) C C#### Allow use of spacegroup number instead of name as in previous C#### versions. READ(SPGNAM,*,ERR=111) LSPGRP 111 call MSYMLB3(1,LSPGRP,SPGNAM,SPGNAM_CIFS, + NAMPG,NSYMPpdbs,NSYMpdbs,RSYMpdbs) C Check SPACEGROUP consistency if(lspgrp.ne.lspgrpPDB .and. lspgrpPDB.gt.0) then CALL CCPERR(2,' You are changing the PDB spacegroup') CALL CCPERR(2, + ' Any Symm ops on the REMARK 290 cards must be deleted') endif lspnam = .true. C........................................................................ C SYMGEN option. generate new chains using supplied symmetry operation(s) 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 write(6,'(a,3i5,4(/,4f8.3))') + ' numsgp,nsym,nsymp,rsym',numsgp,nsym,nsymp, + ((rsym(i,j,1),j=1,4),i=1,4) endif C........................................................................ elseif (key .eq. 'SYMM') then write (6, '(a/a)') $ ' *** SYMMETRY command withdrawn ***', $ ' *** Use SYMGEN command: read write-up ***' C........................................................................ C Chain option. naming of chains 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........................................................................ C BFACTOR option 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 or ZERO) C = 3 reset if .lt. given value (MINIMUM) C = 4 reset if greater than given value (MAXIMUM) C = 5 truncate to given range (RANGE) C = 6 Average Main and Side chain B factors ( cf BAVERAGE) jbfac = 1 n = 2 call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'ALWA') then jbfac = 1 n = 3 elseif (cvalue(2) .eq. 'DEFA' .or.cvalue(2) .eq. 'ZERO') 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 elseif (cvalue(2) .eq. 'AVER') then jbfac = 6 n = 2 go to 100 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........................................................................ C OCCUPANY option 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 or ZERO) C = 3 reset if .lt. given value (MINIMUM) C = 4 reset =0.0 if .lt. given value, else set = 1.0 (RESET) jqfac = 1 n = 2 call ccpupc(cvalue(2)) if (cvalue(2) .eq. 'ALWA') then jqfac = 1 n = 3 elseif (cvalue(2) .eq. 'DEFA' .or.cvalue(2) .eq. 'ZERO') then jqfac = 2 n = 3 elseif (cvalue(2) .eq. 'MINI') then jqfac = 3 n = 3 elseif (cvalue(2) .eq. 'RESE') then jqfac = 4 n = 3 endif call gttrea(n,occ,lflag,ntok,ityp,fvalue) if (lflag .gt. 0) go to 199 C........................................................................ C RENUMBER option 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........................................................................ C SELECT option 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........................................................................ C PICK option 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........................................................................ C ROTATE option 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........................................................................ C SHIFT option 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........................................................................ C TRANSFORM option elseif (key .eq. 'TRAN') then C TRANSFORM define 12 elements of transformation ltrnsf = .true. lfrtrf = .false. 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' .or. cvalue(n) .eq. 'FILE') $ then C Read transformation from O database file or other file if (ntok .gt. n) then filnam = line(ibeg(n+1):) else filnam = 'ODB' endif if (cvalue(n) .eq. 'ODB') then lflag = 0 else lflag = 1 endif call gtrtdb(lunit, filnam, lflag, trnsfm, istat) if (istat .ne. 0) then write (6, '(a)') $ ' **** Error in reading O datablock ****' go to 199 endif lfrtrf = .false. lfrtrn = .false. n = ntok else if (cvalue(n) .eq. 'FRAC') then C If keyword FRACTIONAL present, set flag that translation is c fractional lfrtrf = .true. lfrtrn = .true. n = n+1 else lfrtrf = .false. 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 c lrfsym = 0 orthogonal translation c = 1 fractiona translation if (lfrtrf) then lrfsym(nsym) = 1 else lrfsym(nsym) = 0 endif write (6,'(/a,i3,4(/4f8.3))') $ ' Symmetry Transformation',nsym, $ ((rfsym(i,j,nsym),j=1,4),i=1,4) endif C........................................................................ C XPLOR option 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........................................................................ C SEQUENCE option 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. C Initialize control array do n = 1, maxchn nsqchn(1, n) = 0 nsqchn(2, n) = -1 end do 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........................................................................ C OUTPUT option 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. elseif (cvalue(n) .eq. 'PDB') then filtyp = 'PDB' elseif (cvalue(n) .eq. 'CIF') then filtyp = 'CIF' endif n = n+1 go to 170 endif C........................................................................ C UTOB option elseif (key .eq. 'UTOB') then C UTOB option, convert Us to B lutob = .true. C........................................................................ C ELEMENT option 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........................................................................ C REORTHOGONALIZE option 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 No Output orthogonalization? if (ncdout .le. 0) then write (6, '(/a/)') $ '*** No output REORTHOGONALIZATION given ***' call ccperr(1,'*** REORTHOGONALIZATION TO ??? ***') endif C........................................................................ C REPLACE option 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 contain up to c 3 characters (right-justified into 3 characters), c but atomnames are 4 characters, left-justified c Atom names are LEFT-JUSTIFIED as returned by XYZATOM c Residue names are not tname = line(ibeg(n):iend(n)) if (lflrtp(nreplc) .eq. 1) then c Right-justify residue name into 3 characters call rjstfy(restp3, tname) oldrtp(nreplc) = restp3 else c Left-justify atom name call ljstfy(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 c Right-justify residue name into 3 characters call rjstfy(restp3, tname) newrtp(nreplc) = restp3 else c Don't left-justify new atom name 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........................................................................ C EXCLUDE option 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. elseif (cvalue(n) .eq. 'SIDE') then lnosid = .true. elseif (cvalue(n)(1:3).eq.'WAT'.or.cvalue(n)(1:3).eq.'HOH')then lnowat = .true. endif C........................................................................ C ATRENUMBER option elseif (key .eq. 'ATRE') then latren=.true. C........................................................................ C NOISE option elseif (key .eq. 'NOIS') then lnoise=.true. rnoise=0.2 C initialise the random number generator (call rmarin) C with variants on the time (ustime) call ustime(iseed) iseed2 = iseed + 12345 iseed = mod(iseed,31329) iseed2 = mod(iseed2,30082) call rmarin(iseed,iseed2) n=2 330 if(n .le. ntok) then C check if number (rnoise is only option at this stage) if(ityp(n) .eq. 2 ) then rnoise = fvalue(n) if (rnoise .gt. 0.5) then write(6,'(a)')'Noise factor too large' goto 199 endif else C character input call ccpupc(cvalue(n)) if (cvalue(n) .eq. 'BFAC') then C B-FACTOR section (one non-optional value) lbfmin = .true. if (n .lt. ntok) then call gttrea(n+1,bfmin,lflag,ntok,ityp,fvalue) if (lflag .ne. 0) then write(6,'(a)')'Require B-factor' go to 199 endif n = n+1 endif elseif (cvalue(n) .eq. 'CHAI') then C CHAIN selection (chain values to end of line) k = 0 do 331, 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 nolchn(k) = line(ibeg(i):iend(i)) 331 continue nnlchn=k n=ntok C PICK selection (to end of line) elseif (cvalue(n) .eq. 'PICK') then k = 0 do 332, i = n+1, ntok k = k+1 if (k .gt. maxpat) then write (6, '(a)') ' ** Too many PICKed atoms +SELECT chains**' ierr = ierr+1 go to 100 endif nolatm(k) = line(ibeg(i):iend(i)) 332 continue nnlatm=k n=ntok endif endif n=n+1 goto 330 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 iunout = 0 ifail = 0 call xyzopen('XYZOUT','OUTPUT',FILTYP,iunout, ifail) C lngnmi = .false. lngnmo = .false. C C---- Reorthogonalization option for pdb cell, reset Roin etc using ncdin, if set if (lcellpdb .and. ncdin.gt.-1) + call ortmat(cellpdb, ncdin, roin, rfin) C C---- If no cell in input file and new cell read, set cellpdb to cellnew if(.not.lcellpdb .and. lcellnew) then cellpdb(1) = cellnew(1) cellpdb(2) = cellnew(2) cellpdb(3) = cellnew(3) cellpdb(4) = cellnew(4) cellpdb(5) = cellnew(5) cellpdb(6) = cellnew(6) if(ncdin.gt.-1) then call ortmat(cellpdb, ncdin, roin, rfin) else call ortmat(cellpdb, ncode, roin, rfin) end if end if C C---- If no new cell read and there is a cell in input file, set cellnew to cellpdb C---- rout and rfout may well finish up same as roin.. if (.not.lcellnew .and. lcellpdb) then cellnew(1) = cellpdb(1) cellnew(2) = cellpdb(2) cellnew(3) = cellpdb(3) cellnew(4) = cellpdb(4) cellnew(5) = cellpdb(5) cellnew(6) = cellpdb(6) if(ncdout.gt.-1) then call ortmat(cellnew, ncdout, roout, rfout) else call ortmat(cellpdb, ncode, roout, rfout) end if end if C C---- Get orthogonalisation option for new cell, use reorth option if given, C---- else default ncode=1 if (lcellnew) then if(ncdout.gt.-1) then call ortmat(cellnew, ncdout, roout, rfout) else call ortmat(cellnew, ncode, roout, rfout) end if lorthg = .true. endif c if( lreort)then write (6, '(/a,a,a,a/)') $ ' Coordinates will be reorthogonalized from ', $ codes(ncdin),' to ', codes(ncdout) write(6,'(a,4(/,4f8.3,5x,4f8.3))') ' Input RF and RO ', + ((rfin(i,j),j=1,4),(roin(i,j),j=1,4),i=1,4) write(6,'(a,4(/,4f8.3,5x,4f8.3))') ' Output RF and RO ', + ((rfout(i,j),j=1,4),(roout(i,j),j=1,4),i=1,4) endif C if (lspnam .or. lcellnew .or. lreort) lnewhd = .true. C C---- Exclude headers: reset flag in library to prevent headers being written if (lnohed) call rwnohead C C---- New orthogonalization for header if (lnewhd) then C make PDB header C C set spacegroup name if provided if (lspnam) then call wbspgrp(spgnam) endif if (ncode.le.6) then call wbcell(iunout, cellnew, ncode) elseif (ncode.eq.7) then if (lspnam) then write(iunout,210) cellpdb,spgnam else write(iunout,210) cellpdb,SPGRPpdb endif write(iunout,215) (i,(rfout(i,j),j=1,3),i=1,3) 210 format('CRYST1',3f9.3,3f7.2,1x,a15) 215 format( 2('SCALE',i1,4x,3f10.6,5x,' 0.00000',/), $ 'SCALE',i1,4x,3f10.6,5x,' 0.00000') endif 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 .ne. 0) write (6, '(/a)') $ ' B-factors will be modified - all AnisoUs lost' C 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' else if (jbfac .eq. 6) then write (6, '(/a)') $ ' B-factors will be averaged over main and side chains' 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' else if (jqfac .eq. 4) then write (6, '(/a,f8.2,a/,a)') $ ' Occupancy will be reset to 0.0 if less than ',occ, $ ',',' or to 1.0 if greater 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 CC write(6,*)' before do 240 nsym',nsym 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), rfin) call matmln(4, rfsym(1,1,n), roout, amf) write (6,'(/a,i3,4(/4f8.3))') $ ' Symmetry Transformation',n, $ ((rfsym(i,j,n),j=1,4),i=1,4) 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 if (lrfsym(isym) .gt. 0) then c Multiply translation vector by fractional to orthogonal matrix ro call transfrm(rfsym(1,4,isym), roin) endif 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 C if have transformation, ltrnsf==.true. if (ltrnsf) then C C---- If translation was fractional, convert to orthogonal C if have fractional tranlation lfrtrn==.true. if (lfrtrn) then C must have orthogonalisation matrix, lorthg==.true. if (lorthg) then c Multiply translation vector by fractional to orthogonal matrix ro call transfrm(trnsfm(1,4), roin) else write (6, '(/a,a/)') $ ' **** No orthogonalization matrix available', $ ' for fractional translation ****' call ccperr(1,'** Missing orthogonalization **') endif endif C C if have inversion, linvrt=.true. 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,' )'/) write (6,'(a, F12.5,/,a,/)') $ ' Determinant of requested rotation = ',det, $ ' This should be 1.0, otherwise check precision of input.' c if (lorthg) then c Just in case we are going to transform anisotropic Bs, set up matrices C Multiply rf by Trnsformation matrix and normalise ready for Uf_trn call matmln(4, tmp, rfin, trnsfm) rf_tru(4,4) = 1.0 do i = 1,3 fac= sqrt(tmp(i,1)*tmp(i,1) + tmp(i,2)*tmp(i,2) + + tmp(i,3)*tmp(i,3)) rf_tru(i,1) = tmp(i,1)/fac rf_tru(i,2) = tmp(i,2)/fac rf_tru(i,3) = tmp(i,3)/fac rf_tru(i,4) = 0.0 rf_tru(4,i) = 0.0 end do call rbrinv(rf_tru,ro_tru) endif 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 (lnowat) then write (6,'(a)') $ ' All water atoms removed' endif if (lnosid) then write (6,'(a)') $ ' All side chain past CB and non-protein atoms removed' 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---- NOISE start C testing and output for noise if (lnoise) then write (6,'(/A,F10.4)') +'Random shift introduced with maximum value: ',rnoise if (nnlchn .gt. 0) write (6,'(/A,20(2X,A))') $ ' Perform for chains : ',(nolchn(i), i=1,nnlchn) if (lbfmin) write (6,'(/A,F10.4)') $ ' Perform for atoms with Bfactor .gt. ',bfmin if (nnlatm .gt. 0) write (6,'(/A,20(2X,A))') $ ' Perform for atom types : ',(nolatm(i), i=1,nnlatm) endif C---- ATREN start C Reflect that atoms will be renumbered sequentially on output if (latren) then write (6,'(/A,/A)') $'Atoms from input file will be output with sequential numbering', $'(i.e. no gaps in numbers)' endif C C========================================================================== C Do it C========================================================================== jres = -1000000 natom = 0 natmc = 0 natsc = 0 batmc = 0 batsc = 0 chnold = '?' nqmin = 0 C C------------------------------------------------------------- C Symmetry loop isym = 1 300 continue C------------------------------------------------------------- C Atom loop C note: if not selected go immediately to top of loop (400) 400 continue c if (isym .eq. 1) then c Copy header records only for first symmetry pass icopy = iunout else icopy = 0 endif c If exclude headers then this is necessary to suppress remarks if (lnohed) icopy = 0 c anisflg = .FALSE. C copy header information to output file if 1st symmetry call xyzadvance(iun, icopy, iter, *400, *500) C read PEB atom card from buffer call xyzatom(iun, iser, atnam, restyp, idch, ires, * resno, insert, altcod, segid, iz, atmtyp) call xyzcoord(iun, 'O', 'O', x(1), x(2), x(3), q, biso, anisou) C write(6,'(A,I5,A4,A4,A1,I4,3F8.3,2F6.2,/,6f10.5)') C $ ' input line',iser, atnam, restyp, idch, ires, x(1), x(2), x(3), C $ q, biso, anisou anischk = 0.00 anischk = anisou(2)*anisou(2) + anisou(3)*anisou(3) if(anischk .gt. 0.0001 .or. anisou(2).ne. 0.00 ) then anisflg = .true. C Cant handle resetting Bisos with Aniso Us if(jbfac .ne. 0) then write (6, '(/a/)') $ ' **** PDBSET cannot handle resetting Bfactors if'// $ ' anisotropic Us are present in the file ****' call ccperr(1, '*** AnisoU present: cannot reset B ***') endif Cejd - too difficult to reothogonalize or rotate AnisoUs ! if(lreort ) then write (6, '(/a/)') $ ' **** PDBSET cannot handle reorthogonalization if'// $ ' anisotropic Us are present in the file ****' call ccperr(1, $ '*** AnisoU present: cannot reothogonalize ***') endif C C end if c c chnnam is used for all selections etc c This is set to c (1) idch c (2) if idch blank & segid is set, then set chnnam = segid chnnam = idch if (idch .eq. ' ') chnnam = segid c c---- Omit hydrogens if required if (lhydro) then C if (atnam(1:1).eq.'H' .or. atnam(2:2).eq.'H') go to 400 if (atmtyp(2:2).eq.'H') go to 400 endif c c---- Omit waters if required if (lnowat) then if (restyp(1:3).eq.'WAT'.or.restyp(1:3).eq.'HOH') go to 400 endif c c---- Omit side chain and non-protein atoms if required if (lnosid) then if ( atnam(1:2).ne.'CA' .and. atnam(1:2).ne.'C ' + .and. atnam(1:2).ne.'N ' .and. atnam(1:2).ne.'O ' + .and. atnam(1:2).ne.'CB' ) 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') .or. $ resno(j:j) .eq. '-' .or. resno(j:j) .eq. '+') 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) chnnam = idch resno(j:j) = ' ' endif C C---- Is this a new chain? lnewch = .false. if (chnnam .ne. chnold) then lnewch = .true. chnold = chnnam endif C---- List of chains in input file (after unpacking residue numbers) cc PhilE I don't understand the test on segid here c if (segid .ne. ' ' .and. isym .eq. 1) then c So replace with this if ((lnewch .or. mchain .eq. 0) .and. isym .eq. 1) then if (mchain .gt. 0) then do 402, j = 1,mchain if (chnnam .eq. chains(j)) go to 403 402 continue endif mchain = mchain+1 if (mchain .le. maxchn) then chains(mchain) = chnnam endif 403 continue endif C C---- Select options if (nselch .gt. 0) then do 410, j = 1, nselch if (chnnam .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 (biso .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 biso = utob*biso endif C C---- Reset Bfactor & occupancy C Bfactor option, set B & occupancy 461 if (jbfac .eq. 1) then biso = bfac elseif (jbfac .eq. 2) then if (biso .eq. 0.0) biso = bfac elseif (jbfac .eq. 3) then biso = max(biso, bfac) elseif (jbfac .eq. 4) then biso = min(biso, bfac) elseif (jbfac .eq. 5) then biso = min(max(biso, bfac), bfacmx) elseif (jbfac .eq. 6) then C Does this atom belong to the same residue as previous one? cc call xyzatom(iun, iser, atnam, restyp, idch, ires, cc * resno, insert, altcod, segid, iz, atmtyp) cc call xyzcoord(iun, 'O', 'O', x(1), x(2), x(3), q, biso, anisou) c write(6,'(A,I5,A4,A4,A1,I4,3F8.3,2F6.2,/,6f10.5)') c $ ' input line',iser, atnam, restyp, idch, ires, x(1), x(2), x(3), c $ q, biso, anisou C c write(6,'(a,2(I4,f8.1),2i5,2(/,I5,3(2x,a)))') c $ ' natmc,batmc,natsc,batsc,iresmc(1),iressc(1)', c $ natmc,batmc,natsc,batsc,iresmc(1),iressc(1), c $ iressc(1), restypsc(1),insertsc(1),idchsc(1), c $ iresmc(1), restypmc(1),insertmc(1),idchmc(1) C C Is this part of the same residue? iressc(1) and iresmc(1) = -999 for new start 460 continue if( ( (iressc(1).ne.-999) .and. (jter.eq.1 .or. $ ires .ne.iressc(1) .or. restyp .ne. restypsc(1) $ .or. insert.ne.insertsc(1).or. idch .ne. idchsc(1)) ) .or. $ ( (iresmc(1).ne.-999) .and. (jter.eq.1 .or. $ ires .ne.iresmc(1) .or. restyp .ne. restypmc(1) $ .or. insert.ne.insertmc(1).or. idch .ne. idchmc(1)) ) ) then C C No - then If there is any coordinates stored, output them. c write(6,*)' about to write; natsc,natmc',natsc,natmc C Come here to complete last residue averaging etc. if(natsc.gt.0 .or. natmc.gt.0) then if(natmc.gt.0) then batmc = batmc/natmc anisou(1) = batmc do iat = 1,natmc if (latren) then iserout = natom + 1 else iserout = isermc(iat) endif call xyzatom(iunout, iserout, atnammc(iat), * restypmc(iat), idchmc(iat), iresmc(iat), * resnomc(iat), insertmc(iat), altcodmc(iat), * segidmc(iat), izmc(iat), atmtypmc(iat)) call xyzcoord $ (iunout, 'O', 'O', xmc(1,iat), xmc(2,iat), $ xmc(3,iat), qmc(iat), batmc, anisou) call xyzadvance(iunout, 0, 0, *400, *400) natom = natom+1 end do natmc = 0 batmc = 0.0 iresmc(1)=-999 end if C Write out Side chain now if(natsc.gt.0) then batsc = batsc/natsc anisou(1) = batsc do iat = 1,natsc if (latren) then iserout = natom + 1 else iserout = isersc(iat) endif call xyzatom(iunout, isersc(iat), atnamsc(iat), * restypsc(iat), idchsc(iat), iressc(iat), * resnosc(iat), insertsc(iat), altcodsc(iat), * segidsc(iat), izsc(iat), atmtypsc(iat)) call xyzcoord $ (iunout, 'O', 'O', xsc(1,iat), xsc(2,iat), $ xsc(3,iat), qsc(iat), batsc, anisou) call xyzadvance(iunout, 0, 0, *400, *400) natom = natom+1 end do natsc = 0 batsc = 0.0 iressc(1)=-999 end if end if end if if(jter.gt.0) go to 501 c write(6,*)' about to store and update ; natsc,natmc',natsc,natmc if ( atnam(1:2).ne.'CA' .and. atnam(1:2).ne.'C ' + .and. atnam(1:2).ne.'N ' .and. atnam(1:2).ne.'O ') then natsc = natsc + 1 batsc = batsc + biso xsc(1,natsc) = x(1) xsc(2,natsc) = x(2) xsc(3,natsc) = x(3) isersc(natsc) = iser iressc(natsc) = ires izsc(natsc) = iz qsc(natsc) = q atnamsc(natsc) = atnam restypsc(natsc) = restyp idchsc(natsc) = idch resnosc(natsc) = resno insertsc(natsc) = insert altcodsc(natsc) = altcod segidsc(natsc) = segid(1:1) atmtypsc(natsc) = atmtyp else natmc = natmc + 1 batmc = batmc + biso xmc(1,natmc) = x(1) xmc(2,natmc) = x(2) xmc(3,natmc) = x(3) isermc(natmc) = iser iresmc(natmc) = ires izmc(natmc) = iz qmc(natmc) = q atnammc(natmc) = atnam restypmc(natmc) = restyp idchmc(natmc) = idch resnomc(natmc) = resno insertmc(natmc) = insert altcodmc(natmc) = altcod segidmc(natmc) = segid(1:1) atmtypmc(natmc) = atmtyp end if c write(6,*)' natmc natsc updated',natmc,natsc go to 400 end if C No more B factor fudging.. 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) elseif (jqfac .eq. 4) then if (q .lt. occ) then q = 0.0 else q = 1.0 endif endif C C---- Rename chains C Override chain-ID if CHAIN command given if (nchain .gt. 0) then do 430, j = 1, nchain if (chnnam .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 chnnam = newchn(j, isym) endif endif endif C C---- Residue number C if (jres .ne. ires .or. insert.ne.insold) then if (jres .ne. ires .or. insert.ne.insold .or. lnewch) 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. chnnam) 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) chnnam = chnrnw(j) go to 445 endif endif 440 continue endif 445 continue c C---- All this applies to coordinates c---- Reorthogonalization option if (lreort) then c coordinate -> fractional by input orthogonalization call transfrm(x, rfin) c orthogonalize with output orthogonalization call transfrm(x, roout) endif C C---- Coordinate transformations if (nsym .gt. 0) then C Symmetry generation if(.not.anisflg) then call transfrm(x, rfsym(1,1,isym)) elseif (anisflg) then call transfrm(x, rfsym(1,1,isym)) c Uf_tr = Sym(j) [Uf] SymT(j) C Transform orthog to frac c write(6,'(a,6f10.5)')' orthog1',anisou call cvanisou(anisou,1) uf(1,1)= anisou(1) uf(2,2)= anisou(2) uf(3,3)= anisou(3) uf(1,2)= anisou(4) uf(2,1)= anisou(4) uf(1,3)= anisou(5) uf(3,1)= anisou(5) uf(2,3)= anisou(6) uf(3,2)= anisou(6) c call matmln(4, tmp, rsym(1,1,isym), uf) do i = 1,3 do j = 1,3 tmpt(i,j) = rsym(j,i,isym) enddo enddo call matmln(4, uf, tmp, tmpt) C Transform frac to orthog anisou(1) = uf(1,1) anisou(2) = uf(2,2) anisou(3) = uf(3,3) anisou(4) = uf(1,2) anisou(5) = uf(1,3) anisou(6) = uf(2,3) c write(6,'(a,6f10.5)')' frac1',anisou call cvanisou(anisou,0) c write(6,'(a,6f10.5)')' orthog2',anisou endif endif C---- transformation section if (ltrnsf) then if(.not.anisflg) then call transfrm(x, trnsfm) elseif (anisflg) then call transfrm(x, trnsfm) c Transform anisotropic B factors Cejd - I THINK this might be right! C [Uf_tr_ij] = [RFutr] [Uo_ij] {[RFutr]}T C Then [Uo_tr_ij] = [RFu]-1 [Uf_tr_ij] {[RFu]-1}T C Convert othog anisou to fractional call cvanisou(anisou,1) uf(1,1)= anisou(1) uf(2,2)= anisou(2) uf(3,3)= anisou(3) uf(1,2)= anisou(4) uf(2,1)= anisou(4) uf(1,3)= anisou(5) uf(3,1)= anisou(5) uf(2,3)= anisou(6) uf(3,2)= anisou(6) call matmln(4, tmp, rf_tru, uf) do i = 1,3 do j = 1,3 tmpt(i,j) = rf_tru(j,i) enddo enddo call matmln(4, uf, tmp, tmpt) C Transform frac to orthog anisou(1) = uf(1,1) anisou(2) = uf(2,2) anisou(3) = uf(3,3) anisou(4) = uf(2,1) anisou(5) = uf(1,3) anisou(6) = uf(2,3) call cvanisou(anisou,0) endif endif C C---- start NOISE if (lnoise) then C allowed if selected chain if (nnlchn .gt. 0) then do 333, j = 1, nnlchn if (chnnam .eq. nolchn(j)) go to 334 333 continue go to 339 endif 334 continue C allowed if selected atom type if (nnlatm .gt. 0) then do 335, j = 1, nnlatm if (eqlstr(atnam, nolatm(j))) go to 336 335 continue go to 339 endif 336 continue C allowed if B-factor greater than bfmin if (lbfmin .AND. biso .lt. bfmin) goto 339 C track number of events, nnoise nnoise=nnoise+1 C also track shifts, pxval call noise(x,rnoise,pxval) endif 339 continue C C c---- Residue replace if (nreplc .gt. 0) then do 455, i = 1, nreplc if (lflrtp(i) .eq. 1) then if (eqlstr(restyp, oldrtp(i))) then restyp = newrtp(i) endif elseif (lflrtp(i) .eq. 2) then if (resatp(i) .eq. ' ' $ .or. (eqlstr(resatp(i), restyp))) then if (eqlstr(atnam, oldrtp(i))) then c Atom name should be left-justified for s/r xyzatom: c leading space is picked up from atmtyp if necessary call ljstfy(atnam, newrtp(i)) atmtyp = newrtp(i)(1:2) 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 set starting pointer C Note: mchain is new current chain nsqchn(1, mchain) = nseqn nsqchn(2, mchain) = nseqn else C Not a new chain, so keep updating the end pointer nsqchn(2, mchain) = nseqn endif seqnce(nseqn) = restyp endif endif C---- This used to convert Chain ID to uppercase, but why? C PDB has case-sensitive chain IDs e.g. 1fnt, and MMDB can handle it. idch = chnnam(1:1) if (segid .eq. ' ') then segid = chnnam endif C Get limits here do 450, i = 1,3 xyzmin(i) = min(xyzmin(i), x(i)) xyzmax(i) = max(xyzmax(i), x(i)) xf(i) = rfout(i,1)*x(1) + rfout(i,2)*x(2) + + rfout(i,3)*x(3) + rfout(i,4) xyzminf(i) = min(xyzminf(i), xf(i)) xyzmaxf(i) = max(xyzmaxf(i), xf(i)) 450 continue C---- Write out natom = natom+1 if(.not.anisflg) then anisou(1) = biso anisou(2) = 0.0 anisou(3) = 0.0 anisou(4) = 0.0 anisou(5) = 0.0 anisou(6) = 0.0 end if if (latren) then iserout = natom else iserout = iser endif call xyzatom(iunout, iserout, atnam, restyp, idch, ires, * resno, insert, altcod, segid, iz, atmtyp) call xyzcoord $ (iunout, 'O', 'O', x(1), x(2), x(3), q, biso, anisou) c write(6,'(A,I5,A4,A4,A1,I4,3F8.3,2F6.2,/,6F10.5)') c $ ' output line',iser, atnam, restyp, idch, ires, x(1), x(2), x(3), c $ q, biso, anisou call xyzadvance(iunout, 0, 0, *400, *400) C write coordinates to scratch file for another pass if (com) then call qwrite(1,x,3) xcom (1) = xcom (1) + x (1) xcom (2) = xcom (2) + x (2) xcom (3) = xcom (3) + x (3) end if go to 400 C---- End atom loop 500 continue jter = 1 if(jbfac.eq.6 .and. (natsc.gt.0 .or. natmc.gt.0)) go to 461 501 continue if (isym .lt. nsym) then C---- Symmetry loop isym = isym+1 call xyzrewd(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 C final noise rms value output if (lnoise) then C get rms value of shifts (3 because of XYZ) if (nnoise .gt. 0) then pxval=sqrt(pxval/(3.0*nnoise)) write(6,6509) rnoise,nnoise,pxval 6509 format(//'Random shift'/7x,'maximum: ',f6.4/ + 7x,'atoms: ',2x,i6/7x,'RMS: ',4x,f6.4) else write(6,'(//A)')'WARNING*** Random + shift selected, but no atoms affected' endif endif do 510, i = 1,3 x(i) = 0.5 * (xyzmin(i) + xyzmax(i)) xf(i) = 0.5 * (xyzminf(i) + xyzmaxf(i)) 510 continue if (natom .gt. 0) then C Output: C rf(4,4) orthogonal to fractional matrix C ro(4,4) fractional to orthogonal matrix write (6,6500) $ (xyzmin(i),xyzmax(i),x(i),xyzmax(i)-xyzmin(i), i=1,3), $ (xyzminf(i),xyzmaxf(i),xf(i),xyzmaxf(i)-xyzminf(i), i=1,3) 6500 format(//' Orthogonal 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/, $ //' Fractional 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 qread(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 mwtot=0 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 C Get molecular weight of this chain mwchn(n) =0 do i=nsqchn(1,n),nsqchn(2,n) do k = 1,maxaa if(seqnce(i).eq.seqid3(k)) mwchn(n) = mwchn(n) + + mwseqid(k) end do end do mwtot = mwtot + mwchn(n) write(6,'(a,I3,a,3x,A,I8,A,I8)') + ' ChainID ',n, chains(n), + ' Molecular Weight ',mwchn(n), + ' Cumulative Molecular Weight',mwtot C Only write REMARK card if there are residues in this chain if (nsqchn(2,n).ge.nsqchn(1,n)) 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,2x,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)(1:3), seqnc1(i)) 530 continue do 535, n = 1, mchain C Get molecular weight of this chain mwchn(n) =0 do i=nsqchn(1,n),nsqchn(2,n) do k = 1,maxaa if(seqnc1(i).eq.seqid1(k)) mwchn(n) = mwchn(n) + + mwseqid(k) end do end do mwtot = mwtot + mwchn(n) write(6,'(a,I3,3x,a,3x,A,I8,A,I8)') + ' ChainID ',n, chains(n), + ' Molecular Weight ',mwchn(n), + ' Cumulative Molecular Weight',mwtot C Only write card if there are residues in this chain if (nsqchn(2,n).ge.nsqchn(1,n)) 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 mwtot = 0 do i=1,nseqn do k = 1,maxaa if(seqnce(i).eq.seqid3(k)) mwtot = mwtot + + mwseqid(k) end do end do write(6,*) ' Total Molecular Weight from sequence',mwtot write (9, '(1000(10(A4,2X)/))') (seqnce(i), i = 1,nseqn) endif endif C call xyzclose(iun) call xyzclose(iunout) if (lseqnc) then write(6,*) ' Molecular Weight from sequence ',mwtot,' Da' endif 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 subroutine getorm(iun, rf, ro, cell, lflag) C =========================================== C C Get cell & orthogonlization matrices from PDB file, if available 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 = 3 if cell or sCALEi found, matrices taken from SCALE lines C Options 1 & 2 not now given C C integer iun, lflag real rf(4,4), ro(4,4), cell(6) C integer ncode real vol C C lflag = 0 C call rbcell(cell, vol) C if (cell(1) .gt. 0.000001) then ro(1,1) = 0.0 call rbrorf2(ro, rf, ncode) lflag = 3 else C No cell available lflag = 0 endif return end 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.001/ 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 subroutine ljstfy(a,b) c ====================== c c Left justify string b into a c character*(*) a,b integer j,lenstr external lenstr c do 1, j = 1,lenstr(b) if (b(j:j) .ne. ' ') go to 3 1 continue j = 1 3 a = b(j:) 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, lflag, trnsfm, istat) c ============================================= c c Read an O datablock file or other matrix 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 lflag = 0 O database format c = +1 file containg 3x3 R, t (ie 12 numbers) 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, lflag, istat real trnsfm(4,4) character file*(*) c integer nitems, i, j character blknam*80, type*1, format*80 c istat = 0 if (lflag .eq. 0) then 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) elseif (lflag .eq. +1) then call ccpdpn(lunit, file, 'READONLY', 'F', 0,0) read (lunit, *) ((trnsfm(i,j), j=1,3), i=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) endif 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 C C NOISE routine subroutine noise(x,rmax,pxval) C Inputs: x(3) coordinate C rmax maximum shift C pxval square of shifts C Outputs: x(3) shifted coordinate C pxval updated real x(3),rmax,pxval real rvec(3),xm integer i external ranmar C call ranmar(rvec,3) C do 10 i=1,3 xm=rmax*(rvec(i)+rvec(i)-1.0) pxval=pxval+xm*xm 10 x(i)=x(i)+xm C return end