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  MAKEDICT
c
c  Make dictionary entries, for TNT & PROTIN
c
c  Reads a set of coordinates for the group, in PDB or PROTIN format,
c and writes out dictionary entry in format for either TNT or PROTIN.
c The program can also convert a PROTIN coordinate entry to PDB format.
c Multiple copies of the group may be present (with different residue
c numbers), in which case bond lengths & angles will be averaged (only
c useful for TNT output).
c
c Commands:-
c
c-> INPUT PROTIN | PDB
c     Input format, Protin or PDB
c-> OUTPUT TNT | PROTIN <numres> <rscode> | PDB
c     TNT     Write output dictionary for TNT
c     PROTIN  Write output dictionary for Protin
c             numres    residue number for Protin
c             rscode    1-character residue name for Protin
c     PDB     Output file in PDB format
c                (eg for converting from Protin dictionary form)
c-> RADII [<element_name> <radius>]
c     Set bonding radius for elements (C,N,O,S & P are present by default)
c-> RESIDUE <residue_name>
c     New residue type name, to override input
c-> PLANE <atom list>
c     Define list of atoms in plane
c-> CHIRAL <central atom> <atom2> <atom3> <atom4>
c     List of atoms in chiral group
c     Note the TNT convention for chiral groups:
c       1st atom is central atom
c       Look from central atom towards the lowest-priority atom
c        (nearly always H), then list the other atoms (which are towards you)
c        in a clockwise order
c-> CENTRE <atom name>
c     Atom to consider as centre [default 1st atom]
c---
c--- For TNT output
c-> SD <sd_bond> <sd_angle> <sd_plane> <sd_bfactor>
c     Read sd parameters for TNT dictionary
c---
c--- For Protin output
c-> NONBOND [<type>] <atom pair 1> <atom pair 2> ...
c    List of atomname pairs to be put in the non-bonded contact list
c    Type is either 1 for single-torsion contacts or 2 [default] for
c     multiple-torsion contacts
c
c     Phil Evans
c     MRC Laboratory of Molecular Biology
c     Hills Road
c     Cambridge
c
c     March 1995
c
      integer istat
      character*2 iatm(100)
      external angles,bonds,ccperr,ccpfyp,ccprcs,chirals,contac,
     .         planes,rcards,redmol,wrtout
c
      data iatm/' H','HE','LI','BE',' B',' C',' N',' O',' F','NE',
     *   'NA','MG','AL','SI',' P',' S','CL','AR',' K','CA',
     *   'SC','TI',' V','CR','MN','FE','CO','NI','CU','ZN',
     *   'GA','GE','AS','SE','BR','KR','RB','SR',' Y','ZR',
     *   'NB','MO','TC','RU','RH','PD','AG','CD','IN','SN',
     *   'SB','TE',' I','XE','CS','BA','LA','CE','PR','ND',
     *   'PM','SM','EU','GD','TB','DY','HO','ER','TM','YB',
     *   'LU','HF','TA',' W','RE','OS','IR','PT','AU','HG',
     *   'TL','PB','BI','PO','AT','RN','FR','RA','AC','TH',
     *   'PA',' U','NP','PU','AM','CM','BK','CF','ES','FM'/
c
      call ccpfyp
      call ccprcs(6, 'MAKEDICT', '$Date: 2002/11/04 16:41:23 $')
      call XYZINIT
c
c  Read control commands and set radii
      call rcards(iatm)
c
c  Read in molecule
      call redmol
c
c  Calculate bonds & set up connectivity list
      call bonds
c
c  Calculate bond angles
      call angles
c
      call planes(istat)
      if (istat .ne. 0) then
         call ccperr(1,'**** Failure in planes routine ****')
      endif
c
      call chirals(istat)
      if (istat .ne. 0) then
         call ccperr(1,'**** Failure in chirals routine ****')
      endif
c
      call contac(istat)
      if (istat .ne. 0) then
         call ccperr(1,'**** Failure in contacts routine ****')
      endif
c
c  Write out information in a suitable form
      call wrtout(iatm)
c
      call ccperr(0,'*** Normal termination ***')
      end
c
c
c
      subroutine rcards(iatm)
c     ======================
c
c  Read control input
c
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c Arguments
      character*2 iatm(100)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
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
      character elmtyp*2
      integer lunout, ikey, idum, ierr, i, j, m, n, lflag, jat, jtype
c
      logical ccponl
      external ccponl
c
c Commands
      integer nkeys
      parameter (nkeys=10)
      character*16 keys(nkeys)
c
c
c ---- Data ----
      data keys/'INPUT','OUTPUT','RADII','RESIDUE',
     $   'PLANE','CHIRAL','SD','CENTRE','NONBOND','END'/
c
c
c  Set default atom radii
c cf DISTANG       Default values are: CA 1.5 C 1.5 N 1.5 O 1.5 SG 1.5 OW 1.5
      atmtyp(1) = 'C'
      atmrad(1) = 1.0
      iatmz(1)  =  6
      atmtyp(2) = 'O'
      atmrad(2) = 1.0
      iatmz(2)  =  8
      atmtyp(3) = 'N'   
      atmrad(3) = 1.0
      iatmz(3)  =  7
      atmtyp(4) = 'S'
      atmrad(4) = 1.5
      iatmz(4)  = 16
      atmtyp(5) = 'P'
      atmrad(5) = 1.5
      iatmz(5)  = 15
      natmtyp = 5
c Equivalence between Protin codes & atomnumbers
c  1=C(6)
c  2=N(7)
c  3=O(8)
c  4=S(16)
c  7=P(15)
      iprtyp(1) = 1
      iprtyp(2) = 3
      iprtyp(3) = 2
      iprtyp(4) = 4
      iprtyp(5) = 7
c
c Other defaults
      loutyp = 1
      lintyp = 1
      sdbond = 0.02
      sdbang = 3.0
      sdplane = 0.02
      sdbfac = 5.0
      ierr = 0
      lunout = 6
      residu = ' '
      numres = 1
      rscode = ' '
c
      nplanes  = 0
      nchirals = 0
      cenatm = ' '
      nnonbond = 0
c
c--------- Read input commands
 1    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 1
c
      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 1
       else
        ierr=ierr+1
       endif
      else if(ikey.lt.0) then
       write(lunout,fmt='(a)') ' ** Ambiguous keyword **'
       if(ccponl(idum)) then
        goto 1
       else
        ierr=ierr+1
       endif
      endif
c
      do 11, i=1,ntok
       call ccpupc(cvalue(i))
 11   continue
c
      go to (100, 200, 300, 400, 500, 600, 700, 800, 900, 1000), ikey
c
c============================================================
c
c     INPUT PROTIN | PDB
c  Input type
c
 100  if (cvalue(2) .eq. 'PDB') then
       lintyp = 1
      elseif (cvalue(2) .eq. 'PROT') then
       lintyp = 2
      else
       write (6,*) ' **** Unrecognized input type ****'
       ierr = ierr+1
      endif
      go to 1
