C************************************************************* C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C************************************************************* c c SORTWATER c c Assign waters to appropriate, possibly non-crystallographically c related, molecules c c Based on watersort c c Phil Evans, MRC LMB, January 1995 c c Input commands: c c--> CHAINS c Define all chain IDs of "protein" (ie non-water) chains c c--> WCHAINS c Define chain names for water chains to correspond to "protein" c chains in output file (irrespective of input water chain names) c There must be the same number of water chains defined as "protein" c chains, but the same water chain may be assigned to more than one c protein chain, provided that they are not related by c non-crystallographic symmetry c c--> SYMMETRY ||symmetry c Define crystallographic symmetry c c--> WATER [] c Residue name for waters [default HOH], and atom name [default O] c c--> CARBON ["Yes"|"No"] c No store only non-carbon non-water atoms for contact checking (.true.) c [default] c Yes store all atoms (.false.) c c--> DISTANCE [] [] c maximum distance between putative NCS-related waters to accept c [default 2.0] c maximum distance from non-water atom to accept as belonging to chain c [default 6.0] c c--> NCS | ODB c | MATRIX - c c | IDENTITY c | SAME c c Define NCS operator to transform chain with ID "Chain1" to "Chain2" c c Operators may be given as the filename of an O data block, or as c 12 numbers following the keyword MATRIX (note the ODB file contains c the transposed matrix) c The keyword SAME defines the transformation from "Chain3" to "Chain4" c as being the same as that for "Chain1" to "Chain2". This may be put c at the end of a line defining an operator c Implied operators will be generated automatically (eg B->A from A->B, c and A->C from A->B & B->C) c c--> VERBOSE c Set verbose printing flag c c program sortwater c ================= c INTEGER IXYZIN c c call ccpfyp call ccprcs(6,'SORTWATER','$Date: 2002/06/14 12:11:35 $') call XYZINIT c c Read controls call keyin c c Read coordinate file & store call redmol(IXYZIN) c c Get closest atoms to each water in each chain call closest c c Match them call matchw c c Write them out on end of non-water atoms call wrimol(IXYZIN) c call ccperr(0,'*** Normal termination ***') end c c c subroutine keyin c ================ c c Read controls c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include waters.fh c# waters.fh c Character common c watres*4 water residue name c chnwat water chain names for each water atom c wchain possible water chain names for output c watnam atomname for water c uwchain unique water chain names (nqwchn of them) c ie water chain jch has naem uwchain(jch) c common /waterc/ watres, chnwat(maxwat), wchain(maxchn), $ uwchain(maxchn), watnam character watres*4, chnwat*1, wchain*1, uwchain*1, watnam*4 c Numbers etc c xwater orthogonal coordinates c fwater fractional coordinates c bwater Bfactor c qwater occupancy c distat(,ich) distance to nearest atom in chain ich c ctfat(3,,ich) cell translations to nearest atom in chain ich c xyzcrt(3,,ich) orthogonal coordinates of water moved by crystal c symmetry for protein chain ich c (ie lsymat(,ich), ctfat(3,,ich)) c nearat(,ich) nearest atom in chain ich c lsymat(,ich) symmetry number for nearest atom in chain ich c nearch nearest chain c nrswat residue number (input) c nwater number of waters c nwchain number of water chains c iwchain assigned chain number for water c mltplw multiplicity of this water c mltmax maximum multiplicity c kmswat "master" (ie 1st) water number related to this one c numwat input water number for this allocated water number c kwchain unique water chain number corresponding to protein chain c (allowing for same water chain to match multiple c protein chains) c nqwchn number of unique water chains c common /waters/ xwater(3,maxwat), fwater(3,maxwat), $ bwater(7,maxwat), qwater(maxwat),distat(maxwat,maxchn), $ ctfat(3,maxwat,maxchn), xyzcrt(3,maxwat,maxchn), $ nearat(maxwat,maxchn), lsymat(maxwat,maxchn), nearch(maxwat), $ nrswat(maxwat), nwater, nwchain, iwchain(maxwat), $ mltplw(maxwat), mltmax, kmswat(maxwat), $ numwat(maxwat,maxchn), kwchain(maxchn), nqwchn real xwater, fwater, bwater, qwater, distat, ctfat, xyzcrt integer nwater, nearch, nrswat, nwchain, nearat, lsymat, iwchain, $ mltplw, mltmax, kmswat, numwat, kwchain, nqwchn save /waterc/, /waters/ c## C&&*&& end_include waters.fh C&&*&& include flagcm.fh c# flagcm.fh c c Control flags etc c c dislim maximum distance check c dismol maximum allowed distance of water from "protein" molecule c dissim distance threshold for considering "equivalent" waters c disim2 dissim**2 c lonly .true. to check distances to 'O' & 'N' atoms only c verbose .true. for debug printing c common /flagcm/ dislim, dismol, dissim, disim2, lonly, verbose real dislim, dismol, dissim, disim2 logical lonly, verbose save /flagcm/ c## C&&*&& end_include flagcm.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh C&&*&& include symmcm.fh c# symmcm.fh c c Crystallographic symmetry c integer maxsym parameter (maxsym=96) c common /symmcm/ rsym(4,4,maxsym), rfsym(4,4,maxsym), $ nsym, nsymp, numsgp integer nsym, nsymp, numsgp real rsym, rfsym common /symmcc/ spgnam, pgname character spgnam*10, pgname*10 save /symmcm/, /symmcc/ c## C&&*&& end_include symmcm.fh C&&*&& include ncscm.fh c# ncscm.fh c c NCS operators c c rtncs(4,4,k) k'th NCS operator: operator 1 is always identity c numncs number of NCS operators c ncsopr(i,j) pointer to operator to move chain i on chain j c = 0 if none c common /ncscm/ rtncs(4,4,maxncs), numncs, ncsopr(maxchn,maxchn) real rtncs integer numncs, ncsopr save /ncscm/ c## C&&*&& end_include ncscm.fh c c Things for PARSER ---- integer maxtok parameter (maxtok=100) character key*4,cvalue(maxtok)*4,line*400 integer ibeg(maxtok),iend(maxtok),ityp(maxtok), . idec(maxtok),ntok real fvalue(maxtok) logical lend c c Locals integer ierr, i,j,k,l,m, ikey, lunout, idum, istat logical lsymin c Functions logical ccponl external ccponl c c Commands integer nkeys parameter (nkeys=8) character*16 keys(nkeys) c data keys/'CHAINS','WATER','CARBON','WCHAINS','SYMMETRY', $ 'DISTANCE','NCS','VERBOSE'/ c c Defaults lunout = 6 lonly = .true. lsymin = .false. verbose = .false. watres = 'HOH' watnam = 'O' nchain = 0 nwchain = 0 dismol = 6.0 dissim = 2.0 numncs = 1 call setim4(rtncs(1,1,1)) do 2, j = 1,maxchn do 1, i = 1,maxchn ncsopr(i,j) = 0 1 continue 2 continue c c c ---- Loop to read control input c ierr = 0 10 line=' ' ntok=maxtok call parser( . key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.) if(lend) go to 1000 if(ntok.eq.0) go to 10 call ccpupc(key) call chkkey(key,keys,nkeys,ikey) if(ikey.eq.0) then write(lunout,fmt='(a)') ' ** Unrecognized keyword **' if(ccponl(idum)) then goto 10 else ierr=ierr+1 endif else if(ikey.lt.0) then write(lunout,fmt='(a)') ' ** Ambiguous keyword **' if(ccponl(idum)) then goto 10 else ierr=ierr+1 endif endif c c go to (100,200,300,400,500,600,700,800), ikey c------------------------------------------------------------------------ c c CHAINS c 100 nchain = ntok-1 do 101, i = 1,nchain pchain(i) = cvalue(i+1)(1:1) 101 continue go to 10 c------------------------------------------------------------------------ c c WATER [] c 200 if (ntok .gt. 1) watres = cvalue(2) if (ntok .gt. 2) watnam = cvalue(3) go to 10 c------------------------------------------------------------------------ c c CARBON yes | no c c no store only non-carbon non-water atoms for contact checking (.true.) c [default] c yes store all atoms (.false.) c 300 if (ntok .gt. 1) then call ccpupc(cvalue(2)) if (cvalue(2)(1:1) .eq. 'N') then lonly = .true. elseif (cvalue(2)(1:1) .eq. 'Y') then lonly = .false. else go to 990 endif else lonly = .true. endif go to 10 c------------------------------------------------------------------------ c c WCHAINS c 400 nwchain = ntok-1 do 401, i = 1,nwchain wchain(i) = cvalue(i+1)(1:1) 401 continue go to 10 c------------------------------------------------------------------------ c c SYMMETRY ||symmetry c 500 j = 2 call rdsymm(j,line,ibeg,iend,ityp,fvalue,ntok, . spgnam,numsgp,pgname,nsym,nsymp,rsym) lsymin = .true. go to 10 c------------------------------------------------------------------------ c c DISTANCE c 600 call gttrea(2,dissim,istat,ntok,ityp,fvalue) call gttrea(3,dismol,istat,ntok,ityp,fvalue) go to 10 c c------------------------------------------------------------------------ c c NCS | ODB c | MATRIX - c c | SAME c | IDENTITY c c Define NCS operator to transform chain with ID "Chain1" to "Chain2" c c Operators may be given as the filename of an O data block, or as c 12 numbers following the keyword MATRIX (note the ODB file contains c the transposed matrix) c The keyword SAME defines the transformation from "Chain3" to "Chain4" c as being the same as that for "Chain1" to "Chain2" c Implied operators will be generated automatically (eg B->A from A->B, c and A->C from A->B & B->C) c 700 call setncs(line, ntok, ibeg, iend, ityp, $ fvalue, cvalue, istat) if (istat .ne. 0) then go to 990 endif go to 10 c------------------------------------------------------------------------ c c VERBOSE c Set verbose printing flag c 800 verbose = .true. go to 10 c------------------------------------------------------------------------ c------------------------------------------------------------------------ 1000 continue write (lunout, '(/,a,30(1x,a))') $ ' Non-water chain names: ',(pchain(i), i=1,nchain) write (lunout, '(/,a,30(1x,a))') $ ' Water chain names: ',(wchain(i), i=1,nwchain) write (lunout,'(/,a,a,/)') $ ' Water residue name: ',watres c if (nchain .ne. nwchain) then write (lunout, '(/,a,/)') $ ' **** Number of water & non-water chains is different ****' ierr = ierr+1 endif c if (lsymin) then c Inverse crystallographic symmetry do 1010, i = 1,nsym call rbrinv(rsym(1,1,i), rfsym(1,1,i)) 1010 continue else write (lunout, '(/,a,/)') ' **** No symmetry given ****' ierr = ierr+1 endif c if (ierr .gt. 0) then call ccperr(1,'*** SORTWATER: input error ***') endif c if (numncs .gt. 1) then c Complete NCS operators call ncscpl c if (verbose) then do 1100, m = 1, nchain do 1110, l = 1, nchain if (ncsopr(l,m) .gt. 1) then k = ncsopr(l,m) write (6, '(/,a,a,a,a,/)') $ ' NCS operator from chain ',pchain(l), $ ' to chain ',pchain(m) write (6, 6110) $ ((rtncs(i,j,k), j=1,4), i=1,3) 6110 format(3(3x,'(',3f12.5,' )',5x,'(',f8.3,' )',/)) endif 1110 continue 1100 continue endif endif disim2 = dissim**2 c Number water chains to allow for possibility of the same water chain c belonging to more than one protein chain nqwchn = 1 kwchain(1) = 1 uwchain(1) = wchain(1) if (nchain .gt. 1) then do 1200, l = 2, nchain do 1210, k = 1, l-1 c Has this chain been defined already? if (wchain(k) .eq. wchain(l)) then kwchain(l) = kwchain(k) c Non-unique water chain, check that there is no NCS relation c between these chains if (ncsopr(k,l) .gt. 0) then line = ' **** You may NOT have the same water'// $ 'chain ID for NCS-related "protein" chains ****' call ccperr(2, line) call ccperr(1,'*** Illegal water chain ID ***') endif go to 1200 endif 1210 continue c no, unique nqwchn = nqwchn+1 kwchain(l) = nqwchn uwchain(nqwchn) = wchain(l) 1200 continue endif return c 990 write (lunout,'(A)') ' ** Missing or illegal value **' 999 ierr = ierr+1 go to 10 c end c c c subroutine ncscpl c ================= c c Fill in NCS operators implied by those given explicitly c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh C&&*&& include ncscm.fh c# ncscm.fh c c NCS operators c c rtncs(4,4,k) k'th NCS operator: operator 1 is always identity c numncs number of NCS operators c ncsopr(i,j) pointer to operator to move chain i on chain j c = 0 if none c common /ncscm/ rtncs(4,4,maxncs), numncs, ncsopr(maxchn,maxchn) real rtncs integer numncs, ncsopr save /ncscm/ c## C&&*&& end_include ncscm.fh c integer k,l,m, numlast c numlast = numncs - 1 c Iterate from here as long as no new operators have been added c "do while (numlast .lt. numncs)" 1 if (numlast .lt. numncs) then numlast = numncs c do 10, m = 1, nchain do 20, l = 1, nchain c c Should the operator l->m be created? if (ncsopr(l,m) .eq. 0) then c try: c 1) identity, always operator 1 (preset) if (l .eq. m) then ncsopr(l,m) = 1 c 2) inverse elseif (ncsopr(m,l) .gt. 0) then numncs = numncs+1 if (numncs .gt. maxncs) go to 990 call rbrinv( $ rtncs(1,1,ncsopr(m,l)),rtncs(1,1,numncs)) ncsopr(l,m) = numncs else c 3) chain from somewhere else c We want to find a tranformation from chain l to chain m c Search each path from l->k and check if k->m exists do 30, k = 1,nchain if (ncsopr(l,k) .gt. 0) then if (ncsopr(k,m) .gt. 0) then c Yes, l->k & k->m exist, so [l->m] = [k->m][l->k] numncs = numncs+1 if (numncs .gt. maxncs) go to 990 call matmul4(rtncs(1,1,numncs), $ rtncs(1,1,ncsopr(k,m)), $ rtncs(1,1,ncsopr(l,k))) ncsopr(l,m) = numncs go to 50 endif endif 30 continue 50 continue c endif endif 20 continue 10 continue c end of "do while" go to 1 endif return c 990 write (6, '(/,a,i8,/)') $ ' **** NCS operator list full: ', maxncs call ccperr(1,'*** Too many NCS operators ***') end c c c subroutine setncs(line, ntok, ibeg, iend, ityp, $ fvalue, cvalue, istat) c ================================================================== c c Read NCS operator and store c NCS | ODB c | MATRIX - c c | IDENTITY c | SAME c c On entry: c numncs (/ncscm/) current number of NCS operators defined c line, ntok, ibeg, iend, ityp, fvalue, cvalue from Parser c c On exit: c numncs (/ncscm/) updated number of NCS operators c istat = 0 OK c = -1 syntax error c = -2 operator list full c = -3 illegal matrix c = -4 trying to set identity (ie chain A->A) c = -5 SAME command to undefined matrix c integer ntok, ityp(ntok), istat, ibeg(ntok), iend(ntok) real fvalue(ntok) character line*(*), cvalue(ntok)*(*) c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include ncscm.fh c# ncscm.fh c c NCS operators c c rtncs(4,4,k) k'th NCS operator: operator 1 is always identity c numncs number of NCS operators c ncsopr(i,j) pointer to operator to move chain i on chain j c = 0 if none c common /ncscm/ rtncs(4,4,maxncs), numncs, ncsopr(maxchn,maxchn) real rtncs integer numncs, ncsopr save /ncscm/ c## C&&*&& end_include ncscm.fh c c Locals integer i, j, k, ich1, ich2, ich3, ich4, n character errlin*80, filename*80, ch1*1, ch2*1, ch3*1, ch4*1, + key*4 c c integer lunit c Unit number for reading O files (assigned automatically) data lunit/0/ c c istat = 0 if (ntok .lt. 4) then go to 999 endif c Slot to store transformation numncs = numncs+1 if (numncs .gt. maxncs) then write (errlin, '(a,i8,a)') ' *** More than ', maxncs, $ ' NCS operators ***' call ccperr(2,errlin) istat = -2 return endif c c Chain names ch1 = cvalue(2)(1:1) ch2 = cvalue(3)(1:1) c Chain numbers call chnnum(ch1, ich1) if (ich1 .le. 0) then errlin = ' *** Unrecognized chain: '//ch1 call ccperr(2,errlin) go to 999 endif call chnnum(ch2, ich2) if (ich2 .le. 0) then errlin = ' *** Unrecognized chain: '//ch2 go to 999 endif if (ich1 .eq. ich2) then errlin = ' **** ILLEGAL attempt to set NCS for chain '// $ ch1//' on to itself ****' call ccperr(2,errlin) istat = -4 return endif c n = 4 key = cvalue(n) call ccpupc(key) if (key .eq. 'ODB') then c Read matrix from O data block file filename = line(ibeg(5):iend(5)) call gtrtdb(lunit, filename, rtncs(1,1,numncs), istat) if (istat .ne. 0) then go to 999 endif ncsopr(ich1,ich2) = numncs call chkmt4(rtncs(1,1,numncs), istat) if (istat .ne. 0) then go to 998 endif n = n+2 elseif (key .eq. 'MATR') then c Raed 12 numbers k = 4 do 100, j = 1,3 do 110, i = 1,3 k = k+1 call gttrea(k, rtncs(i,j,numncs), $ istat,ntok,ityp,fvalue) if (istat .ne. 0) go to 999 110 continue rtncs(i,4,numncs) = 0.0 100 continue do 120, i = 1,3 k = k+1 call gttrea(k, rtncs(i,4,numncs), $ istat,ntok,ityp,fvalue) if (istat .ne. 0) go to 999 120 continue rtncs(4,4,numncs) = 1.0 ncsopr(ich1,ich2) = numncs call chkmt4(rtncs(1,1,numncs), istat) if (istat .ne. 0) then go to 998 endif n = n+13 elseif (key .eq. 'IDEN') then call setim4(rtncs(1,1,numncs)) ncsopr(ich1,ich2) = numncs n = n+1 endif c if (ntok .ge. n) then key = cvalue(n) call ccpupc(key) if (key .eq. 'SAME') then if (ntok .ne. n+2) then go to 999 endif c Chain names ch3 = cvalue(n+1)(1:1) ch4 = cvalue(n+2)(1:1) c Chain numbers call chnnum(ch3, ich3) if (ich3 .le. 0) then errlin = ' *** Unrecognized chain: '//ch3 call ccperr(2,errlin) go to 999 endif call chnnum(ch4, ich4) if (ich4 .le. 0) then errlin = ' *** Unrecognized chain: '//ch4 call ccperr(2,errlin) go to 999 endif c c Copy matrix call cpymat(ncsopr(ich1,ich2), numncs, istat) if (istat .ne. 0) then errlin = ' **** Matrix from chain '//ch1//' to '// $ ch2//' is not defined ****' call ccperr(2,errlin) istat = -5 return endif c ncsopr(ich3,ich4) = numncs endif endif return c Some error 999 istat = -1 return 998 call ccperr(2,' **** Matrix is not orthonormal ****') istat = -3 numncs = numncs-1 return end c c c subroutine chnnum(ch, ich) c ========================== c c Lookup chain number in list c c On entry: c ch chain identifier c c On exit: c ich chain number, = 0 if not found c character*1 ch integer ich c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh c integer i do 10, i = 1,nchain if (ch .eq. pchain(i)) then ich = i return endif 10 continue ich = 0 return end c c c subroutine cpymat(nfrom, nto, istat) c =================================== c c Copy NCS matrix number nfrom to nto c c On exit: c istat = 0 OK c = -1 index nfrom or nto out of range c c integer nfrom, nto, istat c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include ncscm.fh c# ncscm.fh c c NCS operators c c rtncs(4,4,k) k'th NCS operator: operator 1 is always identity c numncs number of NCS operators c ncsopr(i,j) pointer to operator to move chain i on chain j c = 0 if none c common /ncscm/ rtncs(4,4,maxncs), numncs, ncsopr(maxchn,maxchn) real rtncs integer numncs, ncsopr save /ncscm/ c## C&&*&& end_include ncscm.fh c integer i,j c istat = 0 if (nfrom .le. 0 .or. nto .le. 0) then istat = -1 return endif c do 10, j = 1,4 do 11, i = 1,4 rtncs(i,j,nto) = rtncs(i,j,nfrom) 11 continue 10 continue return end c c c subroutine redmol(IXYZIN) c ================= c c Read coordinate file, store each chain separately c c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) INTEGER IXYZIN c## C&&*&& end_include parameters.fh C&&*&& include waters.fh c# waters.fh c Character common c watres*4 water residue name c chnwat water chain names for each water atom c wchain possible water chain names for output c watnam atomname for water c uwchain unique water chain names (nqwchn of them) c ie water chain jch has naem uwchain(jch) c common /waterc/ watres, chnwat(maxwat), wchain(maxchn), $ uwchain(maxchn), watnam character watres*4, chnwat*1, wchain*1, uwchain*1, watnam*4 c Numbers etc c xwater orthogonal coordinates c fwater fractional coordinates c bwater Bfactor c qwater occupancy c distat(,ich) distance to nearest atom in chain ich c ctfat(3,,ich) cell translations to nearest atom in chain ich c xyzcrt(3,,ich) orthogonal coordinates of water moved by crystal c symmetry for protein chain ich c (ie lsymat(,ich), ctfat(3,,ich)) c nearat(,ich) nearest atom in chain ich c lsymat(,ich) symmetry number for nearest atom in chain ich c nearch nearest chain c nrswat residue number (input) c nwater number of waters c nwchain number of water chains c iwchain assigned chain number for water c mltplw multiplicity of this water c mltmax maximum multiplicity c kmswat "master" (ie 1st) water number related to this one c numwat input water number for this allocated water number c kwchain unique water chain number corresponding to protein chain c (allowing for same water chain to match multiple c protein chains) c nqwchn number of unique water chains c common /waters/ xwater(3,maxwat), fwater(3,maxwat), $ bwater(7,maxwat), qwater(maxwat),distat(maxwat,maxchn), $ ctfat(3,maxwat,maxchn), xyzcrt(3,maxwat,maxchn), $ nearat(maxwat,maxchn), lsymat(maxwat,maxchn), nearch(maxwat), $ nrswat(maxwat), nwater, nwchain, iwchain(maxwat), $ mltplw(maxwat), mltmax, kmswat(maxwat), $ numwat(maxwat,maxchn), kwchain(maxchn), nqwchn real xwater, fwater, bwater, qwater, distat, ctfat, xyzcrt integer nwater, nearch, nrswat, nwchain, nearat, lsymat, iwchain, $ mltplw, mltmax, kmswat, numwat, kwchain, nqwchn save /waterc/, /waters/ c## C&&*&& end_include waters.fh C&&*&& include flagcm.fh c# flagcm.fh c c Control flags etc c c dislim maximum distance check c dismol maximum allowed distance of water from "protein" molecule c dissim distance threshold for considering "equivalent" waters c disim2 dissim**2 c lonly .true. to check distances to 'O' & 'N' atoms only c verbose .true. for debug printing c common /flagcm/ dislim, dismol, dissim, disim2, lonly, verbose real dislim, dismol, dissim, disim2 logical lonly, verbose save /flagcm/ c## C&&*&& end_include flagcm.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh C&&*&& include cellcm.fh c# cellcm.fh c c Cell etc c c cell unit cell c rcell reciprocal cell (A**-1) C ro(4,4) fractional to orthogonal matrix C rf(4,4) orthogonal to fractional matrix C c common /cellcm/ cell(6), rcell(6), ro(4,4), rf(4,4) real cell, rcell, ro, rf save /cellcm/ C&&*&& end_include cellcm.fh c c Locals integer i, j, n, istat c Atom things integer iser, iresn, iz real biso, x(3), occ, b(6) character errlin*80, atnam*4, restyp*4, chnnam*1, resno*4, id*4, + inscod*1, altcod*1, segid*4 c nwater = 0 c Initialize coordinate input & output ixyzin = 0 c istat = 0 call XYZOPEN('XYZIN','INPUT',' ',IXYZIN,istat) c c c Get cell & orthogonalization matrices from XYZIN file call getorm(rf, ro, cell, istat) if (istat .ne. 3) then call ccperr(1,'*** XYZIN file must contain CRYST & SCALE ***') endif dislim = min(cell(1), cell(2), cell(3)) c c---- Get reciprocal cell c call reccel(cell,rcell) c c c 10 call XYZADVANCE(IXYZIN,0,0,*100,*100) call XYZATOM(IXYZIN,iser,atnam,restyp,chnnam,iresn,resno, $ inscod,altcod,segid,iz,id) call XYZCOORD(IXYZIN,'O','U',x(1),x(2),x(3),occ,biso,b) if (restyp .eq. watres) then c Water c If water atom name unset, copy from first water if (b(2).ne.0.0 .or. b(3).ne.0.0) then if (x(1).eq.0.0 .and. x(2).eq.0.0 .and. x(3).eq.0.0) then do 30 i=nwater,1,-1 if (chnnam .eq. chnwat(i)) then if (iresn .eq. nrswat(i)) then do 40 j=1,6 bwater(j,i) = b(j) 40 continue bwater(7,i) = biso goto 10 endif endif 30 continue goto 10 endif endif if (watnam .eq. ' ') then watnam = atnam endif nwater = nwater+1 if (nwater .gt. maxwat) then call ccperr(1,'*** Too many waters ***') endif do 12 i = 1,3 xwater(i,nwater) = x(i) 12 continue c c---- Store fractional coordinates c call matvc4(fwater(1,nwater),rf,x) c bwater(7,nwater) = biso do 20 i=1,6 bwater(i,nwater) = b(i) 20 continue qwater(nwater) = occ chnwat(nwater) = chnnam nrswat(nwater) = iresn else c Not water, don't store C atoms if lonly .true. if (lonly) then if (iz .eq. 6) then go to 10 endif endif c c---- Which chain? c do 13 n = 1,nchain if (chnnam .eq. pchain(n)) go to 14 13 continue c c errlin = ' **** Chain '//chnnam//' not recognized ****' call ccperr(2,errlin) call ccperr(1,'*** Unrecognized chain ***') c c 14 natmch(n) = natmch(n) + 1 c c---- Store orthogonal and fractional coordinates c do 15 i = 1,3 xprotn(i,natmch(n),n) = x(i) 15 continue c call matvc4(fprotn(1,natmch(n),n),rf,x) c atomnm(natmch(n),n) = atnam//resno//chnnam endif go to 10 100 write (6, '(//,1x,i8,a,/)') $ nwater, ' waters read' write (6, '(/,a,//,50(5x,a,a,a,i8,/))') $ ' Number of non-water atoms stored for each chain:', $ ('Chain ',pchain(i), ': ',natmch(i),i=1,nchain) return c end c c c subroutine closest c ================== c c For each water, find closest atom in each chain c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include waters.fh c# waters.fh c Character common c watres*4 water residue name c chnwat water chain names for each water atom c wchain possible water chain names for output c watnam atomname for water c uwchain unique water chain names (nqwchn of them) c ie water chain jch has naem uwchain(jch) c common /waterc/ watres, chnwat(maxwat), wchain(maxchn), $ uwchain(maxchn), watnam character watres*4, chnwat*1, wchain*1, uwchain*1, watnam*4 c Numbers etc c xwater orthogonal coordinates c fwater fractional coordinates c bwater Bfactor c qwater occupancy c distat(,ich) distance to nearest atom in chain ich c ctfat(3,,ich) cell translations to nearest atom in chain ich c xyzcrt(3,,ich) orthogonal coordinates of water moved by crystal c symmetry for protein chain ich c (ie lsymat(,ich), ctfat(3,,ich)) c nearat(,ich) nearest atom in chain ich c lsymat(,ich) symmetry number for nearest atom in chain ich c nearch nearest chain c nrswat residue number (input) c nwater number of waters c nwchain number of water chains c iwchain assigned chain number for water c mltplw multiplicity of this water c mltmax maximum multiplicity c kmswat "master" (ie 1st) water number related to this one c numwat input water number for this allocated water number c kwchain unique water chain number corresponding to protein chain c (allowing for same water chain to match multiple c protein chains) c nqwchn number of unique water chains c common /waters/ xwater(3,maxwat), fwater(3,maxwat), $ bwater(7,maxwat), qwater(maxwat),distat(maxwat,maxchn), $ ctfat(3,maxwat,maxchn), xyzcrt(3,maxwat,maxchn), $ nearat(maxwat,maxchn), lsymat(maxwat,maxchn), nearch(maxwat), $ nrswat(maxwat), nwater, nwchain, iwchain(maxwat), $ mltplw(maxwat), mltmax, kmswat(maxwat), $ numwat(maxwat,maxchn), kwchain(maxchn), nqwchn real xwater, fwater, bwater, qwater, distat, ctfat, xyzcrt integer nwater, nearch, nrswat, nwchain, nearat, lsymat, iwchain, $ mltplw, mltmax, kmswat, numwat, kwchain, nqwchn save /waterc/, /waters/ c## C&&*&& end_include waters.fh C&&*&& include flagcm.fh c# flagcm.fh c c Control flags etc c c dislim maximum distance check c dismol maximum allowed distance of water from "protein" molecule c dissim distance threshold for considering "equivalent" waters c disim2 dissim**2 c lonly .true. to check distances to 'O' & 'N' atoms only c verbose .true. for debug printing c common /flagcm/ dislim, dismol, dissim, disim2, lonly, verbose real dislim, dismol, dissim, disim2 logical lonly, verbose save /flagcm/ c## C&&*&& end_include flagcm.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh C&&*&& include symmcm.fh c# symmcm.fh c c Crystallographic symmetry c integer maxsym parameter (maxsym=96) c common /symmcm/ rsym(4,4,maxsym), rfsym(4,4,maxsym), $ nsym, nsymp, numsgp integer nsym, nsymp, numsgp real rsym, rfsym common /symmcc/ spgnam, pgname character spgnam*10, pgname*10 save /symmcm/, /symmcc/ c## C&&*&& end_include symmcm.fh C&&*&& include cellcm.fh c# cellcm.fh c c Cell etc c c cell unit cell c rcell reciprocal cell (A**-1) C ro(4,4) fractional to orthogonal matrix C rf(4,4) orthogonal to fractional matrix C c common /cellcm/ cell(6), rcell(6), ro(4,4), rf(4,4) real cell, rcell, ro, rf save /cellcm/ C&&*&& end_include cellcm.fh c c Locals integer i, jw, ich, mat, nch, lsym, nchcl real del, deld, ctf(3), xyz(3), delch c c c Loop waters do 100, jw = 1,nwater c c Loop chains delch = 10000000. nchcl = 0 do 110, ich = 1, nchain c Clear flags for this water nearat(jw, ich) = 0 distat(jw, ich) = dislim deld = dislim del = deld**2 mat = 0 nch = 0 call nearst(ich,fwater(1,jw),natmch(ich), + fprotn(1,1,ich),del,deld,mat,nch,lsym,ctf,xyz) if (lsym .gt. 0) then nearat(jw, ich) = mat distat(jw, ich) = deld lsymat(jw, ich) = lsym do 111, i = 1,3 ctfat(i, jw, ich) = ctf(i) xyzcrt(i, jw, ich) = xyz(i) 111 continue if (deld .lt. delch) then c Update closest chain to water jw delch = deld nchcl = ich endif endif c 110 continue c Store closest chain if (nchcl .eq. 0) then write (6, '(/,a,i8,/)') $ ' **** Help: no closest atom to water', jw call ccperr(1,'*** Closest error ***') endif nearch(jw) = nchcl c c 100 continue end c c c subroutine matchw c ================= c c Find matching waters c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) c## C&&*&& end_include parameters.fh C&&*&& include waters.fh c# waters.fh c Character common c watres*4 water residue name c chnwat water chain names for each water atom c wchain possible water chain names for output c watnam atomname for water c uwchain unique water chain names (nqwchn of them) c ie water chain jch has naem uwchain(jch) c common /waterc/ watres, chnwat(maxwat), wchain(maxchn), $ uwchain(maxchn), watnam character watres*4, chnwat*1, wchain*1, uwchain*1, watnam*4 c Numbers etc c xwater orthogonal coordinates c fwater fractional coordinates c bwater Bfactor c qwater occupancy c distat(,ich) distance to nearest atom in chain ich c ctfat(3,,ich) cell translations to nearest atom in chain ich c xyzcrt(3,,ich) orthogonal coordinates of water moved by crystal c symmetry for protein chain ich c (ie lsymat(,ich), ctfat(3,,ich)) c nearat(,ich) nearest atom in chain ich c lsymat(,ich) symmetry number for nearest atom in chain ich c nearch nearest chain c nrswat residue number (input) c nwater number of waters c nwchain number of water chains c iwchain assigned chain number for water c mltplw multiplicity of this water c mltmax maximum multiplicity c kmswat "master" (ie 1st) water number related to this one c numwat input water number for this allocated water number c kwchain unique water chain number corresponding to protein chain c (allowing for same water chain to match multiple c protein chains) c nqwchn number of unique water chains c common /waters/ xwater(3,maxwat), fwater(3,maxwat), $ bwater(7,maxwat), qwater(maxwat),distat(maxwat,maxchn), $ ctfat(3,maxwat,maxchn), xyzcrt(3,maxwat,maxchn), $ nearat(maxwat,maxchn), lsymat(maxwat,maxchn), nearch(maxwat), $ nrswat(maxwat), nwater, nwchain, iwchain(maxwat), $ mltplw(maxwat), mltmax, kmswat(maxwat), $ numwat(maxwat,maxchn), kwchain(maxchn), nqwchn real xwater, fwater, bwater, qwater, distat, ctfat, xyzcrt integer nwater, nearch, nrswat, nwchain, nearat, lsymat, iwchain, $ mltplw, mltmax, kmswat, numwat, kwchain, nqwchn save /waterc/, /waters/ c## C&&*&& end_include waters.fh C&&*&& include flagcm.fh c# flagcm.fh c c Control flags etc c c dislim maximum distance check c dismol maximum allowed distance of water from "protein" molecule c dissim distance threshold for considering "equivalent" waters c disim2 dissim**2 c lonly .true. to check distances to 'O' & 'N' atoms only c verbose .true. for debug printing c common /flagcm/ dislim, dismol, dissim, disim2, lonly, verbose real dislim, dismol, dissim, disim2 logical lonly, verbose save /flagcm/ c## C&&*&& end_include flagcm.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh C&&*&& include symmcm.fh c# symmcm.fh c c Crystallographic symmetry c integer maxsym parameter (maxsym=96) c common /symmcm/ rsym(4,4,maxsym), rfsym(4,4,maxsym), $ nsym, nsymp, numsgp integer nsym, nsymp, numsgp real rsym, rfsym common /symmcc/ spgnam, pgname character spgnam*10, pgname*10 save /symmcm/, /symmcc/ c## C&&*&& end_include symmcm.fh C&&*&& include ncscm.fh c# ncscm.fh c c NCS operators c c rtncs(4,4,k) k'th NCS operator: operator 1 is always identity c numncs number of NCS operators c ncsopr(i,j) pointer to operator to move chain i on chain j c = 0 if none c common /ncscm/ rtncs(4,4,maxncs), numncs, ncsopr(maxchn,maxchn) real rtncs integer numncs, ncsopr save /ncscm/ c## C&&*&& end_include ncscm.fh c c Locals integer jw, jch, kw, kch, nkch, nscore(maxchn), $ ncs, mkw, njwkw(maxchn,maxchn), n, jwch real x(3), dsim2, score(maxchn), dscore, ds2, dv(3) c real dot external dot c do 1, jw = 1,nwater iwchain(jw) = 0 kmswat(jw) = 0 do 2, jch = 1, nchain numwat(jw,jch) = 0 2 continue 1 continue c mltmax = 0 c Outer loop over waters do 10, jw = 1,nwater-1 if (iwchain(jw) .ne. 0) go to 10 c Loop possible chain allocations for water jw do 20, jch = 1, nchain c Only consider chains closer than dismol if (distat(jw, jch) .lt. dismol) then c+++ c For the hypothesis that water jw "belongs" to chain jch: c see if there are any other equivalent waters on other related chains c Loop chains, testing only those related by NCS nkch = 0 dscore = 0.0 do 30, kch = 1,nchain c mjwkw(jch,kch) water in kch which might match water in jch njwkw(jch,kch) = 0 if (kch .ne. jch .and. ncsopr(kch,jch) .gt. 0) then ncs = ncsopr(kch,jch) c--- Loop waters kw c Which is the "most similar" water attached to chain kch? dsim2 = disim2 mkw = 0 do 40, kw = jw+1, nwater if (iwchain(kw) .eq. 0) then c c Is water kw(chain kch) close enough to protein to be considered? if (distat(kw, kch) .lt. dismol) then c Is water NCS[kw(chain kch)] equivalent to water jw(chain jch)? call matvc4(x, rtncs(1,1,ncs), $ xyzcrt(1,kw,kch)) call vdif(dv, xyzcrt(1,jw,jch), x) ds2 = dot(dv,dv) c ds2 is distance**2 between NCS-transformed kw & jw if (ds2 .le. dsim2) then c Yes dsim2 = ds2 mkw = kw endif endif endif 40 continue c mkw is the water in chain kch most similar to jw/jch c--- End loop waters kw c add in score for chain kch if (mkw .gt. 0) then dscore = dscore + dsim2 nkch = nkch + 1 njwkw(jch,kch) = mkw endif endif 30 continue c+++ End loop chains kch c Score for chain jch for water jw nscore(jch) = nkch if (nkch .gt. 0) then score(jch) = sqrt(dscore/real(nkch)) else score(jch) = 0.0 endif else nscore(jch) = 0 score(jch) = 0.0 endif 20 continue c Which chain should jw belong to? c n = 0 jwch = 0 do 11, jch = 1,nchain if (nscore(jch) .gt. n) then dscore = 100000000. n = nscore(jch) endif if (nscore(jch) .eq. n) then if (score(jch) .lt. dscore) then dscore = score(jch) jwch = jch endif endif 11 continue c if (jwch .gt. 0) then c Assign water jw to "best" chain, & assign matching waters kw iwchain(jw) = jwch mltplw(jw) = n+1 kmswat(jw) = 0 do 12, kch = 1,nchain if (njwkw(jwch,kch) .gt. 0) then iwchain(njwkw(jwch,kch)) = kch mltplw(njwkw(jwch,kch)) = n+1 kmswat(njwkw(jwch,kch)) = jw endif 12 continue else c Waters with no matches are on their own, assign them to nearest c protein chain mltplw(jw) = 1 iwchain(jw) = nearch(jw) c endif mltmax = max(mltmax, mltplw(jw)) 10 continue c What about the last one! if (iwchain(nwater) .eq. 0) then mltplw(nwater) = 1 iwchain(nwater) = nearch(nwater) mltmax = max(mltmax, mltplw(nwater)) endif c c return end c c c subroutine wrimol(IXYZIN) c ========================= c c Copy all non-water atoms to output file, append reclassified waters c C&&*&& include parameters.fh c# parameters.fh c c maxwat maximum number of waters c maxchn maximum number of non-water chains c maxatm maximum number of atoms/chain c maxncs maximum number of NCS operators c integer maxwat, maxchn, maxatm, maxncs parameter (maxwat=2000, maxchn=12, maxatm=5000) parameter (maxncs=360) INTEGER IXYZIN c## C&&*&& end_include parameters.fh C&&*&& include waters.fh c# waters.fh c Character common c watres*4 water residue name c chnwat water chain names for each water atom c wchain possible water chain names for output c watnam atomname for water c uwchain unique water chain names (nqwchn of them) c ie water chain jch has naem uwchain(jch) c common /waterc/ watres, chnwat(maxwat), wchain(maxchn), $ uwchain(maxchn), watnam character watres*4, chnwat*1, wchain*1, uwchain*1, watnam*4 c Numbers etc c xwater orthogonal coordinates c fwater fractional coordinates c bwater Bfactor c qwater occupancy c distat(,ich) distance to nearest atom in chain ich c ctfat(3,,ich) cell translations to nearest atom in chain ich c xyzcrt(3,,ich) orthogonal coordinates of water moved by crystal c symmetry for protein chain ich c (ie lsymat(,ich), ctfat(3,,ich)) c nearat(,ich) nearest atom in chain ich c lsymat(,ich) symmetry number for nearest atom in chain ich c nearch nearest chain c nrswat residue number (input) c nwater number of waters c nwchain number of water chains c iwchain assigned chain number for water c mltplw multiplicity of this water c mltmax maximum multiplicity c kmswat "master" (ie 1st) water number related to this one c numwat input water number for this allocated water number c kwchain unique water chain number corresponding to protein chain c (allowing for same water chain to match multiple c protein chains) c nqwchn number of unique water chains c common /waters/ xwater(3,maxwat), fwater(3,maxwat), $ bwater(7,maxwat), qwater(maxwat),distat(maxwat,maxchn), $ ctfat(3,maxwat,maxchn), xyzcrt(3,maxwat,maxchn), $ nearat(maxwat,maxchn), lsymat(maxwat,maxchn), nearch(maxwat), $ nrswat(maxwat), nwater, nwchain, iwchain(maxwat), $ mltplw(maxwat), mltmax, kmswat(maxwat), $ numwat(maxwat,maxchn), kwchain(maxchn), nqwchn real xwater, fwater, bwater, qwater, distat, ctfat, xyzcrt integer nwater, nearch, nrswat, nwchain, nearat, lsymat, iwchain, $ mltplw, mltmax, kmswat, numwat, kwchain, nqwchn save /waterc/, /waters/ c## C&&*&& end_include waters.fh C&&*&& include flagcm.fh c# flagcm.fh c c Control flags etc c c dislim maximum distance check c dismol maximum allowed distance of water from "protein" molecule c dissim distance threshold for considering "equivalent" waters c disim2 dissim**2 c lonly .true. to check distances to 'O' & 'N' atoms only c verbose .true. for debug printing c common /flagcm/ dislim, dismol, dissim, disim2, lonly, verbose real dislim, dismol, dissim, disim2 logical lonly, verbose save /flagcm/ c## C&&*&& end_include flagcm.fh C&&*&& include coorcm.fh c# coorcm.fh c c Coordinate store for non-waters c c xprotn orthogonal coordinates c fprotn fractional coordinates c natmch number of atoms in each chain c atomnm atom name//residue number//chain ID c pchain chain ID for each chain c nchain number of chains c common /coorcm/ xprotn(3,maxatm,maxchn),fprotn(3,maxatm,maxchn), $ natmch(maxchn), nchain real xprotn, fprotn integer natmch, nchain c common /coorcc/ atomnm(maxatm,maxchn), pchain(maxchn) character atomnm*12, pchain*1 save /coorcm/, /coorcc/ c## C&&*&& end_include coorcm.fh c c Atom things integer iser, iresn, iz real biso, btmp(6) character atnam*4, restyp*4, chnnam*1, id*4, resno*4, segid*4, $ inscod*1, altcod*1 c integer ifail, jwater, m, iw, jw, kw, jch, nw, j, k, kch, $ jwch, kwch, IXYZOUT, ndupl real dist, dismin c c call XYZREWD(IXYZIN) c ifail = 0 IXYZOUT = 0 call XYZOPEN('XYZOUT','OUTPUT',' ',IXYZOUT,ifail) c iser = 0 c 10 call XYZADVANCE(IXYZIN,IXYZOUT,0,*100,*100) call XYZATOM(IXYZIN,j,atnam,restyp,chnnam,iresn,resno, + inscod,altcod,segid,iz,id) if (restyp(1:3) .ne. watres) then iser = iser + 1 c call XYZATOM(IXYZOUT,iser,atnam,restyp,chnnam,iresn, + resno,inscod,altcod,segid,iz,id) call XYZADVANCE(IXYZOUT,0,0,*100,*100) endif go to 10 c c Assign water numbers c Loop in order of decreasing multiplicity c c Water serial number 100 continue c Check through waters to remove duplicates c Minimum allowed distance between waters to eliminate as duplicates dismin = 0.1 ndupl = 0 do jw = 2, nwater do iw = 1, jw-1 c Calculate distance between water call disclc(xyzcrt, maxwat, maxchn, iwchain, iw, jw, dist) if (dist .lt. dismin) then c These two waters are coincident, remove the second one mltplw(iw) = 0 write (6, 6005) $ nrswat(iw), chnwat(iw), nrswat(jw), chnwat(jw) 6005 format (' Water ',i4,1x,a, $ ' removed as it duplicates water ',i4,1x,a) ndupl = ndupl+1 endif enddo enddo write (6, '(/1x,i6,a/)') ndupl, ' duplicate waters removed' c jwater = 0 do 20, m = mltmax,1,-1 do 30, jw = 1, nwater if (kmswat(jw) .eq. 0 .and. mltplw(jw).eq. m) then jwater = jwater+1 jch = iwchain(jw) jwch = kwchain(jch) numwat(jwater,jwch) = jw j = nearat(jw,jch) biso = bwater(7,jw) c List principal contacts write (6, 6020) jwater, m 6020 format (/,' Water residue ',i8,' multiplicity ',i4) write (6, 6030) jch,wchain(jch), $ nrswat(jw),chnwat(jw),biso, $ atomnm(j,jch),distat(jw,jch) 6030 format (7x,'Chain ',i4,'(',a,') original',i7,a, $ ' B-factor',f8.2,' nearest atom ',a, $ 'distance',f8.2) if (jw .lt. nwater) then do 40, kw = jw+1, nwater if (kmswat(kw) .eq. jw) then kch = iwchain(kw) kwch = kwchain(kch) numwat(jwater,kwch) = kw k = nearat(kw,kch) biso = bwater(7,kw) write (6, 6030) kch,wchain(kch), $ nrswat(kw),chnwat(kw),biso, $ atomnm(k,kch),distat(kw,kch) endif 40 continue endif endif 30 continue 20 continue c write (6,'(//,1x,i8,a,/)') jwater,' residues assigned' c c c Append waters, output by chain c nw = 0 id = ' O ' segid = ' ' inscod = ' ' altcod = ' ' resno = ' ' do 50, jch = 1, nqwchn do 60, kw = 1,jwater jw = numwat(kw,jch) if (jw .gt. 0) then iser = iser + 1 do 70 j = 1,6 btmp(j) = bwater(j,jw) 70 continue biso = bwater(7,jw) c call XYZATOM(IXYZOUT,iser,watnam,watres, $ uwchain(jch),kw,resno,inscod,altcod,segid,iz,id) call XYZCOORD(IXYZOUT,'O','U', $ xyzcrt(1,jw,iwchain(jw)),xyzcrt(2,jw,iwchain(jw)), $ xyzcrt(3,jw,iwchain(jw)),qwater(jw),biso,btmp) call XYZADVANCE(IXYZOUT,0,0,*100,*100) nw = nw+1 endif 60 continue 50 continue write (6,'(/,1x,i8,a,/)') nw, ' waters output' call XYZCLOSE(IXYZIN) call XYZCLOSE(IXYZOUT) return end c c c c subroutine getorm(rf, ro, cell, lflag) c =========================================== c c 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 c c c integer lflag real rf(4,4), ro(4,4), cell(6) c real vol c c lflag = 0 c ro(1,1) = 0.0 rf(1,1) = 0.0 cell(1) = 0.0 call rbrorf(ro,rf) call rbcell(cell,vol) if (cell(1) .gt. 0.0) lflag = 1 if (ro(1,1) .gt. 0.0) then lflag= lflag + 2 endif c return end c c subroutine nearst(n,fwater,jat,fprotn,del,deld,mat,nch, + lsym,ctf,xyz) c ============================================================== c c c c Search all atoms of chain N to find nearest to water XWATER (FWATER) c c Returns MAT = atom number, NCH = chain number of nearest so far, c LSYM crystallographic symmetry number c c On entry: c n chain number c fwater fractional coordinates of water c jat number of atoms in chain c fprotn(3,jat) fractional coordinates of protein c del (maximum) distance**2 (= deld**2) c deld (maximum) distance c mat atom number, = 0 c nch (closest chain so far) c c On exit: c del distance**2 (= deld**2) c deld distance to nearest atom c mat closest atom number in this chain c nch closest chain so far c lsym crystallographic symmetry number to get close c ctf(3) cell translations to get close c xyz(3) orthogonal coordinates of water after applying c crystallographic symmetry lflag/ctf c c c .. Scalar Arguments .. real del,deld integer jat,lsym,mat,n,nch c .. c .. Array Arguments .. real ctf(3),fprotn(3,*),fwater(3),xyz(3) c .. c .. Local Scalars .. integer lflag,m c .. c .. Local Arrays .. real fdel(3) c .. c .. External Subroutines .. external frclim,tstdis c .. c c---- Set fractional limits equivalent to DELD c c ***************** call frclim(deld,fdel) c ***************** c c---- Loop all protein atoms in this chain c do 10 m = 1,jat c c ***************************************************** call tstdis(fprotn(1,m),del,deld,fdel,fwater,lflag,ctf,xyz) c ***************************************************** c if (lflag.gt.0) then mat = m nch = n lsym = lflag end if 10 continue c c---- End atom loop c end c c c subroutine tstdis(fprotn,del,deld,fdel,fwater,lflag,ctf,xyz) c ============================================================ c c Test if atom at fractional coordinate FPROTN is nearer to atom c FWATER than current closest defined by DEL, DELD, FDEL c c Returns LFLAG = LSYM and CTF = cell translation from water to protein c if new atom found, else = 0 c c On entry: c fprotn(3) fractional coordinates of protein atom c del (maximum) distance**2 (= deld**2) c deld (maximum) distance c fdel(3) fractional limits equivalent to deld c fwater fractional coordinates of water c c On exit: c del distance**2 (= deld**2) c deld distance to nearest atom c fdel(3) fractional limits equivalent to deld c lflag crystallographic symmetry number applied if closest, c else = 0 c ctf(3) cell translations to get close c xyz(3) orthogonal coordinates of water after applying c crystallographic symmetry lflag/ctf c c c .. Scalar Arguments .. real del,deld integer lflag c .. c .. Array Arguments .. real ctf(3),fdel(3),fprotn(3),fwater(3),xyz(3) c .. c .. Local Scalars .. real a,d,fd,fr integer i,js c .. c .. Local Arrays .. real ct(3),x(3),y(3),xprotn(3) c .. c .. External Subroutines .. external frclim,vsum c .. c .. Intrinsic Functions .. intrinsic abs,mod,sqrt c .. c .. Common blocks .. C&&*&& include symmcm.fh c# symmcm.fh c c Crystallographic symmetry c integer maxsym parameter (maxsym=96) c common /symmcm/ rsym(4,4,maxsym), rfsym(4,4,maxsym), $ nsym, nsymp, numsgp integer nsym, nsymp, numsgp real rsym, rfsym common /symmcc/ spgnam, pgname character spgnam*10, pgname*10 save /symmcm/, /symmcc/ c## C&&*&& end_include symmcm.fh C&&*&& include cellcm.fh c# cellcm.fh c c Cell etc c c cell unit cell c rcell reciprocal cell (A**-1) C ro(4,4) fractional to orthogonal matrix C rf(4,4) orthogonal to fractional matrix C c common /cellcm/ cell(6), rcell(6), ro(4,4), rf(4,4) real cell, rcell, ro, rf save /cellcm/ C&&*&& end_include cellcm.fh c .. c c Orthogonalize protein coordinates call matvc4(xprotn, ro, fprotn) c lflag = 0 c c---- Loop crystallographic symmetry c do 40 js = 1,nsym c c---- Apply symmetry to fractional coordinates of water c c **************************** call matvc4(x,rsym(1,1,js),fwater) c **************************** c c---- Check for atom in fractional box of size +-DELD (FDEL) c do 10 i = 1,3 fr = x(i) - fprotn(i) + fdel(i) fd = mod(fr,1.0) if (fd.lt.0.0) fd = fd + 1.0 c c---- Jump out immediately if too far away c if (fd.gt.2.0*fdel(i)) then go to 40 else c c---- Store cell translation used c ct(i) = fd - fr end if 10 continue c c---- Atom is within current fractional box (DELD, FDEL) c Check distance in orthogonal space c Apply cell translation and convert to orthogonal coordinates c call vsum(x,x,ct) call matvc4(y, ro, x) c a = 0.0 c c do 20 i = 1,3 d = y(i) - xprotn(i) c c if (abs(d).gt.deld) then go to 40 else a = d*d + a end if 20 continue c c if (a.lt.del) then c c---- New closest atom found c del = a deld = sqrt(a) lflag = js c c---- Reset fractional limits c c ***************** call frclim(deld,fdel) c ***************** c c---- store cell translation from water to protein c do 30 i = 1,3 ctf(i) = +ct(i) xyz(i) = y(i) 30 continue c c end if 40 continue c c---- End symmetry loop c end c c c subroutine frclim(deld,fdel) c ============================ c c Set fractional limits FDEL equivalent to distance DELD c = DELD . a*, .b*, .c* c c c .. Scalar Arguments .. real deld c .. c .. Array Arguments .. real fdel(3) c .. c .. Local Scalars .. integer i c .. c .. Common blocks .. C&&*&& include cellcm.fh c# cellcm.fh c c Cell etc c c cell unit cell c rcell reciprocal cell (A**-1) C ro(4,4) fractional to orthogonal matrix C rf(4,4) orthogonal to fractional matrix C c common /cellcm/ cell(6), rcell(6), ro(4,4), rf(4,4) real cell, rcell, ro, rf save /cellcm/ C&&*&& end_include cellcm.fh c .. c c do 10 i = 1,3 fdel(i) = rcell(i)*deld 10 continue c c end c c c subroutine reccel(cell,rcell) c ============================= c c Calculation of reciprocal cell parameters c c c .. Array Arguments .. real cell(6),rcell(6) c .. c .. Local Scalars .. real alpha,ar,beta,br,const,cosa,cosas,cosb,cosbs,cosg,cosgs,cr, + gamma,s,sina,sinb,sinbs,sing,vol c .. c .. Intrinsic Functions .. intrinsic cos,sin,sqrt c .. c c const = 3.1415926/180.0 alpha = cell(4)*const beta = cell(5)*const gamma = cell(6)*const cosa = cos(alpha) sina = sin(alpha) cosb = cos(beta) sinb = sin(beta) cosg = cos(gamma) sing = sin(gamma) cosas = (cosb*cosg-cosa)/ (sinb*sing) cosbs = (cosa*cosg-cosb)/ (sina*sing) cosgs = (cosa*cosb-cosg)/ (sina*sinb) sinbs = sqrt(1.0-cosbs*cosbs) ar = 1.0/ (cell(1)*sinb*sqrt(1.0-cosgs*cosgs)) br = 1.0/ (cell(2)*sing*sqrt(1.0-cosas*cosas)) cr = 1.0/ (cell(3)*sina*sinbs) rcell(1) = ar rcell(2) = br rcell(3) = cr rcell(4) = cosas rcell(5) = cosbs rcell(6) = cosgs c c---- Cell volume c s = (alpha+beta+gamma)*0.5 vol = cell(1)*2.0*cell(2)*cell(3)* + sqrt(sin(s-alpha)*sin(s)*sin(s-beta)*sin(s-gamma)) c c end c c c 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 c 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 = +1 blank line c c integer lunit, nitems, istat character*(*) file, blknam, type, format 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 c c istat = 0 call ccpdpn(lunit, file, 'READONLY', 'F', 0,0) c 1 read (lunit, '(a)') line if (line(1:1) .eq. '!' .or. line(1:1) .eq. '#') go to 1 c ntok=maxtok call parser( . key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok, $ lend,.false.) if(lend) then istat = -1 go to 900 endif if(ntok.eq.0) then istat = +1 go to 900 endif c blknam = line(ibeg(1):iend(1)) type = line(ibeg(2):iend(2)) nitems = nint(fvalue(3)) format = line(ibeg(4):) c 900 return end c c subroutine chkmt4(a, istat) c =========================== c c Check 4x4 matrix for orthogonality, c ie that 3x3 rotation part is orthonormal c c Returns istat = 0 OK, -1 not c real a(4,4) integer istat c real r(3,3), s(3,3), p(3,3), eps integer i,j external matmul data eps/0.0001/ c c r is 3x3 bit, s is r(transpose) do 1, j = 1,3 do 2, i = 1,3 r(i,j) = a(i,j) s(i,j) = a(j,i) 2 continue 1 continue c call matmul(p,r,s) c istat = 0 do 10, j = 1,3 do 11, i = 1,3 if(i .eq. j) then if (abs(p(i,j)-1.0) .gt. eps) istat = -1 else if (abs(p(i,j)) .gt. eps) istat = -1 endif 11 continue 10 continue return end c c c subroutine setim4(a) c ==================== c c Set 4x4 matrix a to identity c real a(4,4) c integer i,j c do 1, j = 1,4 do 2, i = 1,4 a(i,j) = 0.0 2 continue a(j,j) = 1.0 1 continue return end c c c subroutine disclc(xyzcrt, maxwat, maxchn, iwchain, iw, jw, dist) c =============================================================== c c c Calculate distance between 2 waters c c On entry: c c xyzcrt(3,maxwat,maxchn) array of coordinates c maxwat, maxchn dimensions of array c iwchain(maxwat) array of chains for each water c iw number of first water c jw number of second water c c On exit: c dist distance between waters iw & jw c integer maxwat, maxchn, iw, jw real xyzcrt(3, maxwat, maxchn), dist integer iwchain(maxwat) c integer i real d c d = 0.0 do 10, i = 1,3 d = d + $ (xyzcrt(i,iw,iwchain(iw)) - xyzcrt(i,jw,iwchain(jw)))**2 10 continue c if (d .gt. 0.0) then dist = sqrt(d) else dist = 0.0 endif return end