c
c============================================================
c
c     OUTPUT TNT | PROTIN <numres> <rscode> | PDB
c
 200  n = 2
      if (cvalue(n) .eq. 'TNT') then
c
c     TNT
c     write output dictionary in TNT form
         loutyp = 1
c
      elseif (cvalue(n) .eq. 'PROT') then
c     PROTIN <numres> <rscode>
c  Write output dictionary for Protin
c     numres    residue number for Protin
c     rscode    1-character residue name for Protin

         loutyp = 2
         if (ntok .gt. n) then
            n = n+1
            call gttint(n, numres, lflag, ntok,ityp,fvalue)
            if (lflag .ne. 0) then
               ierr = ierr+1
            endif
         endif
         if (ntok .gt. n) then
            n = n+1
            rscode = line(ibeg(n):ibeg(n))
         endif
      elseif (cvalue(n) .eq. 'PDB') then
c     PDB
c        Write PDB output file to XYZOUT
         loutyp = 3
      endif
      go to 1
c
c============================================================
c
c     RADII [<element_name> <radius>]
c
 300  do 301, n = 2,ntok,2
       if (iend(n)-ibeg(n) .eq. 0) then
c  Element name is uppercase
c  Single character name, right justify
        elmtyp = ' '//cvalue(n)(1:1)
       else
c    else take 2 characters
        elmtyp = cvalue(n)(1:2)
       endif
c Find atomic number from table
       do 302, i = 1, 100
        if (elmtyp .eq. iatm(i)) go to 303
 302   continue
       write (6, '(/,a,a,a,/)') ' Element name ', elmtyp,
     $    ' not recognized'
       ierr = ierr+1
       go to 1
c
 303   do 304, j = 1, natmtyp
        if (iatmz(j) .eq. i) then
         jat = j
         go to 305
        endif
 304   continue 
       natmtyp = natmtyp+1
       if (natmtyp .gt. maxtyp) then
        write (6, '(/,a,i6,a,/)')
     $     ' **** Too many atom types, maximum ',maxtyp,
     $     ' ****'
        ierr = ierr+1
        go to 1
       endif
       jat = natmtyp
c
 305   atmtyp(jat) = elmtyp
       iatmz(jat)  = i
       call gttrea(n+1, atmrad(jat), lflag, ntok,ityp,fvalue)
 301  continue 
      go to 1
c============================================================
c
c     RESIDUE <residue_name>
c  Override residue type name
c
 400  residu = line(ibeg(2):iend(2))
      go to 1
c============================================================
c
c     PLANE <atom list>
 500  nplanes = nplanes+1
      if (nplanes .gt. maxpln) then
         call ccperr(1,'**** Too many planes ****')
      endif
      do 501, i = 2, ntok
         atmpln(i-1, nplanes) = line(ibeg(i):iend(i))
 501  continue
      natpln(nplanes) = ntok-1
      go to 1
c============================================================
c
c     CHIRAL <central atom> <atom2> <atom3> <atom4>
c     List of atoms in chiral group
 600  nchirals = nchirals+1
      if (nchirals .gt. maxcrl) then
         call ccperr(1,'**** Too many chirals ****')
      endif
      if (ntok .ne. 5) then
         call ccperr(2,' *** CHIRAL needs 4 atoms defined ***')
         ierr = ierr+1
         go to 1
      endif
      do 601, i = 2, ntok
         atmcrl(i-1, nchirals) = line(ibeg(i):iend(i))
 601  continue
      go to 1
c============================================================
c
c     SD <sd_bond> <sd_angle> <sd_plane> <sd_bfactor>
c     Read sd parameters for TNT dictionary
c
 700  call gttrea(2, sdbond, lflag, ntok,ityp,fvalue)
      call gttrea(3, sdbang, lflag, ntok,ityp,fvalue)
      call gttrea(4, sdplane, lflag, ntok,ityp,fvalue)
      call gttrea(5, sdbfac, lflag, ntok,ityp,fvalue)
      go to 1
c============================================================
c
c     CENTRE <atom name>
c  Name of "central" atom
c
 800  if (ntok .ne. 2) then
         call ccperr(2,' *** No CENTRAL atom defined ***')
         ierr = ierr+1
         go to 1
      endif
      cenatm = line(ibeg(2):iend(2))
      go to 1
c
c============================================================
c
c     NONBOND [<type>] <atomname pairs>
c
 900  m = 2
      if (ityp(2) .eq. 2) then
c  number read
         call gttint(2, jtype, lflag, ntok,ityp,fvalue)
         if (jtype .ne. 1 .and. jtype .ne. 2) then
            call ccperr(2,' **** Type must be 1 or 2 ****')
            ierr = ierr+1
            go to 1
         endif
         m = 3
      endif
      if (mod(ntok,2) .eq. mod(m,2)) then
         call ccperr(2,
     +            ' *** Must give pairs of atoms in nonbond list ***')
         ierr = ierr+1
         go to 1
      endif
      do 901, n = m, ntok, 2
         nnonbond = nnonbond+1
         if (nnonbond .gt. maxnnb) then
            call ccperr(1,'*** Too many nonbonds ***')
         endif
         atnonb(1,nnonbond) = line(ibeg(n):iend(n))
         atnonb(2,nnonbond) = line(ibeg(n+1):iend(n+1))
         ktnonb(nnonbond) = jtype
 901  continue 
      go to 1
c
c============================================================
c
c     END
c
 1000 continue

      if (ierr .gt. 0) then
       call ccperr(1,'*** Input error ***')
      endif
c
c
      if (lintyp .eq. 1) then
         write (6,'(/,a)')
     $      ' Input file is in PDB format'
      elseif (lintyp .eq. 2) then
         write (6,'(/,a)')
     $      ' Input file is in Protin format'
      endif
c
      if (loutyp .eq. 1) then
         write (6, '(/,a)') ' Output dictionary is in TNT format'
c
         write (6,'(/,a,2(/,a,f10.3),/)')
     $      ' Overall standard deviations for TNT dictionary:',
     $      '    Bond lengths:   ',sdbond,
     $      '    Bond angles :   ',sdbang,
     $      '    B-factors   :   ',sdbfac
c
      elseif (loutyp .eq. 2) then
         write (6, '(/,a)') ' Output dictionary is in PROTIN format'
      elseif (loutyp .eq. 3) then
         write (6, '(/,a)') ' Output file in PDB format'
      endif
      return
      end
c
c
c
      subroutine redmol
c     =================
c
c  Read molecule
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
c
c  Atom things
      integer iser, iresn, iz
      real xx, yy, zz, occ, biso, b(6)
      character atnam*4, resnam*4, chnnam*1, resno*4, id*4, inscod*1,
     $          altcod*1, segid*4
c
      character*80 errlin
c
c  Locals
      integer i, ifail, j , lres, ires, iun, ixyzin, natres,
     $   jatom
      logical newatm
c===================================
c
      kcenat = 1
c
      natom = 0
      lres  = -1000
      rsname = ' '
      ifail = 0
      ixyzin = 0
      if (lintyp .eq. 1) then
       call XYZOPEN('XYZIN','INPUT',' ',IXYZIN,IFAIL)
      elseif (lintyp .eq. 2) then
       iun = 0
       call ccpdpn(iun, 'XYZIN','READONLY','F',0,ifail)
       read (iun, '(50x,i5,10x,A3)') natres, resnam
      endif
c
      do 8, j = 1,maxres
       do 9, i = 1,maxatm
        latom(i,j) = .false.
 9     continue 
 8    continue 
c
c Read coordinate records
C
 10   if (lintyp .eq. 1) then
       call XYZADVANCE(IXYZIN,0,0,*100,*100)
       call XYZATOM(IXYZIN,iser,atnam,resnam,chnnam,iresn,resno,
     $                inscod,altcod,segid,iz,id)
       call XYZCOORD(IXYZIN,'O','U',xx,yy,zz,occ,biso,b)
       if (b(2).eq.0.0 .or. b(3).eq.0.0) then
         if (occ.eq.0.0) go to 10
C        if (xx.eq.0.0 .and. yy.eq.0.0 .and. zz.eq.0.0) goto 10
       endif
      elseif (lintyp .eq. 2) then
       read (iun, '(f8.3,2x,f8.3,2x,f8.3,12x,i5,30x,a3)',err=100)
     $    xx,yy,zz,iz,atnam
       iresn = 1
       do 11, i = 1, natmtyp
          if (iz .eq. iprtyp(i)) then
             iz = iatmz(i)
             go to 12
          endif
 11    continue
       write (6, '(/,a,i6,a,/)') ' **** Atom type ',iz,' not known ***'
       iz = 1
 12    continue
      endif
c
      if (lres .gt. 0) then
       if (lres .ne. iresn) then
c  Increment residue number
        ires = ires+1
        if (ires .gt. maxres) then
         write(errlin,fmt='(a,i10,a)')
     $      ' **** Maximum number of residues allowed is: ',maxres,
     $      ' Too many residues ***'
         call ccperr(1,errlin)
        endif
        lres = iresn
       endif
      else
       lres = iresn
       ires = 1
      endif
c      
c
      if (rsname .ne. ' ') then
       if (resnam .ne. rsname) then
        call ccperr(3,'*** WARNING: residue name has changed ***')
       endif
       rsname = resnam
      else
       rsname = resnam
      endif
c     
      if (ires .eq. 1) then
c  First residue
       newatm = .true.
      else
c  Not first residue, check for atom already seen
       do 5, i = 1, natom
        if (atnam .eq. atname(i)) then
         jatom = i
         newatm = .false.
         go to 6
        endif
 5     continue
       newatm = .false.
      endif
c
 6    if (newatm) then
c  This is a new atom, store it
       natom = natom+1
       if (natom .gt. maxatm) then
        write(errlin,fmt='(a,i10,a)')
     +       ' **** Too many atoms, maximum ',maxatm,' ****'
        call ccperr(1,errlin)
       endif
       jatom = natom
       atname(jatom) = atnam
       attype(jatom) = atnam(1:2)
       do 1, i = 1, natmtyp
        if (iz .eq. iatmz(i)) then
         lattyp(jatom) = i
         go to 3
        endif
 1     continue
       call ccperr(1,' Atom type not recognized: '//atnam)
      endif
c
c  Store coordinates
 3    continue
      xyz(1,jatom,ires) = xx
      xyz(2,jatom,ires) = yy
      xyz(3,jatom,ires) = zz

      latom(jatom,ires) = .true.
c Check for central atom
      if (cenatm .ne. ' ') then
         if (cenatm .eq. atnam) then
            kcenat = jatom
         endif
      endif
c
      if (lintyp .eq. 2) then
       if (natom .lt. natres) go to 10
      endif
      go to 10
c
c================= end of coordinate read
 100  nresid = ires
      write (6, '(/,1x,i8,a,/)') nresid, ' residues read', 
     $   natom, ' atoms/residue'
      if (kcenat .gt. 0) then
         write (6, '(/,a,a,/)')
     $      ' Central atom is ',atname(kcenat)
      endif
      if (lintyp .eq. 1) then
       call XYZCLOSE(IXYZIN)
      else
       close (unit=iun,status='keep')
      endif
      return
      end
c
c
c
      subroutine bonds
c     ================
c
c  Calculate bonds & set up connectivity list
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
c
c
      integer i, j, ires, jbond
      real ri, rj, r, d2
      logical newbnd
c
c
      nbond = 0
      do 1, j = 1,nresid
       do 2, i = 1,maxbnd
        ldist(i,j) = .false.
 2     continue 
 1    continue 
c
      do 5, ires = 1,nresid
       do 10, j = 1, natom-1
c     
        if (latom(j,ires)) then
         rj = atmrad(lattyp(j))
c     
         do 20, i = j+1, natom
          if (latom(i,ires)) then
           ri = atmrad(lattyp(i))
           r = ri + rj
c     
           d2 = (xyz(1,i,ires)-xyz(1,j,ires))**2 +
     $        (xyz(2,i,ires)-xyz(2,j,ires))**2 +
     $        (xyz(3,i,ires)-xyz(3,j,ires))**2 
c     
           if (d2 .lt. r**2) then
c     Yes it's a bond
            if (ires .eq. 1) then
             newbnd = .true.
            else
             call fbond(i,j,jbond)
             if (jbond .eq. 0) then
              newbnd = .true.
             else
              newbnd = .false.
             endif
            endif
c     
            if (newbnd) then
             nbond = nbond+1
             if (nbond .gt. maxbnd) then
              call ccperr(1,'*** Too many bonds ***')
             endif
             jbond = nbond
             kbond(1,jbond) = i
             kbond(2,jbond) = j
            endif
            distrs(jbond,ires) = sqrt(d2)
            ldist(jbond,ires) = .true.
           endif
          endif
 20      continue 
        endif
 10    continue 
 5    continue 
c Now average the bond lengths
      do 100, jbond = 1,nbond
       distnc(jbond) = 0.0
       sddist(jbond) = 0.0
       ndist(jbond) = 0
       do 110, i = 1,nresid
        if (ldist(jbond,i)) then
         ndist(jbond) = ndist(jbond) + 1
         distnc(jbond) = distnc(jbond) + distrs(jbond,i)
         sddist(jbond) = sddist(jbond) + distrs(jbond,i)**2
        else
         distrs(jbond,i) = 0.0
        endif        
 110   continue
       if (ndist(jbond) .gt. 0) then
        distnc(jbond) = distnc(jbond)/ndist(jbond)
        if (ndist(jbond) .gt. 1) then
         sddist(jbond) =
     $      sqrt(sddist(jbond)/float(ndist(jbond))
     $      - distnc(jbond)**2)
        else
         sddist(jbond) = 0.0
        endif
       endif
 100  continue 
      return
      end
c
c
c
      subroutine angles
c     =================
c
c  Calculate bond angles
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
c
c Locals
      integer i,j,k,l,m,ires,jangle
      real    rtod
      logical newang
c
      real bangle
      external bangle
c
      nangle = 0
      rtod = 45./atan(1.)
      do 1, j = 1,nresid
       do 2, i = 1,maxang
        langl(i,j) = .false.
 2     continue 
 1    continue 
c
c Loop residues
      do 5, ires = 1,nresid
c     Loop each atom as candidate for central atom in bond angle
c     
       do 10, k = 1,natom
        if (latom(k,ires)) then
c     
c     Search bond list for a bond involving this atom
         do 20, i = 1, nbond-1
           if (kbond(1,i) .eq. k .or. kbond(2,i) .eq. k) then
c     Bond i involves atom k, search for a second one
            do 30, j = i+1, nbond
              if (kbond(1,j) .eq. k .or. kbond(2,j) .eq. k) then
c     Two bonds, i & j from atom k to another atom, this defines a new angle
               if (kbond(1,i) .eq. k) then
                l = kbond(2,i)
               else
                l = kbond(1,i)
               endif
               if (kbond(1,j) .eq. k) then
                m = kbond(2,j)
               else
                m = kbond(1,j)
               endif
               if (ires .eq. 1) then
                newang = .true.
               else
c     Try to find angle in previous list
                call fangle(l,k,m,jangle)
                if (jangle .eq. 0) then
                 newang = .true.
                else
                 newang = .false.
                endif
               endif
               if (newang) then
                nangle = nangle+1
                if (nangle .gt. maxang) then
                 write (6, '(/a,i8,a/)')
     $              ' **** Too many angles, maximum: ',
     $              maxang,' ****'
                 call ccperr(1,'*** Too many angles ***')
                endif
                jangle = nangle
                kangle(1,jangle) = l
                kangle(2,jangle) = k
                kangle(3,jangle) = m
               endif
               anglrs(jangle,ires) =
     $            rtod*bangle(
     $            xyz(1,l,ires), xyz(1,k,ires), xyz(1,m,ires))
               langl(jangle,ires) = .true.
              endif
 30         continue 
           endif
 20      continue 
        endif
 10    continue 
 5    continue 
c Average the angles
      do 100, jangle = 1,nangle
       angleb(jangle) = 0.0
       sdangl(jangle) = 0.0
       nangl(jangle) = 0.0
       do 110, i = 1,nresid
        if (langl(jangle,i)) then
         nangl(jangle) = nangl(jangle) + 1
         angleb(jangle) = angleb(jangle) + anglrs(jangle,i)
         sdangl(jangle) = sdangl(jangle) + anglrs(jangle,i)**2
        else
         anglrs(jangle,i) = 0.0
        endif
 110   continue 
       if (nangl(jangle) .gt. 0) then
        angleb(jangle) = angleb(jangle)/float(nangl(jangle))
        if (nangl(jangle) .gt. 1) then
         sdangl(jangle) =
     $      sqrt(sdangl(jangle)/float(nangl(jangle)) -
     $      angleb(jangle)**2)
        else
         sdangl(jangle) = 0.0
        endif
       endif
 100  continue 
      return
      end
c
c
c
      real function bangle(x1,x2,x3)
c     ==============================
c
c  Calculate bond angle between coordinates x1,x2,x3
c     ie angle at x2, answer in radians
c
      real x1(3), x2(3), x3(3)
c
c
      real v1(3), v2(3), v3(3), cp, sp
c
      real dot
      external dot
c
      call vdif(v1, x2, x1)
      call vdif(v2, x2, x3)
c
      cp = dot(v1,v2)
      call cross(v3,v1,v2)
      sp = sqrt(dot(v3,v3))
      bangle = atan2(sp,cp)
      return
      end
c
c
c    
      subroutine wrtout(iatm)
c     =======================
c
c  Write out dictionary stuff
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Arguments
      character*2 iatm(100)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
c
c
c
      if (loutyp .eq. 1) then
         call wrttnt
      elseif (loutyp .eq. 2) then
c Protin output
         call wrtprt
      elseif (loutyp .eq. 3) then
         call wrtpdb(iatm)
      endif
      return
      end
c
c
c
      subroutine wrttnt
c     =================
c
c TNT output
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
c
c
      integer i,idict,j,k,l,n, nd, nb
      character rjatn1*4, rjatn2*4, rjatn3*4, rjatpl(maxapl)*4,
     $   line*1024, plntyp*8
      real sddmean, sdbmean, sd, sdwt, bdiff
      integer lenstr
      external lenstr
c
c  sdwt is relative weight of individual sd compared to overall weight
c
      sdwt = 0.1
c
      if (residu .ne. ' ') rsname = residu
c
      idict = 0
      call ccpdpn(idict, 'DICT', 'NEW', 'F', 0,0)
c
      nb = 0
      nd = 0
      sddmean = 0.0
      sdbmean = 0.0
      bdiff = 0.0
c
c==> Bonds
      do 10, k = 1, nbond
c  Weighted mean sd
         if (ndist(k) .gt. 1) then
            sd = sdwt * sddist(k) + (1.0-sdwt) * sdbond
         else
            sd = sdbond
         endif
         call rjstfy(rjatn1, atname(kbond(1,k)))
         call rjstfy(rjatn2, atname(kbond(2,k)))
         if (nresid .gt. 1) then
            write (idict, 6001) rsname, distnc(k), sd,
     $         rjatn1, rjatn2, (distrs(k,i),i=1,nresid)
            
 6001       format('GEOMETRY ',a,' BOND ',2f9.3,3x,a,', ',a,' # ',
     $         20f7.3)
         else
c  If only one, don't write out individual components
            write (idict, 6011) rsname, distnc(k), sd,
     $         rjatn1, rjatn2
 6011       format('GEOMETRY ',a,' BOND ',2f9.3,3x,a,', ',a)
         endif
         nd = nd+1
         sddmean = sddmean + sd
 10   continue
c
c==> Angles
      do 20, k = 1, nangle
         if (nangl(k) .gt. 1) then
c  Weighted mean sd
            sd = sdwt * sdangl(k) + (1.0-sdwt) * sdbang
         else
            sd = sdbang
         endif
         call rjstfy(rjatn1, atname(kangle(1,k)))
         call rjstfy(rjatn2, atname(kangle(2,k)))
         call rjstfy(rjatn3, atname(kangle(3,k)))
         if (nresid .gt. 1) then
            write (idict, 6002) rsname, angleb(k), sd,
     $         rjatn1, rjatn2, rjatn3, (anglrs(k,i),i=1,nresid)
 6002       format('GEOMETRY ',a,' ANGLE ',2f9.1,3x,a,', ',a,', ',a,
     $         ' # ',20f7.1)
         else
            write (idict, 6012) rsname, angleb(k), sd,
     $         rjatn1, rjatn2, rjatn3
 6012       format('GEOMETRY ',a,' ANGLE ',2f9.1,3x,a,', ',a,', ',a)
         endif
         nb = nb+1
         sdbmean = sdbmean + sd
 20   continue
c
c==> Planes
      if (nplanes .gt. 0) then
         do 30, k = 1, nplanes
            do 31, j = 1, natpln(k)
               call rjstfy(rjatpl(j), atname(kplane(j,k)))
 31         continue
            if (natpln(k) .eq. 4) then
c  TRIGONAL is a special case of PLANE
               plntyp = 'TRIGONAL'
               n = 0
            else
               plntyp = 'PLANE'
               n = natpln(k)
            endif
            write (line, 6004) rsname, plntyp, n, sdplane,
     $         (rjatpl(j), j=1,natpln(k))
 6004       format('GEOMETRY ',a,1x,a,i5,f7.2,100(a,','))
c  Strip final ","
            l = lenstr(line)
            if (line(l:l) .eq. ',') then
               line(l:l) = ' '
               l = l-1
            endif
            write (idict, '(a)') line(1:l)
 30      continue 
      endif
c
c Chirals
      if (nchirals .gt. 0) then
         do 35, k = 1, nchirals
            do 36, j = 1, 4
               call rjstfy(rjatpl(j), atname(kchiral(j,k)))
 36         continue
            write (line, 6005) rsname, 
     $         (rjatpl(j), j=1,4)
 6005       format('GEOMETRY ',a,'CHIRAL  1  1 ',100(a,','))
c  Strip final ","
            l = lenstr(line)
            if (line(l:l) .eq. ',') then
               line(l:l) = ' '
               l = l-1
            endif
            write (idict, '(a)') line(1:l)
 35      continue 
      endif
c  Blanks for CHIRAL
c       write (idict, 6003) rsname, rsname
c 6003  format(
c     $    '#GEOMETRY ',a,' PLANE  <Natom>  0.02 a1, a2, a3, . . ',/,
c     $    '#GEOMETRY ',a,' TRIGONAL  0     0.02 a1, a2, a3, a4',/,
c     $    '#GEOMETRY ',a,' CHIRAL    1      1   a1, a2, a3, a4')
c
c==> Bfactors
      do 40, k = 1, nbond
         call rjstfy(rjatn1, atname(kbond(1,k)))
         call rjstfy(rjatn2, atname(kbond(2,k)))
         write (idict, 6010) rsname, bdiff, sdbfac,
     $      rjatn1, rjatn2

 6010    format('GEOMETRY ',a,' BCORR ',2f9.3,3x,a,', ',a)
 40   continue
c
      if (nb .gt. 0) then
         write (6, '(/,a,f10.4,/,a,f10.4,/)')
     $      ' Mean sd(bond length) = ', sddmean/nb,
     $      ' Mean sd(bond angle)  = ', sdbmean/nb
      endif
      close (unit=idict,status='keep')
      return
      end
c
c
c
      subroutine wrtprt
c     =================
c
c Protin output
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
c
c
      integer i,j,k, i1, i2, nt, nc,
     $   j1(10), j2(10), ikwt(10), ibwt(10), ihand, istat
      real centre(3), x(3), y(3)
      character line*256
      integer lenstr
      external lenstr
c
c
      if (residu .ne. ' ') rsname = residu
c
      idict = 0
      call ccpdpn(idict, 'DICT', 'NEW', 'F', 0,0)
c
c Coordinates of centre, to subtract from each
c  [defaults to 1st atom]
      k = max(kcenat, 1)
      call getxyz(kcenat, centre, istat)
      if (istat .ne. 0) then
         call ccperr(1,'**** No centre atom ****')
      endif
c
      write (idict, '(a,a,/,a)')
     $   '## Dictionary entry in Protin format: edit out sections ',
     $   'separated by "##"',
     $   '##  Coordinates'
      nc = 1
      write (idict, 6100) nc, numres, rsname, rscode
 6100 format (47x,i3,i5,10x,a3,1x,a1)
      nc = -1
      do 100, j = 1,natom
         call getxyz(j, x, istat)
         if (istat .ne. 0) then
            call ccperr(1,'**** No centre atom ****')
         endif
c  Shift coordinates to centre
         call vdif(y, x, centre)
         write (idict, 6101) y, iprtyp(lattyp(j)),
     $      j, nc, atname(j)
 6101    format (3f10.5,10x,i5,5x,2i5,15x,a4)
 100  continue
c Bonds
      write (idict, '(a)') '## Distances (bonds & angles)'
      nt = nbond+nangle
      write (idict, 6110) rsname, numres, nt
 6110 format (A4,6x,2i5)
      i2 = 0
 110  i1 = i2 + 1
c   do while (i1 .le. nt)
      if (i1 .le. nt) then
         i2 = min(nt, i1 + 7)
c
c Distances (bonds & angles), write out 8 / line
c  For each istance, write atom1, atom2, ikwt, ibwt
c     ikwt, ibwt  are 1 1  for bond between mainchain atoms
c                     2 2  for angle distance (1-3) between mainchain atoms
         k = 0
         do 111, i = i1, i2
            k = k+1
            if (i .le. nbond) then
c  Bond distance
               j1(k) = kbond(1,i)
               j2(k) = kbond(2,i)
               ikwt(k) = 1
               ibwt(k) = 1
            else
c  Angle distance
               j1(k) = kangle(1,i-nbond)
               j2(k) = kangle(3,i-nbond)
               ikwt(k) = 2
               ibwt(k) = 2
            endif
 111     continue
         write (idict, 6111) (j1(i), j2(i), ikwt(i), ibwt(i), i=1,k)
 6111    format (8(2i3,2i2))
         go to 110
      endif
c Planes
      if (nplanes .gt. 0) then
         write (idict, '(a)') '## Planes'
         do 200, k = 1, nplanes
c
            write (line, 6201) rsname, numres, natpln(k), natpln(k),
     $         (kplane(i,k), i =  1, natpln(k))
 6201       format (a4,i2,2i3,17i4)
c   Residues with more than one plane are flagged in iat(17) as
c "multiplane" groups
            if (nplanes .gt. 1) then
               write (line(77:), 6202) numres
 6202          format (i4)
            endif
            write (idict, '(a)') line(1:lenstr(line))
 200     continue
      endif
c Chirals
      if (nchirals .gt. 0) then
         write (idict, '(a)') '## Chirals (all intrinsic, ihand = 1)'
         ihand = 1
         do 300, k = 1, nchirals
            write (idict, 6301) rsname, numres, ihand,
     $         (kchiral(i,k), i = 1,4), (atmcrl(i,k), i = 1,4)
 6301       format (a4,2i3,4i5,'    ! ',4a4)
 300     continue
      endif
c
      if (nnonbond .gt. 0) then
         write (idict, '(a)') '## Nonbonded contacts, types 1 & 2'
         write (idict, 6401) rsname, numres, nnonbond
 6401    format (A4,6x,2i5)
         i2 = 0
 400     i1 = i2+1
c   do while (i1 .lt. nnonbond) 
         if (i1 .le. nnonbond) then
            i2 = min(nnonbond, i1+9)
            k = 0
            do 410, i = i1, i2
               k = k+1
               j1(k) = knonbn(1,i)
               j2(k) = knonbn(2,i)
               ikwt(k) = ktnonb(i)
 410        continue
            write (idict, 6402) (j1(i), j2(i), ikwt(i), i=1,k)
 6402       format (10(2I3,I2))
            go to 400
         endif
      endif
      close (unit=idict,status='keep')
      return
      end
c
c
c
      subroutine getxyz(j, x, istat)
c     =============================
c
c  Get coordinates of atom j, from first residue in which it is present
c
c On entry:
c   j        atom number required
c
c On exit:
c   x(3)     coordinates
c   istat    =  0  OK
c            = +1  not available (shouldn't happen)
c
       integer j, istat
      real x(3)
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
      integer i, k
c
      istat = 0
      do 1, k = 1, nresid
         if (latom(j, k)) then
            do 2, i = 1,3
               x(i) = xyz(i,j,k)
 2          continue
            return
         endif
 1    continue 
      istat = +1
      return
      end
c
c
c
      subroutine wrtpdb(iatm)
c     =======================
c
c PDB output
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Arguments
      character*2 iatm(100)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
C&&*&& include  control.fh
c== control.fh ====
c
c   sdbond    overall bond sd
c   sdbang    overall angle sd
c   sdplane   plane sd
c   sdbfac    bfactor standard deviations
c   loutyp    = 1  write tnt dictionary
c             = 2  write Protin dictionary
c             = 3  write PDB file
c   lintyp    = 1  input PDB file
c             = 2 input Prolsq file

      common /control/ sdbond, sdbang, sdplane, sdbfac,
     $   loutyp, lintyp
      integer loutyp, lintyp
      real sdbond, sdbang, sdplane, sdbfac
c<== control.fh
C&&*&& end_include  control.fh
C&&*&& include  radii.fh
c== radii.fh ====
c
c  Bonding radii & atomtypes
c
      integer maxtyp
      parameter (maxtyp=30)
c
      common /cradii/ atmtyp(maxtyp)
      character atmtyp*2
c
      common /rradii/ atmrad(maxtyp), natmtyp, iatmz(maxtyp),
     $   iprtyp(maxtyp)
      real atmrad
      integer natmtyp, iatmz, iprtyp
c
      save /cradii/, /rradii/
c
c<===
C&&*&& end_include  radii.fh
c
      integer i, ifail, ixyzout, iz, j
      real xx, yy, zz, q, b(6)
      character chnnam*1,id*4,atnm*4,resno*4,inscod*1,altcod*1,segid*4
c
      ifail = 0
      ixyzout = 0
      call XYZOPEN('XYZOUT','OUTPUT',' ',IXYZOUT,ifail)
      segid = ' '
      resno = ' '
      inscod = ' '
      altcod = ' '
      iz = 0
      q = 1.0
      b(1) = 20.0
      do 10 i=2,6
       b(i) = 0.0
   10 continue
      chnnam = ' '
      do 200, j = 1,natom
        atnm = atname(j)
        xx = xyz(1,j,1)
        yy = xyz(2,j,1)
        zz = xyz(3,j,1)
C--- ID is (element symbol (atomic number (atomic type (atom j))))
        id = iatm(iatmz(lattyp(j)))//'  '
        call XYZATOM(IXYZOUT,j,atnm,rsname,chnnam,numres,resno,
     $               inscod,altcod,segid,iz,id)
        call XYZCOORD(IXYZOUT,'O','U',xx,yy,zz,q,biso,b)
        call XYZADVANCE(IXYZOUT,0,0,*200,*200)
 200  continue
      call XYZCLOSE(IXYZOUT)
      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 fbond(i,j,jbond)
c     ===========================
c
c  Find bond between atoms i & j in bondlist
c
c On entry:
c     i,j     atom numbers
c
c On exit:
c     jbond   bond number, = 0 if not found
c
      integer i,j,jbond
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
c
      do 1, jbond = 1, nbond
       if (kbond(1,jbond) .eq. i) then
        if (kbond(2,jbond) .eq. j) then
         go to 10
        endif
       elseif (kbond(2,jbond) .eq. i) then
        if (kbond(1,jbond) .eq. j) then
         go to 10
        endif
       endif
 1    continue 
      jbond = 0
 10   return
      end
c
c
c     
      subroutine fangle(l,k,m,jangle)
c     ===============================
c
c  Find angle between atoms l.k.m in angle list
c
c On entry:
c     l,k,m     atom numbers
c
c On exit:
c     jangle    angle number, = 0 if not found
c
      integer l,k,m,jangle
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
      do 1, jangle = 1, nangle
       if (kangle(2,jangle) .eq. k) then
        if (kangle(1,jangle) .eq. l) then
         if (kangle(3,jangle) .eq. m) then
          go to 10
         endif
        elseif (kangle(3,jangle) .eq. l) then
         if (kangle(1,jangle) .eq. m) then
          go to 10
         endif
        endif
       endif
 1    continue
      jangle = 0
 10   return
      end
c
c
c
      subroutine chirals(istat)
c     ========================
c
c  Get chiral pointers kchiral
c
c On exit:
c     istat      =  0  OK
c                = +1  atom in chirals not found in group
c                = +2  duplicate atoms
c
      integer istat
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
      integer i, j, k, l
c
      istat = 0
      if (nchirals .le. 0) return
c Loop chirals
      do 10, k = 1, nchirals
c  Loop atoms
         do 20, j = 1, 4
c   Find this atom
            do 30, i = 1, natom
               if (atmcrl(j,k) .eq. atname(i)) then
c  Set pointer
                  kchiral(j,k) = i
                  go to 21
               endif
 30         continue
            write (6, '(/a,a,a,i4,a/)')
     $         ' **** Atom ', atmcrl(j,k), ' of chiral ', k,
     $         ' not found in group ***'
            istat = +1
            return
c  Check for duplicates
 21         if (j .gt. 1) then
               do 22, i = 1, j-1
                  if (kchiral(i,k) .eq. kchiral(j,k)) then
                     write (6, '(/a,a,a,i4,a/)')
     $                  ' *** Duplicate atoms ',atmcrl(j,k),
     $                  ' in chiral ',k,' ***'
                     istat = +2
                     go to 10
                  endif
 22            continue
            endif
 20      continue
c
c Pointers to chiral k defined
c End chiral loop
 10   continue
      return
      end
c
c
c
      subroutine planes(istat)
c     ========================
c
c  Get plane pointers kplane
c
c On exit:
c     istat      =  0  OK
c                = +1  atom in plane not found in group
c                = +2  duplicate atoms
c                = +3  some plane in some residue could not be fitted
c
      integer istat
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
      integer i, j, k, l, iat, istat1
      real xyzpln(3, maxapl), wtpln(maxapl), plane_normal(3), cofg(3),
     $   thickness(3), radii(3), radius, rmat(3,3)
c
      istat = 0
      if (nplanes .le. 0) return
c Loop planes
      do 10, k = 1, nplanes
c  Loop atoms
         do 20, j = 1, natpln(k)
c   Find this atom
            do 30, i = 1, natom
               if (atmpln(j,k) .eq. atname(i)) then
c  Set pointer
                  kplane(j,k) = i
                  go to 21
               endif
 30         continue
            write (6, '(/a,a,a,i4,a/)')
     $         ' **** Atom ', atmpln(j,k), ' of plane ', k,
     $         ' not found in group ***'
            istat = +1
            return
c  Check for duplicates
 21         if (j .gt. 1) then
               do 22, i = 1, j-1
                  if (kplane(i,k) .eq. kplane(j,k)) then
                     write (6, '(/a,a,a,i4,a/)')
     $                  ' *** Duplicate atoms ',atmpln(j,k),
     $                  ' in plane ',k,' ***'
                     istat = +2
                     go to 10
                  endif
 22            continue
            endif
 20      continue
c
c Pointers to plane k defined
c   For each residue, fit plane & get deviation
         do 40, l = 1,nresid
c     Get coordinates, set unit weight (or 0 is absent)
            do 41, j = 1, natpln(k)
               iat = kplane(j,k)
               do 42, i=1,3
                  xyzpln(i,j) = xyz(i,iat,l)
 42            continue 
               if (latom(iat,l)) then
                  wtpln(j) = 1.0
               else
                  wtpln(j) = 0.0
               endif
 41         continue
c     Fit plane
            call fit_plane
     $         (xyzpln, wtpln, natpln(k),
     $         plane_normal, cofg, rmat, thickness, radii, radius,
     $         istat1)
            if (istat1 .ne. 0) then
               write (6, '(/a,i4,a,i4,a/)')
     $            ' *** Can''t calculate plane,',j,
     $            ' from residue', l,', too few atoms ***'
               istat = +3
               devpln(k,l) = 0.0
            else
c     rms deviation is "thickness" along z
               devpln(k,l) = thickness(3)
               write (6,'(a,i4,a,f8.3,a,i4,a)')
     $            ' RMS deviation from plane', k, ' is ',
     $            devpln(k,l),', ', natpln(k),' atoms'
            endif
c  End residue loop
 40      continue
c End plane loop
 10   continue
      return
      end
c
c
c
      subroutine fit_plane
     $   (xyz, weight, nat,
     $   plane_normal, cofg, rmat, thickness, radii, radius, istat)
c     =============================================================
c
c  Fit atoms to plane, from principle axes of inertia matrix
c
c On entry:
c   xyz(3,nat)     atom coordinates
c   weight(nat)    weight for fit
c   nat            number of atoms
c
c On exit:
c   plane_normal(3)  plane normal vector (unit)
c                  sign of normal is defined to be approximately parallel
c                  to the vector (xyz2-xyz1) x (xyz3-xyz1), where
c                  xyz1, xyz2, xyz3 are the coordinates of atoms 1,2,3
c   cofg(3)        centre of gravity of plane
c   rmat(3,3)      rotation matrix to rotate into plane coordinates
c                  (plane normal along z)
c   thickness(3)   "thickness" = weighted rms deviation along principle axes
c   radii(3)       "radii" = weighted rms deviation in plane perpendicular
c                  to principle axes
c   radius         rms radius
c   istat          = 0 OK, = -1 failure (too few atoms)
c
      integer nat, istat
      real xyz(3,nat), weight(nat), plane_normal(3), cofg(3),
     $   thickness(3), radii(3), radius, rmat(3,3)
c
      integer i,j,k,n
      real wt, dik, djk, a(3,3), v1(3), v2(3), vn(3), resid(3),
     $   vmat(3,3), det
c
      real dot
      external dot
c
      istat = 0
      n = 0
c  Weighted centre of gravity cofg
      wt = 0.0
      call ccpzr(cofg, 3)
c
      do 1, k = 1, nat
         if (weight(k) .gt. 0.0) then
            n = n+1
            wt = wt + weight(k)
            do 2, i = 1,3
               cofg(i) = cofg(i) + weight(k) * xyz(i,k)
 2          continue
         endif
 1    continue
      if (n .lt. 3) then
         istat = -1
         return
      endif
      do 3, i = 1,3
         cofg(i) = cofg(i)/wt
 3    continue 
c
c---
      call ccpzr(a, 9)
c  Inertia matrix 
      wt = 0.0
      do 10, k = 1,nat
         if (weight(k) .gt. 0.0) then
            wt = wt + weight(k)
            do 11, j = 1,3
               djk = xyz(j,k) - cofg(j)
               do 12, i = 1,3
                  dik = xyz(i,k) - cofg(i)
                  a(i,j) = a(i,j) + weight(k) * dik * djk
 12            continue 
 11         continue 
         endif
 10   continue 
c
      call eign3(a, vmat, resid)
c resid      eigenvalues of residual matrix, in descending order
c            = sums of weighted squared residuals along principle directions
c vmat(i,j)  eigenvector matrix
c            column vmat(.,3) is the plane normal
c
c  Get RMS residuals = thickness
      do 20, i = 1,3
         thickness(i) = sign(sqrt(abs(resid(i)/wt)), resid(i))
 20   continue
c Radii
      radii(1) = sqrt(thickness(2)**2 + thickness(3)**2)
      radii(2) = sqrt(thickness(3)**2 + thickness(1)**2)
      radii(3) = sqrt(thickness(1)**2 + thickness(2)**2)
      radius   = (radii(1)*radii(2)*radii(3))**(1./3.)
c
c Sign of plane_normal defined by approximate normal along vector
c (xyz2-xyz1) x (xyz3-xyz1),
c where  xyz1, xyz2, xyz3 are the coordinates of atoms 1,2,3
c  v1 = xyz2 - xyz1
      call vdif(v1, xyz(1,2), xyz(1,1))
c  v2 = xyz3 - xyz1
      call vdif(v2, xyz(1,3), xyz(1,1))
c  vn = v1 x v2
      call cross(vn, v1, v2)
c
      if (dot(vn, vmat(1,3)) .lt. 0.0) then
c Invert normal to make parallel to vn
         do 21, i = 1,3
            vmat(i,3) = -vmat(i,3)
 21      continue 
      endif
c Invert vmat(.,2) if necessary to make RH set
      call cross(v1, vmat(1,1), vmat(1,2))
      if (dot(v1, vmat(1,3)) .lt. 0.0) then
         do 22, i = 1,3
            vmat(i,2) = -vmat(i,2)
 22      continue 
      endif
c Normal
      call ccpmvr(plane_normal, vmat(1,3), 3)
c Transformation matrix is inverse of vmat
      call minv3(rmat, vmat, det)
c
      return
      end
C
C - - - - - - - - - - - -
      SUBROUTINE EIGN3(AMAT,VECS,VALS)
      real AMAT(3,3),VECS(3,3),VALS(3)
     $   ,WORK(5,3) 
C ONLY LOWER TRIANGLE OF AMAT(I,J) (IE,I.GE.J)   NEED BE PROVIDED 
C RETURNS EIG-VALS IN DESC ORDER IN VALS; 
C CORRESP  VECS AS COLS OF VECS 
C RETURNS AMAT AS SYMMETRIC FULL MATRIX 
C REQUIRES LIBRARY=HARWELL.FORTLIB FOR EA06C
c
      integer n, i1p, i1, ih, i, j
      real v
c
      CALL EA06C(AMAT,VALS,VECS,3,3,3,WORK) 
C RETURNS VALS & VECS IN ARBITRARY ORDER & SCREWS UP 
C UPPER TRIANGLE OF AMAT 
C ORDER VALS,VECS, THEN SYMMETRIZE AMAT
      N=3
      I1P=1
 10   I1=I1P
      I1P=I1P+1
      IF(I1P.GT.N) GO TO 60
      IH=I1
      V=VALS(IH)
      DO 20 I=I1P,N
	 IF(VALS(I).LE.V) GO TO 20
	 IH=I
	 V=VALS(IH)
 20   CONTINUE
      IF(IH.EQ.I1) GO TO 10
      V=VALS(I1)
      VALS(I1)=VALS(IH)
      VALS(IH)=V
      DO 30 J=1,N
	 V=VECS(J,I1)
	 VECS(J,I1)=VECS(J,IH)
 30   VECS(J,IH)=V
      GO TO 10

 60   CONTINUE
      DO 70 I=1,N
         DO 70 J=1,I
            AMAT(J,I)=AMAT(I,J)
 70   CONTINUE 
      RETURN
      END
c
c
c
      subroutine contac(istat)
c     ========================
c
c  Get atom pointers for contact list
c
c On exit:
c  istat =  0  OK
c        = +1  some atom name not found
c
      integer istat
c
C&&*&& include  cmolec.fh
c== cmolec.fh ====
c
c  Molecule stores
c
      integer maxatm, maxres
      parameter (maxatm=1000, maxres=100)
      integer maxbnd, maxang, maxapl, maxpln, maxcrl, maxnnb
      parameter (maxbnd=3*maxatm)
      parameter (maxang=2*maxbnd)
      parameter (maxpln=20, maxapl=100)
      parameter (maxcrl=30)
      parameter (maxnnb=1000)
c
c Residue:
c     rsname        residue name
c     rscode        1-character name (for Protin)
c     numres        residue type number (for Protin)
      common /nmolec/
c Atoms:
c     natom         number of atoms (maximum number in any residue)
c     xyz           coordinates
c     latom         .true. is this atom present in this residue
c     lattyp        atom type number
c     nresid        number of copies of residue
c     atname        atom name
c     cenatm        central atom name
c     kcenat        central atom number
     $     xyz(3,maxatm,maxres),
c Bonds:
c     nbond         number of bonds
c     distrs        distances in each residue
c     ldist         .true. if this distance defined for this residue
c     kbond         atom pointers for this bond
c     distnc        distance (averaged)
c     sddist        standard deviation of distance
c     ndist         number of observations of this distance
     $     distrs(maxbnd,maxres), 
     $     distnc(maxbnd), sddist(maxbnd),
c Angles:
c     nangle        number of angles
c     anglrs        angles in each residue
c     langl         .true. if angle defined in this residue
c     angleb        angle (averaged)
c     sdangl        standard deviation of angle
c     nangl         number of observations of this angle
c     kangle        atom pointers for this angle
     $     anglrs(maxang,maxres),
     $     angleb(maxang), sdangl(maxang),
c Planes
c     nplanes       number of planes
c     kplane        atom pointers for this plane
c     natpln        number of atoms in this plane
c     atmpln        atom names for this plane
c     devpln        rms deviation from plane, for this residue
     $     devpln(maxpln,maxres),
c
c Chirals
c     nchirals      number of chirals
c     kchiral       atom pointers for this chiral group (central atom 1st)
c     atmcrl        atom names for this chiral group
c   
c Nonbonded contacts
c     nnonbond      number of non-bonded pairs
c     atnonb        atomname pairs
c     knonbn        atom pointers to pair
c     ktnonb        type, = 1 single torsion angle, =2 multiple torsion angle
c-----
c integers
     $     lattyp(maxatm),
     $     natom, kcenat, nbond, kbond(2,maxbnd),
     $     nangle, kangle(3,maxang),
     $     nresid, ndist(maxbnd), nangl(maxang), numres,
     $     nplanes, kplane(maxapl,maxpln), natpln(maxpln),
     $     nchirals, kchiral(4,maxcrl),
     $     nnonbond, knonbn(2,maxnnb), ktnonb(maxnnb),
c logicals
     $     latom(maxatm,maxres), ldist(maxbnd,maxres),
     $     langl(maxang,maxres)
c
      real xyz, distrs, distnc, sddist, anglrs, angleb, sdangl, devpln
      integer lattyp, natom, kcenat, nbond, kbond, nangle, kangle,
     $     nresid, ndist, nangl, numres,
     $     nplanes, kplane, natpln, nchirals, kchiral,
     $     nnonbond, knonbn, ktnonb
      logical latom, ldist, langl
c
      common /cmolec/ atname(maxatm), rsname, residu, attype(maxatm),
     $   rscode, atmpln, atmcrl, cenatm, atnonb
      character atname*4, rsname*4, residu*4, attype*2, rscode*1,
     $   atmpln(maxapl,maxpln)*4, atmcrl(4,maxcrl)*4, cenatm*4,
     $   atnonb(2,maxnnb)*4
c
      save /nmolec/, /cmolec/
c
c<== cmolec.fh
C&&*&& end_include  cmolec.fh
c
      integer i, j, k
      logical fail
c
      istat = 0
      if (nnonbond .le. 0) return
c
c  Find atom pointers
      fail = .false.
      do 1, k = 1, nnonbond
         do 2, j = 1,2
            knonbn(j,k) = 0
            do 3, i = 1, natom
               if (atnonb(j,k) .eq. atname(i)) then
                  knonbn(j,k) = i
                  go to 2
               endif
 3          continue
            write (6, '(/a,a,a/)')
     $         ' *** Non-bonded Atom ', atnonb(j,k),
     $         ' not recognized ***'
            fail = .true.
c
 2       continue
 1    continue
c
      if (fail) then
         istat = +1
         return
      endif
c
      end