C
C
C     This code is distributed under the terms and conditions of the
C     CCP4 licence agreement as `Part ii)' software.  See the conditions
C     in the CCP4 manual for a copyright statement.
C
C********************************************************************
C
C  This program finds solvent extended accessible area of atoms
C  in a Brookhaven coordinate file and writes an extra column to the file.
C  Program generates symmetry related molecules if Smode = 'I'.
C
C  Updated to use Keyworded input
C
C  There are four possible modes of operation, keywords are
C
C  1) MODE ALL | NOHOH | HOH | HOHALL
C
C  ALL    Calculate accessible area for all atoms, including waters
C         if present in file. Water atoms are treated as solvent when
C         calculating accessible area.
C  NOHOH  All waters (residue type HOH or WAT) are ignored, and are not
C         written to the output file.
C  HOH    The accessible area will only be calculated for waters (HOH or
C         WAT), treating other waters as solvent. Only waters will be
C         written to output file.
C  HOHALL As HOH, but waters are treated as protein, and consequently
C         more waters will have low solvent accessiblity.
C
C  2) SMODE IMOL | OFF
C
C  IMOL   generate symmetry related molecules from coordinates in
C         Brookhaven file before calculating accessible areas.
C  OFF    symmetry related molecules are not generated
C
C  3) SYMMETRY <spgrp.name> | <spgrp.no.> | <operations>
C
C         Standard keyword ie. input can be spacegroup name, spacegroup
C         number, or a list of operations in the format -X,-Y,-Z etc (in
C         which case operations are separated by '*').
C         SYMMETRY is compulsory if SMODE = IMOL, and optional otherwise.
C
C  4) EXCLUDE <res1> <res2> ... <res3>
C
C         EXCLUDE is followed by a list of residues; atoms belonging to
C         residues with residue name res1, res2 etc (3 character format)
C         will be ignored in the area calculations, and will NOT be
C         written to the output Brookhaven file, up to 30 residues.
C         NOTE: the EXCLUDE keyword can appear any number of times.
C
C  New keywords added:
C
C  5) PROBE <x>
C
C         Sets the solvent radius (used as a "probe") to be equal to <x>
C         angstroms (where x is a real number). The default value is 1.4
C         (Added 29/5/98)
C
C  6) ATOM <name> <number> <radius>
C
C         Allows user to add/overwrite atom type information used in the
C         area calculation.
C
C  7) PNTDEN <no_of_points>
C
C         Sets the density of points used in the area calculation equal to
C         <no_of_points> per square angstrom. The precision of the
C         calculation is then the reciprocal of this value.
C
C  October 1998:
C
C  AREAIMOL is now combined with the programs DIFFAREA, WATERAREA and RESAREA
C  The area calculations are automatically piped into WATERAREA or RESAREA for
C  analysis. Additionally, differences in area can be analysed through
C  DIFFAREA, controlled by a new keyword DIFFMODE:
C
C  8) DIFFMODE OFF | IMOL | WATER | COMPARE | CAVITY
C
C     This must be the first keyword, if it appears at all, and controls what
C     input the program expects and how it handle the data.
C
C         OFF: equivalent to the original AREAIMOL program.
C                This is the default.
C
C         IMOL: works on one input file, and finds the difference in area
C                when two different sets of symmetry operators are applied.
C                No SMODE keyword is needed, if one SYMMETRY keyword is
C                given then the second set of operators effectively consists
C                only of the identity, otherwise two SYMMETRY keywords are
C                needed to specify both sets of operators.
C                All other keywords have their usual function.
C
C         WATER: works on one input pdb file, and finds the difference in area
C                between the HOH and HOHALL modes. No MODE keyword is needed,
C                other keywords function as described.
C
C         COMPARE: requires two input files. Does area calculation for each
C                and then gets the differences for only those atoms which are
C                common to both files. All other keywords have their usual
C                function.
C
C         CAVITY: works on a single input file. Does multiple ASA calculations
C                with increasing probe radius to search for cavitis and surface
C                features.
C                WARNING THIS IS A DEVELOPMENTAL FEATURE AND SHOULD NOT BE USED!
C
C  9) VERBOSE
C
C     Switch on additional output:
C         List of recognised atom types
C         SCALE matrices (for conversion btwn orth. and frac. coords)
C         Symmetry operation matrices
C         
C
C  10) TRANS NONE | 1 | 2 | BOTH
C
C     Short for "translations mode"; this switches on the generation of
C     lattice vectors during the generation of symmetry-related atoms as
C     follows:
C
C         NONE: no translations are applied
C         1   : translations are applied when symmetry-related atoms are
C               generated during the first of two area calculations only
C         2   : as "1" but for the second area calculation only
C         BOTH: translations are applied during both of the area calcs.
C
C     The TRANS keyword is only compatible with DIFFMODE IMOL or SMODE IMOL.
C
C  11) CAV_PROBE <step> <max>
C
C     Set extra parameters for DIFFMODE CAVITIES
C     <step> is the probe radius increment which will be applied at each
C            iteration
C     <max>  is the maximum probe radius to be used - the calculation
C            terminates once this is exceeded
C     Both values are in Angstroms
C
C     WARNING THIS IS A DEVELOPMENTAL FEATURE AND SHOULD NOT BE USED!
C
C  12) MATCHUP ALL | NOCOORDS
C
C     Switch for DIFFMODE COMPARE, to control the comparision criteria
C     used when doing comparision of XYZIN and XYZIN2:
C
C         ALL     : (default) atoms are compared only if atom name, residue
C                   name/number, chain id and atomic coordinates are the
C                   same between both files.
C         NOCOORDS: use the same criteria as ALL, except for the atomic
C                   coordinates.
C
C  13) OUTPUT
C
C     Switch on output of areas (or area differences) to a pseudo pdb file.
C
C  February 2001:
C
C  TRANS mode modified so that up to +/-2 lattice vector space is searched
C  when generating symmetry related molecules (cf BIGSEARCH mode in CONTACT).
C  This is because searching only +/-1 lattice vectors may not find all the
C  crystallographic molecules in contact with the original molecule.
C
C  7 Nov 2001:
C
C  Added mask array SELECT_CALC to allow internal atom selection.
C  Added DIFFMODE CAVITIES and CAV_PROBE for developmental cavity search.
C  WARNING THESE ARE DEVELOPMENTAL AND SHOULD NOT BE USED!
C
C  20 Dec 2001:
C
C  Allow large probe radii (required some redimensioning).
C  Added MATCHUP keyword.
C  Added calculation of contact areas.
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER MAXTYPE
      PARAMETER (MAXTYPE=20)
      INTEGER MAXRES
      PARAMETER (MAXRES=30)
C     ..
C     ..
C     .. Local Scalars ..
      REAL HOHRAD,PNTDEN,PROBE_STEP,PROBE_MAX
      INTEGER NSYM1,NSYM2,NREJ,NATOMTYP
      INTEGER NREADIN1,NREADIN2,IDIFF,ITRANS
      LOGICAL ALL,HOH,HOHALL,NOHOH
      LOGICAL SOFF,SIMOL,SMAT,SGEN,VERB,NOTRANS,LOUTPUT,LCLOSED
      LOGICAL USE_COORDS
C     ..
C     .. Local Arrays ..
      REAL TR1(4,4,12),TR2(4,4,12)
      REAL X1(3,MAXATOM),X2(3,MAXATOM),XYZLIM(2,3),CELL(6)
      REAL RADIUS(0:MAXTYPE),BRKMAT(4,4),IBRKMAT(4,4)
      INTEGER NPOINT1(MAXATOM),NPOINT2(MAXATOM)
      INTEGER ITYPE1(MAXATOM),ITYPE2(MAXATOM)
      INTEGER ATOMN(0:MAXTYPE),NOINTYP1(0:MAXTYPE),NOINTYP2(0:MAXTYPE)
      INTEGER SELECT_CALC(MAXATOM),SURFACE(MAXATOM)
      CHARACTER ATYPE(0:MAXTYPE)*2
      CHARACTER ATOMSTORE1(MAXATOM)*10,ATOMSTORE2(MAXATOM)*10
      CHARACTER REJTYP(MAXRES)*3
      LOGICAL*1 RESTYP1(MAXATOM),RESTYP2(MAXATOM)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPFYP,CCPRCS,CCPERR,ATOMSETUP,KEYINPUT,SETRADII,READIN,
     +         SETFLAGS,GETAREAS,MATCHUP,DIFFAREA,RESAREA,WATERAREA,
     +         WRITEOUT
C
C          ******
      CALL CCPFYP
C          ******
      call ccp4h_init()  
      call ccp4h_header('AREAIMOL','areaimol',2)
      CALL CCPRCS(6, 'AREAIMOL', '$Date: 2002/07/30 16:19:25 $')  

      call ccp4h_summary_beg()
      call ccp4h_toc_beg()
      call ccp4h_toc_ent('Command Input','#command')
      call ccp4h_toc_ent('Mode','#mode')
      call ccp4h_toc_ent('Data From Input Co-ordinate File','#data')
      call ccp4h_toc_ent('Total Area Calculation','#area')
      call ccp4h_toc_end()
      call ccp4h_rule()
      call ccp4h_header('Command Input','command',3)
      call ccp4h_header('   Main Keywords','main',3)
      call ccp4h_link('DIFFMODE   ','areaimol.html#diffmode')
      call ccp4h_link('MODE ','areaimol.html#mode')
      call ccp4h_link('SMODE ','areaimol.html#smode')
      call ccp4h_link('SYMMETRY  ','areaimol.html#symmetry')
      call ccp4h_link('TRANS','areaimol.html#trans')
      call ccp4h_header('   Secondary Keywords','secondary',3)
      call ccp4h_link('ATOM   ','areaimol.html#atom')
      call ccp4h_link('EXCLUDE ','areaimol.html#exclude')
      call ccp4h_link('MATCHUP ','areaimol.html#matchup')
      call ccp4h_link('PNTDEN ','areaimol.html#exclude')
      call ccp4h_link('PROBE  ','areaimol.html#probe')
      call ccp4h_header('   Auxiliary Keywords','auxiliary',3)
      call ccp4h_link('VERBOSE  ','areaimol.html#verbose')  
      call ccp4h_link('OUTPUT  ','areaimol.html#output')
      call ccp4h_link('END','areaimol.html#end')
      call ccp4h_summary_end()
C
C --- Print a warning about this new version
C
      WRITE (6,fmt='(a,/,a,/,a,//,a,//)')
     +    ' NOTE: this is an updated version of AREAIMOL which also',
     +    ' includes the functions of the old RESAREA, WATERAREA and',
     +    ' DIFFAREA programs',
     +    ' See program documentation for details.'
C
C --- Set up internal database of recognised atom types and radii
C
      CALL ATOMSETUP(ATYPE,RADIUS,ATOMN,NATOMTYP,MAXTYPE)
C
C --- Get keyworded input
C
      CALL KEYINPUT(ATYPE,RADIUS,ATOMN,NATOMTYP,MAXTYPE,ALL,NOHOH,HOH,
     +              HOHALL,SOFF,SIMOL,HOHRAD,PNTDEN,IDIFF,NSYM1,TR1,
     +              NSYM2,TR2,NREJ,REJTYP,VERB,ITRANS,LOUTPUT,LCLOSED,
     +              PROBE_STEP,PROBE_MAX,USE_COORDS)
C
C --- Setup atom selection
C
      CALL SETUP_SELECTION(MAXATOM,SELECT_CALC)
C
C---- Increment radii by solvent radius
C
      CALL SETRADII(MAXTYPE,NATOMTYP,HOHRAD,RADIUS)
C
C --- Read in data from XYZIN
C
      CALL READIN('XYZIN',NREADIN1,ATOMSTORE1,RADIUS,SOFF,NOHOH,NREJ,
     +             REJTYP,X1,CELL,XYZLIM,RESTYP1,ITYPE1,ATOMN,ATYPE,
     +             NATOMTYP,NOINTYP1,BRKMAT,IBRKMAT,VERB)
C  
      call ccp4h_summary_beg()
      call ccp4h_header('Total Area Calculation','area',3)
      call ccp4h_summary_end()
C --- Set the MODE and SMODE flags
C
      CALL SETFLAGS(1,IDIFF,NSYM1,SOFF,ALL,NOHOH,HOH,HOHALL,SGEN,SMAT,
     +              ITRANS,NOTRANS)
C
C --- Obtain the areas and output the analysis
C
      CALL GETAREAS(NREADIN1,X1,ATOMSTORE1,MAXATOM,XYZLIM,HOHRAD,
     +              SOFF,SGEN,SMAT,BRKMAT,IBRKMAT,NSYM1,TR1,ITYPE1,
     +              RESTYP1,ALL,HOHALL,HOH,NOHOH,RADIUS,NATOMTYP,PNTDEN,
     +              NOINTYP1,MAXTYPE,NPOINT1,NOTRANS,VERB,LCLOSED,
     +              SELECT_CALC,SURFACE,PROBE_STEP,PROBE_MAX)
C
C
C --- In DIFFMODE COMPARE there are two files; read in XYZIN2
C
      IF (IDIFF .EQ. 3) THEN
C
        CALL READIN('XYZIN2',NREADIN2,ATOMSTORE2,RADIUS,SOFF,NOHOH,NREJ,
     +              REJTYP,X2,CELL,XYZLIM,RESTYP2,ITYPE2,ATOMN,ATYPE,
     +              NATOMTYP,NOINTYP2,BRKMAT,IBRKMAT,VERB)
        CALL SETFLAGS(2,IDIFF,NSYM2,SOFF,ALL,NOHOH,HOH,HOHALL,SGEN,SMAT,
     +                ITRANS,NOTRANS)
        CALL GETAREAS(NREADIN2,X2,ATOMSTORE2,MAXATOM,XYZLIM,HOHRAD,
     +                SOFF,SGEN,SMAT,BRKMAT,IBRKMAT,NSYM2,TR2,ITYPE2,
     +                RESTYP2,ALL,HOHALL,HOH,NOHOH,RADIUS,NATOMTYP,
     +                PNTDEN,NOINTYP2,MAXTYPE,NPOINT2,NOTRANS,VERB,
     +                LCLOSED,SELECT_CALC,SURFACE,PROBE_STEP,PROBE_MAX)
C
      ELSE IF (IDIFF .EQ. 2) THEN
C
        CALL SETFLAGS(2,IDIFF,NSYM1,SOFF,ALL,NOHOH,HOH,HOHALL,SGEN,SMAT,
     +                ITRANS,NOTRANS)
        CALL GETAREAS(NREADIN1,X1,ATOMSTORE1,MAXATOM,XYZLIM,HOHRAD,
     +                SOFF,SGEN,SMAT,BRKMAT,IBRKMAT,NSYM1,TR1,ITYPE1,
     +                RESTYP1,ALL,HOHALL,HOH,NOHOH,RADIUS,NATOMTYP,
     +                PNTDEN,NOINTYP1,MAXTYPE,NPOINT2,NOTRANS,VERB,
     +                LCLOSED,SELECT_CALC,SURFACE,PROBE_STEP,PROBE_MAX)
C
      ELSE IF (IDIFF .EQ. 1) THEN
C
        CALL SETFLAGS(2,IDIFF,NSYM2,SOFF,ALL,NOHOH,HOH,HOHALL,SGEN,SMAT,
     +                ITRANS,NOTRANS)
        CALL GETAREAS(NREADIN1,X1,ATOMSTORE1,MAXATOM,XYZLIM,HOHRAD,
     +                SOFF,SGEN,SMAT,BRKMAT,IBRKMAT,NSYM2,TR2,ITYPE1,
     +                RESTYP1,ALL,HOHALL,HOH,NOHOH,RADIUS,NATOMTYP,
     +                PNTDEN,NOINTYP1,MAXTYPE,NPOINT2,NOTRANS,VERB,
     +                LCLOSED,SELECT_CALC,SURFACE,PROBE_STEP,PROBE_MAX)
      END IF
C
C --- If DIFFMODE is not OFF or CAVITY then do the area difference analysis
C
      IF (IDIFF .NE. 0 .AND. IDIFF .NE. 4) THEN
C
C     ... and if DIFFMODE COMPARE then match up the atoms
C
        IF (IDIFF .EQ. 3) THEN
          CALL MATCHUP(NREADIN1,NREADIN2,X1,ATOMSTORE1,X2,ATOMSTORE2,
     +                 NPOINT1,NPOINT2,ITYPE1,ITYPE2,RESTYP1,RESTYP2,
     +                 USE_COORDS)
        END IF
C
        CALL DIFFAREA(NREADIN1,NPOINT1,NPOINT2,PNTDEN,ATYPE,NATOMTYP,
     +                ITYPE1,RESTYP1,IDIFF)
C
C --- Do the final analysis on the differenced areas
C
        IF (IDIFF .EQ. 2) THEN
          CALL WATERAREA(NREADIN1,ATOMSTORE1,NPOINT1,PNTDEN,
     +                   HOH,HOHALL,0)
        ELSE
          IF (ALL .OR. NOHOH) THEN
            CALL RESAREA(NREADIN1,ATOMSTORE1,NPOINT1,PNTDEN,RESTYP1,
     +                   NOHOH,0)
          ELSE
            CALL WATERAREA(NREADIN1,ATOMSTORE1,NPOINT1,PNTDEN,
     +                     HOH,HOHALL,0)
          END IF
        END IF
C
      END IF
C
C---- Write out a psuedo-pdb file if required
C
      IF (LOUTPUT) CALL WRITEOUT('XYZOUT',CELL,X1,ATOMSTORE1,ATYPE,
     +                           ITYPE1,NPOINT1,PNTDEN,MAXTYPE,MAXATOM,
     +                           NREADIN1,IDIFF,HOH,HOHALL,RESTYP1)
C
C          ************************************************
      CALL CCPERR(0,' Normal termination')
C          ************************************************
C
      END
C     ....

C
C     ================================================================
      SUBROUTINE GETAREAS(NREADIN,X,ATOMSTORE,MAXATOM,XYZLIM,HOHRAD,
     +                    SOFF,SGEN,SMAT,BRKMAT,IBRKMAT,NSYM,TR,ITYPE,
     +                    RESTYP,ALL,HOHALL,HOH,NOHOH,RADIUS,NATOMTYP,
     +                    PNTDEN,NOINTYP,MAXTYPE,NPOINT,NOTRANS,VERB,
     +                    LCLOSED,SELECT_CALC,SURFACE,PROBE_STEP,
     +                    PROBE_MAX)
C     ================================================================
C
C---- Given the atomic coordinates etc return the array NPOINT
C     containing the no. of surface points for each atom.
C
C---- On entry: NREADIN - number of real atoms
C               X contains atomic coordinates (in A, orthog axies)
C               ATOMSTORE - seq. id from pdb file
C               MAXATOM - maximum number of atoms
C               XYZLIM - limits of grid for searching
C               HOHRAD - probe radius
C               SOFF - flag: if set then no symmetry gen. is peformed
C               SGEN - flag: if set then symm. generation is allowed
C               SMAT - flag: if set then symm. matrix gen. is allowed
C               BRKMAT - matrix derived from CRYST and SCALE cards
C               IBRKMAT - inverse of BRKMAT
C               NSYM - no of symmetry ops. in TR
C               TR - array holding symmetry operation matrices 
C               ITYPE - atom type (ie number) for each atom
C               RESTYP - whether atom is protein or water
C               ALL - flag: if set then MODE ALL
C               HOHALL - flag: if set then MODE HOHALL
C               HOH - flag: if set then MODE HOH
C               NOHOH - flag: if set then MODE NOHOH
C               RADIUS contains radius for each atom type
C               NATOMTYP - number of recognised atom types
C               PNTDEN - point density (no. of dot per sqr.A)
C               NOINTYP - no of atoms of each type
C               MAXTYPE - maximum no of allowed atom types
C               VERB - flag for verbose output
C               LCLOSED - flag to select closed surface search
C               SELECT_CLAC - array: flags atom selection for asa calc.
C               SURFACE - array: one element for each atom (used in
C                         isolated surface searching & analysis
C               PROBE_STEP,PROBE_MAX : increment for probe radius (and
C                         maximum value) in cavity analysis
C
C---- on exit:  NPOINT - no of surface points on each atom
C
C---- Local variables: 
C               BINSIZE - physical length of cubic cell edge in NET
C               NET - grid with atoms listed in each spatical cell
C               MINX,MAXX,MINY,MAXY,MINZ,MAXZ - limits of NET
C               TRVEC - array holding translation vectors
C               S - array holding transformed symmetry matrices
C               NATOM - total (ie real + symm.gener'ed) no. of atoms
C
C
      IMPLICIT NONE
C
C   ..Parameters..
      INTEGER MAXNET
      PARAMETER (MAXNET=5000000)
      INTEGER MAXBIN
      PARAMETER (MAXBIN=1000)
C
C   ..Arguments ..
      INTEGER MAXATOM,MAXTYPE,NREADIN,NSYM,NATOMTYP
      INTEGER ITYPE(MAXATOM),NOINTYP(0:MAXTYPE),NPOINT(MAXATOM)
      INTEGER SELECT_CALC(MAXATOM),SURFACE(MAXATOM)
      REAL PNTDEN,HOHRAD,PROBE_STEP,PROBE_MAX
      REAL X(3,MAXATOM),XYZLIM(2,3),BRKMAT(4,4),IBRKMAT(4,4)
      REAL TR(4,4,12),RADIUS(0:MAXTYPE)
      CHARACTER ATOMSTORE(MAXATOM)*10
      LOGICAL SOFF,SMAT,SGEN,ALL,HOHALL,HOH,NOHOH,NOTRANS,VERB
      LOGICAL LCLOSED,LCONTACT
      LOGICAL*1 RESTYP(MAXATOM)
C
C   ..Local Scalars ..
      INTEGER I,J,NSAREA,NINSRF
      INTEGER NATOM,MAXX,MINX,MAXY,MINY,MAXZ,MINZ,NSURFACES
      REAL    BINSIZE,PROBE_DECR,SAREA
      CHARACTER RESNAM*3,CHNAME*1,SEQID*6,STRING*10
C
C   ..Local Arrays ..
      INTEGER NET(MAXNET)
      REAL S(4,4,12),TRVEC(3,125),COM(3)
C
C   ..External subroutines ..
      EXTERNAL BRICKSETUP,SYMMAT,SYMGEN,BINATOMS,AREA,RESAREA,WATERAREA
C
C---- Initialise NATOM - this records the total no of atoms
C     i.e. real + symmetry generated
C
      NATOM = NREADIN
C
      LCONTACT = .TRUE.
C
C---- Set up bricking array for searching
C
      CALL BRICKSETUP(BINSIZE,HOHRAD,XYZLIM,MINX,MAXX,MINY,MAXY,MINZ,
     +                MAXZ,MAXNET,MAXBIN,NET,VERB)
C
C---- Generate symmetry-related atoms
C
      IF (.NOT. SOFF) THEN
         IF (SGEN) THEN
           CALL SYMMAT(BRKMAT,IBRKMAT,NSYM,S,TR,TRVEC,VERB)
           CALL SYMGEN(NATOM,X,XYZLIM,BINSIZE,ITYPE,
     +                 RESTYP,NSYM,S,TRVEC,NOTRANS,VERB)
         END IF
      END IF
C
C---- Put atoms in bins
C
      CALL BINATOMS(X,NATOM,NET,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,
     +          BINSIZE,RESTYP,HOHALL,VERB)
C
C---- Calculate area for each atom
C
      CALL AREA(X,ITYPE,RADIUS,NATOMTYP,PNTDEN,NPOINT,NET,MINX,
     +          MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,MAXATOM,BINSIZE,RESTYP,
     +          NOINTYP,NREADIN,HOH,HOHALL,SELECT_CALC)
C
C---- Analysis of total accessible area by residue (MODEs ALL, NOHOH)
C
      IF (ALL .OR. NOHOH) THEN
        CALL RESAREA(NREADIN,ATOMSTORE,NPOINT,PNTDEN,RESTYP,NOHOH,1)
C
C---- Analysis of contact areas by residue
        IF (LCONTACT) THEN
          CALL CONTACTAREA(NREADIN,ATOMSTORE,NPOINT,PNTDEN,RESTYP,
     +         NOHOH,HOHRAD,ITYPE,NATOMTYP,RADIUS,1)
        END IF
C
C---- Search for closed surfaces (clusters of atoms with non-zero ASA)
C
        IF (LCLOSED) THEN
          CALL SEARCHISOLATED(SELECT_CALC,NREADIN,X,RADIUS,HOHRAD,
     +         NPOINT,ITYPE,NET,BINSIZE,MAXATOM,NATOMTYP,MINX,MAXX,MINY,
     +         MAXY,MINZ,MAXZ,MAXBIN,PNTDEN,NSURFACES,SURFACE,1)


C
C---- Write out results
C
          WRITE(6,FMT='(/,A,/,A,/)')
     +         ' RESULTS OF ISOLATED-SURFACE SEARCH',
     +         ' =================================='
          IF (NSURFACES.EQ.1) THEN
            call ccp4h_summary_beg()
            WRITE(6,FMT='(A,/)')
     +           ' All surface atoms belong to a single surface'
            call ccp4h_summary_end()
          ELSE
            call ccp4h_summary_beg()
            WRITE(6,FMT='(I3,A)') NSURFACES,
     +           ' isolated surfaces were found'
            call ccp4h_summary_end()
            WRITE(6,FMT='(/,A,/)') ' List of surfaces follows:'
          END IF
C
C---- Loop over surfaces and calculate number of atoms, centre of mass
C     and total ASA
C
          DO I=1,NSURFACES
            NSAREA = 0
            NINSRF = 0
            COM(1) = 0.0
            COM(2) = 0.0
            COM(3) = 0.0
            DO J=1,NREADIN
              IF (SURFACE(J).EQ.I) THEN
                NSAREA = NSAREA + NPOINT(J)
                NINSRF = NINSRF + 1
                COM(1) = COM(1) + X(1,J)
                COM(2) = COM(2) + X(2,J)
                COM(3) = COM(3) + X(3,J)
              END IF
            END DO
            IF (NINSRF.LE.0) CALL CCPERR (1,
     +           ' Fatal error in s/r GETAREAS: empty surface!')
            COM(1) = COM(1)/REAL(NINSRF)
            COM(2) = COM(2)/REAL(NINSRF)
            COM(3) = COM(3)/REAL(NINSRF)
            SAREA  = NSAREA/PNTDEN
            WRITE(6,FMT='(A,I3)')' SURFACE ',I
            WRITE(6,FMT='(A,I8)')'    No of surface atoms : ',NINSRF
            WRITE(6,FMT='(A,F10.2)')'    Total area  : ',SAREA
            WRITE(6,FMT='(A,3F8.2,/)')'    COM (x,y,z) : ',
     +           (COM(J),J=1,3)
            IF (NINSRF.LE.10) THEN
              DO J=1,NREADIN
                IF (SURFACE(J).EQ.I) THEN
                  STRING = ATOMSTORE(J)
                  RESNAM = STRING(1:3)
                  CHNAME = STRING(5:5)
                  SEQID  = STRING(5:10)
                  WRITE(6,FMT='(5X,A,1X,A,2X,A,2X,3F8.2)')
     +                 RESNAM,CHNAME,SEQID,X(1,J),X(2,J),X(3,J)
                END IF
              END DO
              WRITE(6,*)
            END IF
          END DO
C
C---- Further analysis of surfaces
C
C---- PROBE_STEP is the increase in probe radius for each step
C     PROBE_MAX  is the maximum probe radis
          IF (PROBE_STEP.GT.0.0) THEN
C
            DO I=1,NSURFACES
              DO J=1,NREADIN
                IF (SURFACE(J).EQ.I) THEN
                  SELECT_CALC(J) = 1
                ELSE
                  SELECT_CALC(J) = -1
                ENDIF
              END DO
              WRITE(6,*) 'Starting PROBE_FATTENING for surface ',I
              CALL PROBE_FATTENING(I,NREADIN,SELECT_CALC,HOHRAD,
     +             PROBE_STEP,PROBE_MAX,X,ITYPE,RADIUS,NATOMTYP,PNTDEN,
     +             NPOINT,NET,MAXNET,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,
     +             MAXBIN,BINSIZE,XYZLIM,RESTYP,NOINTYP,HOH,HOHALL)
C---- Reset the expanded VdW radii the calculations
              PROBE_DECR = HOHRAD - PROBE_MAX
              CALL SETRADII(NATOMTYP,NATOMTYP,PROBE_DECR,RADIUS)
            END DO
          END IF
C
        END IF
C
C---- Depending on the MODE, do analysis of waters
C
      ELSE IF (HOH .OR. HOHALL) THEN
        CALL WATERAREA(NREADIN,ATOMSTORE,NPOINT,PNTDEN,HOH,HOHALL,1)
C     
      END IF
C
      RETURN
      END
C
C
C     ================================================================
      SUBROUTINE PROBE_FATTENING(SURF_NO,NATOM,SELECT_CALC,PROBE_RAD,
     +     PROBE_STEP,PROBE_MAX,X,ITYPE,RADIUS,NATOMTYP,PNTDEN,NPOINT,
     +     NET,MAXNET,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,BINSIZE,
     +     XYZLIM,RESTYP,NOINTYP,HOH,HOHALL)
C     ================================================================
C
C Perform the atomic fattening procedure for selected atoms
C Atomic fattening based on the algorithm in Kleywegt and Jones
C Acta Cryst D (1994) 50 pp178-185
C In this case the size of the probe sphere is increased in steps
C and the areas calculated each time
C
C Arguments
C     SURF_NO     : identifier for the surface being analysed
C     NATOM       : no of atoms (used for sizing arrays)
C     SELECT_CALC : atom selection array, one element per atom
C     PROBE_RAD   : initial probe sphere radius
C     PROBE_STEP  : increment for probe radius for each round of
C                   calculations
C     PROBE_MAX   : on entry, upper limit on probe radius
C                   on exit, actual maximum value of probe radius
C     
C
      IMPLICIT NONE
C
C     ..Parameters
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER MAXRUNS
      PARAMETER (MAXRUNS=30)
C
C     ..Arguments..
      INTEGER SURF_NO,NATOM,NATOMTYP
      INTEGER MAXX,MINX,MAXY,MINY,MAXZ,MINZ,MAXBIN,MAXNET
      INTEGER SELECT_CALC(NATOM),ITYPE(NATOM),NPOINT(NATOM)
      INTEGER NET(0:MAXBIN,MINX:MAXX,MINY:MAXY,MINZ:MAXZ)
      INTEGER NOINTYP(0:NATOMTYP)
      REAL    PROBE_RAD,PROBE_STEP,PROBE_MAX,PNTDEN,BINSIZE
      REAL    X(3,NATOM),RADIUS(0:NATOMTYP),XYZLIM(2,3)
      LOGICAL HOHALL,HOH
      LOGICAL*1 RESTYP(NATOM)
C
C     ..Local scalars..
      INTEGER I,J,N_BURIED,N_CLUSTER,N_BURIED_CLUSTER
      INTEGER IRUN,N_IN_CLUSTER
      REAL    PROBE,TOTAL_ASA
      LOGICAL VERB
C
C     ..Local arrays..
      INTEGER CLUSTER(MAXATOM),RUN_N_BURIED(MAXRUNS)
      INTEGER RUN_N_CLUSTER(MAXRUNS),RUN_N_BURIED_CLUSTER(MAXRUNS)
      REAL    RUN_PROBE(MAXRUNS),RUN_TOTAL_ASA(MAXRUNS)
C
C---- Initial probe size
C
      PROBE = PROBE_RAD
C
C---- Run index
C
      IRUN = 0
C
C---- Loop over probe radius
C
 10   CONTINUE
      IRUN = IRUN + 1
C
C---- Re-brick and re-bin the atoms
      VERB = .FALSE.
      CALL BRICKSETUP(BINSIZE,PROBE_RAD,XYZLIM,MINX,MAXX,MINY,MAXY,MINZ,
     +     MAXZ,MAXNET,MAXBIN,NET,VERB)
      CALL BINATOMS(X,NATOM,NET,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,
     +     BINSIZE,RESTYP,HOHALL,VERB)

C
C---- Get the areas on all the selected atoms
C
      CALL AREA(X,ITYPE,RADIUS,NATOMTYP,PNTDEN,NPOINT,NET,MINX,
     +          MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,NATOM,BINSIZE,RESTYP,
     +          NOINTYP,NATOM,HOH,HOHALL,SELECT_CALC)
C
C---- Do analysis of the results
C
C---- Sum total ASA and count number of buried atoms
      TOTAL_ASA = 0.0
      N_BURIED  = 0
      DO I = 1,NATOM
        TOTAL_ASA = TOTAL_ASA + NPOINT(I)/PNTDEN
        IF (SELECT_CALC(I).EQ.1 .AND. NPOINT(I).EQ.0) THEN
          N_BURIED = N_BURIED + 1
        END IF
      END DO
C
C---- Clustering
C
C     1. Cluster atoms with positive ASA
      CALL SEARCHISOLATED(SELECT_CALC,NATOM,X,RADIUS,PROBE,NPOINT,
     +         ITYPE,NET,BINSIZE,MAXATOM,NATOMTYP,MINX,MAXX,MINY,
     +         MAXY,MINZ,MAXZ,MAXBIN,PNTDEN,N_CLUSTER,CLUSTER,1)
C
C     2. Cluster atoms with zero ASA
      CALL SEARCHISOLATED(SELECT_CALC,NATOM,X,RADIUS,PROBE,NPOINT,
     +         ITYPE,NET,BINSIZE,MAXATOM,NATOMTYP,MINX,MAXX,MINY,
     +         MAXY,MINZ,MAXZ,MAXBIN,PNTDEN,N_BURIED_CLUSTER,CLUSTER,2)
C
C---- Write out some information
C
      WRITE(6,FMT='(1X,A,I2,A,I2,A)')
     +     'Run ',IRUN,': there are now ',N_BURIED_CLUSTER,
     +     ' clusters of buried atoms'
      DO I = 1,N_BURIED_CLUSTER
        WRITE(6,FMT='(5X,A,I2)') 'Cluster ',I
        N_IN_CLUSTER = 0
        DO J = 1,NATOM
          IF (CLUSTER(J).EQ.I) N_IN_CLUSTER = N_IN_CLUSTER + 1
        END DO
        WRITE(6,FMT='(5X,A,I5,A)') 'This cluster contains ',
     +       N_IN_CLUSTER,' buried atoms'
      END DO
C
C---- Store the results of this run to be displayed at the end
C
      RUN_PROBE(IRUN) = PROBE
      RUN_TOTAL_ASA(IRUN) = TOTAL_ASA
      RUN_N_BURIED(IRUN) = N_BURIED
      RUN_N_CLUSTER(IRUN) = N_CLUSTER
      RUN_N_BURIED_CLUSTER(IRUN) = N_BURIED_CLUSTER
C
C---- Increment probe radius and generate the new
C     values for the expanded VdW radii
      CALL SETRADII(NATOMTYP,NATOMTYP,PROBE_STEP,RADIUS)
      PROBE = PROBE + PROBE_STEP
C
C---- Repeat if required
      IF (PROBE.LE.PROBE_MAX) THEN
        IF (IRUN.LE.MAXRUNS) GO TO 10
      END IF
C
C---- Finished - Write out the loggraph
C
C---- Write the header
C
      WRITE(6,6000) SURF_NO
 6000 FORMAT (/,
     +' $TABLE:Analysis of closed surface no.',I3,':',/,
     +' $GRAPHS',/,
     +' :Total ASA as a function of probe radius:A:1,2:',/,
     +' :Number of buried atoms as function of probe radius:A:1,3:',/,
     +' :Number of clusters of buried and exposed atoms:A:1,4,5:',/,
     +' $$',/,
     +' Probe_Radius Total_ASA N_buried N_cluster N_cluster_buried',/,
     +' $$ $$ ')
C
C---- Write the results in the loggraph
C
      DO I = 1,IRUN
        WRITE(6,FMT='(F10.2,1X,F10.2,1X,I10,1X,I10,1X,I10)')
     +       RUN_PROBE(I),RUN_TOTAL_ASA(I),RUN_N_BURIED(I),
     +       RUN_N_CLUSTER(I),RUN_N_BURIED_CLUSTER(I)
      END DO
C
C---- End the loggraph and return
C
      WRITE(6,*) '$$'
      WRITE(6,*)
      PROBE_MAX = PROBE
      RETURN
      END
C
C
C     ===============================================================
      SUBROUTINE SETFLAGS(CYCLE,IDIFF,NSYM,SOFF,ALL,NOHOH,HOH,HOHALL,
     +                    SGEN,SMAT,ITRANS,NOTRANS)
C     ===============================================================
C
C---- ON input, CYCLE  = 1 or 2 (corresponding to 1st or 2nd
C                        calculation)
C               IDIFF  = 0,1,2,3 or 4 (DIFFMODE)
C               NSYM   = no. of symmetry operators
C               ITRANS = 0,1,2 or 3 (TRANSlations mode)
C
C---- On output, ALL/NOHOH/HOH/HOHALL = set correctly in DIFFMODE
C                                       WATERS (otherwise unchanged)
C                SOFF = set correctly for DIFFMODE IMOL (otherwise
C                                       unchanged)
C                SMAT = allow symmetry matrices to be built if TRUE
C                SGEN = allow symmetry-related atoms to be generated
C                       if TRUE
C                NOTRANS = set TRUE or FALSE according to ITRANS,
C                       (and possibly CYCLE) - see code below.
C
C     NB: even if SMAT and SGEN are set TRUE, the matrices/symm-rel.
C     atoms will only be generated if SOFF is FALSE. See s/r GETAREAS.     
C
      IMPLICIT NONE
C
C   ..Arguments
      INTEGER IDIFF,CYCLE,NSYM,ITRANS
      LOGICAL SOFF,ALL,NOHOH,HOH,HOHALL,SGEN,SMAT,NOTRANS
C
C   ..External subroutines
      EXTERNAL CCPERR
C
C
C---  Set up flags for diffarea
C
C---- DIFFMODE 0 = AREAIMOL
C     DIFFMODE 4 = CAVITY
C
      IF (IDIFF .EQ. 0 .OR. IDIFF .EQ. 4) THEN
          SGEN = .TRUE.
          SMAT = .TRUE.
          WRITE(6,FMT=1000)
 1000     FORMAT(/,1X,'===================== STARTING AREA CALCULATION',
     +               ' =====================')
C
C---- DIFFMODE 1 = IMOL
C
      ELSE IF (IDIFF .EQ. 1) THEN
        SOFF = .FALSE.
        IF (NSYM .EQ. 0) SOFF = .TRUE.
        SGEN = .TRUE.
        SMAT = .TRUE.
        WRITE(6,FMT=1010) CYCLE
 1010   FORMAT(/,1X,'===================== STARTING AREA CALCULATION ',
     +           I1,' =====================')
C
C---- DIFFMODE 2 = WATER
C
      ELSE IF (IDIFF .EQ. 2) THEN
        IF (CYCLE. EQ. 1) THEN
          ALL    = .FALSE.
          NOHOH  = .FALSE.
          HOH    = .TRUE.
          HOHALL = .FALSE.
          SGEN   = .TRUE.
          SMAT   = .TRUE.
        ELSE IF (CYCLE .EQ. 2) THEN
          ALL    = .FALSE.
          NOHOH  = .FALSE.
          HOH    = .FALSE.
          HOHALL = .TRUE.
          SGEN   = .TRUE.
          SMAT   = .FALSE.
        ELSE
          CALL CCPERR(1, 'BAD VALUE OF CYCLE IN S/R SETFLAGS')
        END IF
        WRITE(6,FMT=1020) CYCLE
 1020   FORMAT(/,1X,'===================== STARTING AREA CALCULATION ',
     +           I1,' =====================')
C
C---- DIFFMODE 3 = COMPARE (same as AREAIMOL above)
C
      ELSE IF (IDIFF .EQ. 3) THEN
        SGEN = .TRUE.
        SMAT = .TRUE.
        WRITE(6,FMT=1030) CYCLE
 1030   FORMAT(/,1X,'===================== STARTING AREA CALCULATION ',
     +           I1,' =====================')
      END IF
C
C---- Set NOTRANS flag
C
C---- DEFAULT: the translations are switched off
C
      NOTRANS = .TRUE.
C
C---- Now check to see if we have to switch them on
C
C     In DIFFMODE IMOL:
C
      IF (IDIFF .EQ. 1) THEN
        IF ((ITRANS .EQ. CYCLE).OR.(ITRANS .EQ. 3)) NOTRANS = .FALSE.
      ELSE IF (ITRANS .EQ. 1) THEN
        NOTRANS = .FALSE.
      END IF
C
C---- Print out summary information
C
C     MODE:
C
      IF (ALL)    WRITE(6,FMT='(/,A)')' MODE set to ALL'
      IF (HOHALL) WRITE(6,FMT='(/,A)')' MODE set to HOHALL'
      IF (NOHOH)  WRITE(6,FMT='(/,A)')' MODE set to NOHOH'
      IF (HOH   ) WRITE(6,FMT='(/,A)')' MODE set to HOH'
C
C     SYMMETRY GENERATED ATOMS:
C
      IF (SOFF) THEN
        WRITE(6,FMT='(/,2A)')' No symmetry related atoms',
     +                     ' will be generated'
      ELSE
        IF (SGEN) WRITE(6,FMT='(/,2A)')' Symmetry related atoms will',
     +                           ' be generated for this calculation'
        IF (NOTRANS) THEN
          WRITE(6,FMT='(A)')' using symmetry operations only'
        ELSE
          WRITE(6,FMT='(2A)')' using symmetry operations plus lattice',
     +                          ' translations'
        END IF
      END IF
C
C
      RETURN
C
      END
C
C
C     ==================================================================
      SUBROUTINE READIN(NAME,NREAD,ATOMSTORE,RADIUS,
     +                  SOFF,NOHOH,NREJ,REJTYP,X,CELL,XYZLIM,
     +                  RESTYP,ITYPE,ATOMN,ATYPE,NATOMTYP,NOINTYP,
     +                  BRKMAT,IBRKMAT,VERB)
C     ==================================================================
C
C  Subroutine to extract information from PDB file to be used in the
C  area calculations.
C
C  Converted to use rwbrook routines Oct 1999
C
C  Arguments:
C
C  INPUT:
C
C    NAME (character)   logical name of file to read in, eg 'XYZIN'
C    SOFF (logical)     flag = TRUE if no symmetry-related atoms are
C                       required
C    NOHOH (logical)    flag = TRUE if MODE is NOHOH (ie. reject waters)
C    NREJ (integer)     number of residue names specified for rejection
C    REJTYP (character*3 array)   list of residue names to be rejected
C    ATOMN (integer array)   list of recognised atomic numbers
C    ATYPE (character*2 array)    list of recognised atomic names
C    NATOMTYP (integer)      number of recognised atom types
C    VERB (logical)     switch on extra output
C    RADIUS (real)      VdW radii for each atom type
C
C  OUTPUT:
C
C    NREAD (integer)         number of atoms read in
C    ATOMSTORE (character*10 array)  chunk of ATOM record containing seq.id.
C    X (real array)          x, y and z position of atom
C    XYZLIM (real array)     highest & lowest values of x, y and z co-ords.
C    RESTYP (logical array)  flag = TRUE if residue is not water
C    ITYPE (integer array)   atom type
C    BRKMAT (real array)     transformation matrix from SCALE cards
C    IBRKMAT (real array)    inverse of BRKMAT above.
C    CELL (real array)       cell parameters from input file
C
      IMPLICIT NONE
C
C   ..
C   ..Parameters
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER MAXTYPE
      PARAMETER (MAXTYPE=20)
      INTEGER MAXRES
      PARAMETER (MAXRES=30)
C   ..
C   ..Scalar arguments
      INTEGER NREAD, NREJ, NATOMTYP
      LOGICAL SOFF, NOHOH, VERB
      CHARACTER NAME*(*)
C   ..
C   ..Array arguments
      INTEGER ITYPE(MAXATOM), ATOMN(0:MAXTYPE), NOINTYP(0:MAXTYPE)
      REAL X(3,MAXATOM),XYZLIM(2,3),BRKMAT(4,4),IBRKMAT(4,4),CELL(6)
      REAL RADIUS(0:MAXTYPE)
      LOGICAL*1 RESTYP(MAXATOM)
      CHARACTER ATOMSTORE(MAXATOM)*10, REJTYP(MAXRES)*3
      CHARACTER ATYPE(0:MAXTYPE)*2
C   ..
C   ..Local scalars
      INTEGER I, IFAIL, IXYZIN, IRESN, ISER, IZ
      INTEGER NOCCREJ, NEXREJ, NWATREJ, NHREJ
      REAL OCC, VOL, X1, Y1, Z1, BISO
      LOGICAL LOCC, WATER, LANISO
      CHARACTER STRING*80, ID*4, INSCOD*1, ALTCOD*1, ATMCHK*20
      CHARACTER ATNAM*4, RESNAM*3 ,CHNAM*1, RESNO*4, SEGID*5
C
C   ..Local arrays
      INTEGER NOINRES(MAXRES)
      REAL U(6)
C
C   ..External Subroutines ..
      EXTERNAL XYZINIT,XYZOPEN,RBCELL,RBRORF,XYZADVANCE,XYZATOM
      EXTERNAL XYZCOORD,CCPERR,XYZCLOSE
C
C
C --- Initialise some arrays and variables
C
      ATMCHK = ' '
      NREAD = 0
C
      DO 10 I = 1,MAXATOM
        RESTYP(I) = .FALSE.
   10 CONTINUE
C
      DO 20 I = 0,MAXTYPE
        NOINTYP(I) = 0
   20 CONTINUE
C
      LOCC = .FALSE.
C
      NOCCREJ = 0
      NEXREJ  = 0
      NWATREJ = 0
      NHREJ   = 0
C
      DO 25 I = 1,MAXRES
        NOINRES(I) = 0
 25   CONTINUE
C
   30 DO 40 I = 1,3
        XYZLIM(1,I) = 100000.0
        XYZLIM(2,I) = -100000.0
   40 CONTINUE
C
      NREAD = 0
C
C -- Open the input file
C
C    nb: there is another occurance of xyzinit in
C    s/r writeout
      CALL XYZINIT
      IFAIL  = 0
      IXYZIN = 0
C     *******************************************
      CALL XYZOPEN(NAME,'INPUT',' ',IXYZIN,IFAIL)
C     *******************************************
C
C -- Recover the cell and orthogonalising matrices
C
      CALL RBCELL(CELL,VOL)
      IBRKMAT(1,1) = 0.0
      BRKMAT(1,1)  = 0.0
      CALL RBRORF(IBRKMAT,BRKMAT)
C
C -- Read in the atom data
C
      call ccp4h_summary_beg()
      call ccp4h_header('Data From Input Coordinate File','data',3) 
      WRITE (6,FMT='(//,2A)')' READING IN ATOM DATA FROM INPUT',
     +                       ' COORDINATE FILE'
C
 110  CALL XYZADVANCE(IXYZIN,0,0,*201,*201)
      CALL XYZATOM(IXYZIN,ISER,ATNAM,RESNAM,CHNAM,IRESN,RESNO,INSCOD,
     +             ALTCOD,SEGID,IZ,ID)
C
      WATER = (RESNAM(1:3).EQ.'HOH' .OR. RESNAM(1:3).EQ.'WAT')
      IF (WATER .AND. NOHOH) THEN
        NWATREJ = NWATREJ + 1
        GO TO 110
      END IF
C
C -- Reset aniso U factor before reading in
      DO I = 1,6
        U(I) = 0.0
      END DO
C
C -- Get atomic coordinates etc
      CALL XYZCOORD(IXYZIN,'O','U',X1,Y1,Z1,OCC,BISO,U)
C
C -- Check that we don't treat an ANISOU card as an atom
C    If 
      LANISO = .TRUE.
      DO I = 1,6
        LANISO = (U(I).NE.0.0 .AND. LANISO)
      END DO
      IF (LANISO) THEN
        IF (X1.EQ.0.0 .AND. Y1.EQ.0.0 .AND. Z1.EQ. 0.0) GO TO 110
      END IF
C
C -- Check the occupancy - if it's less than 1.0 then set the flag
      IF ((.NOT.LOCC) .AND. (OCC .LT. 1.0)) LOCC = .TRUE.
C
C -- If it's zero then reject the atom
      IF (OCC .EQ. 0.0) THEN
        NOCCREJ = NOCCREJ + 1
        GOTO 110
      ENDIF
C
C -- Reject any EXCLUDEd residues
C
      IF (NREJ.NE.0) THEN
        DO 130 I = 1,NREJ
          IF (RESNAM.EQ.REJTYP(I)) THEN
            NOINRES(I) = NOINRES(I) + 1
            NEXREJ = NEXREJ + 1
            GO TO 110
          ENDIF
  130   CONTINUE
      END IF
C
C----  Reject hydrogens
      IF (ID(1:2).EQ.' H'.OR.IZ.EQ.1) THEN
        NHREJ = NHREJ + 1
        GO TO 110
      END IF
C
C -- Store the atom information and increment counter
C
      NREAD = NREAD + 1
      IF (NREAD .GT. MAXATOM)
     +          CALL CCPERR(1,'TOO MANY ATOMS: INCREASE MAXATOM')
      IF (.NOT.WATER) RESTYP(NREAD) = .TRUE.     
C
C -- Coordinates:
      X(1,NREAD) = X1
      X(2,NREAD) = Y1
      X(3,NREAD) = Z1
C
C -- Store residue name (1-3)
C          chainid (5)
C          residue sequence number (6-9)
C          insertion code (10)
      STRING = ' '
      WRITE(STRING(1:3),FMT='(A)')  RESNAM
      WRITE(STRING(5:5),FMT='(A)')  CHNAM
      WRITE(STRING(6:9),FMT='(I4)') IRESN
      WRITE(STRING(10:10),FMT='(A)')INSCOD
      ATOMSTORE(NREAD) = STRING(1:10)
C
C----  Find limits of molecule
C
      DO 140 I = 1,3
        XYZLIM(1,I) = MIN(X(I,NREAD),XYZLIM(1,I))
        XYZLIM(2,I) = MAX(X(I,NREAD),XYZLIM(2,I))
  140 CONTINUE
C
C----  Find atom type
C
      ITYPE(NREAD) = 0
C
C      First check atomic no.
      DO 145 I= 1,NATOMTYP
       IF(IZ .EQ. ATOMN(I))THEN
        ITYPE(NREAD) = I
        NOINTYP(I) = NOINTYP(I) + 1
        GOTO 110
       ENDIF
  145 CONTINUE
C
C      Use atom name to get atom type
  149 DO 150 I = 1,NATOMTYP
        IF (ID(1:2).EQ.ATYPE(I)) THEN
          ITYPE(NREAD) = I
          NOINTYP(I) = NOINTYP(I) + 1
          GO TO 110
        END IF
  150 CONTINUE
C
C ---  If ITYPE = 0 then atom type is unknown (default radius 1.8A)
C
      IF (ITYPE(NREAD).EQ.0) THEN
        WRITE (6,FMT=6026) STRING(1:10),ID(1:2)
 6026   FORMAT (' atom type found and not recognised = ',/1X,A,1X,A)
C
C          ********************************
        CALL CCPERR(2,
     +  'Unknown atom type: adding to list with default radius 1.8 A')
C          ********************************
C
        CALL ADDATOMTYPE(MAXTYPE,ATYPE,RADIUS,ATOMN,ATMCHK,
     +       ID(1:2),1.8,IZ,NATOMTYP)
C
        NOINTYP(NATOMTYP) = NOINTYP(NATOMTYP) + 1
        ITYPE(NREAD) = NATOMTYP
C
C       Also count number of new atom types
        NOINTYP(0) = NOINTYP(0) + 1
        GO TO 110
      ENDIF
C
  160 CONTINUE
C
C----  File read over - close the file
C
  201 CONTINUE
C
      CALL XYZCLOSE(IXYZIN)
C
C----  Do a bit of messaging
C
      WRITE (6,FMT=6027) NREAD
 6027 FORMAT (//' Read in ',I5,' atoms from pdb file',
     +          ' of which there were:',/)
C
      DO 200 I=1,MAXTYPE
        IF (NOINTYP(I).GT.0) THEN
          WRITE (6,FMT=6035) NOINTYP(I),ATYPE(I),ATOMN(I)
 6035     FORMAT (3X,I5,' atoms of type ',A2,' (number ',I2,')')
        ENDIF
  200 CONTINUE
C
      IF (NOINTYP(0).EQ.1) THEN
        WRITE (6,FMT=6040) NOINTYP(0)
        WRITE (6,FMT=6055)
      ELSE IF (NOINTYP(0).GT.0) THEN
        WRITE (6,FMT=6045) NOINTYP(0)
        WRITE (6,FMT=6055)
      ELSE
        WRITE (6,FMT=6050)
      ENDIF
 6040 FORMAT (/,I4,' atom type was unidentified and will be treated as',
     +         /,' having default radius 1.8A')
 6045 FORMAT (/,I4,' atom types were unidentified and will be treated',
     +         /,' as having default radius 1.8A')
 6050 FORMAT (/,' All atoms were identified',/)
 6055 FORMAT (' Use the ATOM keyword to explicitly define new atom',
     +         ' types',/)
C
      IF (NHREJ.GT.0) THEN
        WRITE (6,FMT=7110) NHREJ
      ENDIF
 7110 FORMAT (/,I4,' atoms were rejected as hydrogens',/)
C
      IF (LOCC) THEN
        WRITE (6,FMT=7050)
        IF (NOCCREJ.GT.0) THEN
          WRITE (6,FMT=7060) NOCCREJ
        ENDIF
      ENDIF
 7050 FORMAT (/,' One or more atoms had occupancy < 1.0 in',
     +          ' input file.',/)
 7060 FORMAT (I5,' atoms with occupancy = 0 have been excluded.',/)
C
      IF (NREJ.GT.0) THEN
C
        IF (NEXREJ.GT.0) THEN
          WRITE (6,FMT=7070) NEXREJ
C
          DO 310 I = 1,NREJ
            IF (NOINRES(I).GT.0) WRITE (6,FMT=7080) NOINRES(I),REJTYP(I)
 310      CONTINUE
C
        ELSE
          WRITE (6,FMT=7090)
        ENDIF
C
      ENDIF
C          
 7070 FORMAT (/,I5,' atoms were rejected on the basis of residue ',
     +             'type',/)
 7080 FORMAT (3X,I5,' atoms in residue type ',A3)
 7090 FORMAT (/,' No atoms were identified as belonging to the',
     +          ' EXCLUDED residues')
C
      IF (NWATREJ.GT.0) WRITE (6,FMT=7100) NWATREJ
 7100 FORMAT (//,I5,' atoms were rejected as being either WAT or HOH')
      call ccp4h_summary_end()
C
      RETURN
C
      END
C
C
C     =================================================================
      SUBROUTINE MATCHUP(NREADIN1,NREADIN2,X1,ATOMSTORE1,X2,ATOMSTORE2,
     +                   NPOINT1,NPOINT2,ITYPE1,ITYPE2,RESTYP1,RESTYP2,
     +                   USE_COORDS)
C     =================================================================
C
C---- Remove unwanted atoms and reorder so that the lists correspond
C
C     On entry: NREADIN* - no of real atoms in each list
C               X*       - x,y,z coordinates of each atom
C               ATOMSTORE* - seq. id. of each atom
C               NPOINT*   - no of surface points for each atom
C               ITYPE*    - atom type (ie atomic no) of each atom
C               RESTYP*   - TRUE if atom is protein (ie not WAT or HOH)
C
C               USE_COORDS - TRUE if coordinates are also to be used
C                            as a matchup criterion
C
C     On exit:  NREADIN1 - no of matching atoms
C               X1 - sorted/matched xyz coordinates
C               ATOMSTORE1 - sorted/matched seq. id.
C               NPOINT1 - sorted/matched no of surface points
C               ITYPE1  - sorted/matched atomic numbers
C               RESTYP1 - sorted/matched WAT/HOH flags
C
      IMPLICIT NONE
C
C   ..Parameters
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
C
C   ..Arguments
      INTEGER NREADIN1,NREADIN2
      INTEGER NPOINT1(MAXATOM),NPOINT2(MAXATOM)
      INTEGER ITYPE1(MAXATOM),ITYPE2(MAXATOM)
      REAL X1(3,MAXATOM),X2(3,MAXATOM)
      CHARACTER*10 ATOMSTORE1(MAXATOM),ATOMSTORE2(MAXATOM)
      LOGICAL*1 RESTYP1(MAXATOM),RESTYP2(MAXATOM)
      LOGICAL USE_COORDS
C
C   ..Local scalars
      INTEGER IC,JC,JSTART,NATOMS
      LOGICAL DEL
C   ..Local arrays
      INTEGER DELFLAG1(MAXATOM),NP1(MAXATOM),NP2(MAXATOM)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C
C---- Tell the user what's happening (since this procedure may take
C     some time)
C
      WRITE(6,FMT=9010)

      IF (USE_COORDS) THEN
        WRITE(6,FMT=9020)
      END IF
      WRITE(6,FMT=9030)
 9010 FORMAT(//,'Attempting to remove extra (ie unduplicated) atoms',
     +        /,'The following information will be used:'
     +        /,'     Atom type, residue name and number, chain id')
 9020 FORMAT('     Atomic coordinates')
 9030 FORMAT(/,'Only atoms with matching attributes will be considered',
     +       /,'in the comparison.',
     +       //,'Matching may take a while - please wait ....')
C
C---- Set up the flagging arrays
C
      DO 100 IC = 1,MAXATOM
        DELFLAG1(IC) = 0
 100  CONTINUE
C
C---- PART 1: Check the two lists for unique entries
C     ie. entries appearing in one but not the other
C
      JSTART = 1
C
C---- First list
C
C---- Loop over all entries in the first list
      DO 200 IC = 1,NREADIN1
C
C----   JC loops over all entries in the second list, starting
C       one on from the last matched item in that list
        JC  = JSTART
C
C----   Jump out if the second list has been exhausted
        IF (JC .GT. NREADIN2) GO TO 220
C
 210    CONTINUE
C
C----   Assume that the entry of interest in 1st list
C       matches the current entry in the 2nd list
        DEL = .FALSE.
C
C----   Check all the data and set flags.
C
C----   Use coordinates as a match criterion
        IF (USE_COORDS) THEN
          IF (X1(1,IC) .NE. X2(1,JC)) DEL = .TRUE.
          IF (X1(2,IC) .NE. X2(2,JC)) DEL = .TRUE.
          IF (X1(3,IC) .NE. X2(3,JC)) DEL = .TRUE.
        ENDIF
C
C-----  Match on atom type and residue information
        IF (ITYPE1(IC) .NE. ITYPE2(JC)) DEL = .TRUE.
        IF (ATOMSTORE1(IC) .NE. ATOMSTORE2(JC)) DEL = .TRUE.
C
C----   If DEL is TRUE then the two entries differ
C       Check the next entry in the 2nd list, or (if the
C       2nd list is exhausted) flag this entry as
C       unmatched
C
        IF (DEL) THEN
          JC = JC + 1
          IF (JC .LE. NREADIN2) THEN
            GO TO 210
          END IF
C
C----   Otherwise DEL is FALSE and the entries are the
C       same; flag IC'th element of DELFLAG1 to point to JC
C
        ELSE
          DELFLAG1(IC) = JC
          JSTART = JC + 1
        END IF
C
C---- Do the next entry in the first list
 200  CONTINUE
C
C---- Jump to here to escape the loop if second list is
C     exhausted before the first
 220  CONTINUE
C
C
C---- PART 2: Resort the arrays so that all unmatched entries are
C     discarded.
C
      NATOMS = 0
C
C---- Loop over entries in the first list
      DO 300 IC = 1,NREADIN1
C
        JC = DELFLAG1(IC)
C
        IF (JC .NE. 0) THEN
C
          NATOMS = NATOMS + 1
          X1(1,NATOMS)       = X2(1,JC)
          X1(2,NATOMS)       = X2(2,JC)
          X1(3,NATOMS)       = X2(3,JC)
          ITYPE1(NATOMS)     = ITYPE2(JC)
          ATOMSTORE1(NATOMS) = ATOMSTORE2(JC)
          RESTYP1(NATOMS)    = RESTYP2(JC)
C
          NP1(NATOMS) = NPOINT1(IC)
          NP2(NATOMS) = NPOINT2(JC)
C
        END IF
C
 300  CONTINUE
C
C---- Check for number of actual matches (better late than never)
C
      IF (NATOMS .EQ. 0)
     +   CALL CCPERR(1, 'S/R MATCHUP : Files have no common entries')
C
      DO 400 IC = 1,NATOMS
        NPOINT1(IC) = NP1(IC)
        NPOINT2(IC) = NP2(IC)
 400  CONTINUE
C
C---- Last bits of messaging and tidying up
C
      IF (NREADIN1 .EQ. NATOMS .AND. NREADIN2.EQ.NATOMS) THEN
        WRITE(6,FMT=9040)
 9040   FORMAT(/,'No atoms were removed: area differences may be zero')
      END IF
C
      NREADIN1 = NATOMS
C
      WRITE(6,FMT=9050) NREADIN1
 9050 FORMAT(/,'Finished sorting atoms: total number now = ',I8,//)
C
      RETURN
C
      END
C
C     =========================================================
      SUBROUTINE ATOMSETUP(ATYPE,RADIUS,ATOMN,NATOMTYP,MAXTYPE)
C     =========================================================
C
C     Put data for atom types into temporary arrays, and then read this
C     into the actual arrays. Means I can use data statement without
C     "incomplete initialisers"...
C
      IMPLICIT NONE
C     ..
C     .. Scalars in argument list
      INTEGER MAXTYPE,NATOMTYP
C     ..
C     .. Arrays in argument list
      INTEGER ATOMN(0:MAXTYPE)
      REAL RADIUS(0:MAXTYPE)
      CHARACTER ATYPE(0:MAXTYPE)*2
C     ..
C     .. Local scalars
      INTEGER ICOUNT
C     ..
C     .. Local arrays
      INTEGER TMPN(0:8)
      REAL TMPRAD(0:8)
      CHARACTER TMPTYP(0:8)*2
C
C     ..
C     .. Data statements ..
      DATA TMPTYP/'UN',' C',' N',' O','MG',' S',' P','CL','CO'/
      DATA TMPRAD/1.8,1.8,1.65,1.6,1.6,1.85,1.90,1.80,1.80/
      DATA TMPN / 0, 6, 7, 8, 12, 16, 15, 17, 27 /
C
      NATOMTYP = 8
C
      DO 10 ICOUNT = 0,NATOMTYP
C
        ATOMN(ICOUNT) = TMPN(ICOUNT)
        RADIUS(ICOUNT) = TMPRAD(ICOUNT)
        ATYPE(ICOUNT) = TMPTYP(ICOUNT)
C
 10   CONTINUE
C
      RETURN
C
      END
C
C     ==================================================================
      SUBROUTINE KEYINPUT(ATYPE,RADIUS,ATOMN,NATOMTYP,MAXTYPE,ALL,NOHOH,
     +                    HOH,HOHALL,SOFF,SIMOL,HOHRAD,PNTDEN,IDIFF,
     +                    NSYM,TR,NSYM2,TR2,NREJ,REJTYP,VERB,ITRANS,
     +                    LOUTPUT,LCLOSED,PROBE_STEP,PROBE_MAX,
     +                    USE_COORDS)
C     ==================================================================
C
C     Reads in keyworded input for the following keywords:
C --- 1. MODE ALL | NOHOH | HOH | HOHALL
C --- 2. SMODE IMOL | OFF
C --- 3. SYMMETRY <no> | <name> | <operators>
C --- 4. EXCLUDE <res1,res2, ... ,res30>
C --- 5. PROBE <radius>
C --- 6. ATOM <name> <no> <radius>
C --- 7. PNTDEN <no_of_points>
C --- 8. DIFFMODE OFF | IMOL | WATER | COMPARE | CAVITY
C --- 9. VERBOSE
C --- 10.TRANS NONE | 1 | 2 | BOTH
C --- 11.CAV_PROBE <step> <max>
C --- 12.MATCHUP ALL | NOCOORDS
C --- 13.OUTPUT
C --- There is also an 'END' keyword (automatically supplied by
C     KEYPARSE)
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXRES
      PARAMETER (MAXRES=30)
C     ..
C     .. Scalars in argument list
      INTEGER NATOMTYP,MAXTYPE,IDIFF,ITRANS,NSYM,NSYM2,NREJ
      REAL    HOHRAD,PNTDEN,PROBE_STEP,PROBE_MAX
      LOGICAL ALL,HOH,HOHALL,NOHOH
      LOGICAL SOFF,SIMOL,VERB,LOUTPUT,LCLOSED,USE_COORDS
C     ..
C     .. Arrays in argument list
      REAL RADIUS(0:MAXTYPE),TR(4,4,12),TR2(4,4,12)
      INTEGER ATOMN(0:MAXTYPE)
      CHARACTER ATYPE(0:MAXTYPE)*2,REJTYP(MAXRES)*3
C     ..
C     .. Local scalars ..
      REAL    VDWR
      INTEGER NSYMP,IC,NSPGRP,MODE,SMODE,ATNUM,LENSTR
      LOGICAL SYMM,OLD,FAIL,CONT,EXCEED,LATOM,DIFFMODE,LTRANS
      LOGICAL NO_COORDS
C     ..
C     .. Local arrays ..
      CHARACTER RTYPE*3,ATOMNAME*2,ATMCHECK*20
      CHARACTER NAMSPG*10,PGNAME*10,REST*10
C     ..
C     .. External subroutines/functions ..
      EXTERNAL CCPERR,MEMOPARSE,PARSESYMM,PARSEDIAGNOSE,PARSESUBCHAR
      EXTERNAL PARSEKEY,PARSESUBKEY,PARSEKEYARG,PARSEREAL,CCPUPC,LENSTR
      EXTERNAL PARSESUBINT,PARSESUBREAL
C     ..
C     ..
C
C --- Set up initial values
C
      DIFFMODE = .FALSE.
      IDIFF    = 0
      ALL = .FALSE.
      NOHOH = .FALSE.
      HOH = .FALSE.
      HOHALL = .FALSE.
      SOFF = .FALSE.
      SIMOL = .FALSE.
      SYMM = .FALSE.
      NSYM  = 0
      NSYM2 = 0
      NREJ = 0
      FAIL = .FALSE.
      OLD = .FALSE.
      EXCEED = .FALSE.
      HOHRAD = 1.4
      PNTDEN = 1.0
      ATMCHECK = ' '
      LATOM = .FALSE.
      VERB = .FALSE.
      LTRANS  = .FALSE.
      ITRANS  = 0
      LOUTPUT = .FALSE.
      LCLOSED = .FALSE.
      PROBE_STEP = -1.0
      PROBE_MAX  = 1.4
      USE_COORDS = .TRUE.
      NO_COORDS  = .FALSE.
C
C --- Use KEYPARSE routines to get keywords
C
 10   CALL MemoParse(.TRUE.)
C
C --- First check for "keywords" ALL,HOH,NOHOH,HOHALL
C     If these are detected then it suggests an out of date script
C     (ie non-keyworded input script)
C
      CALL ParseKey('ALL',OLD)
      CALL ParseKey('HOH',OLD)
      CALL ParseKey('HOHA',OLD)
      CALL ParseKey('NOHO',OLD)
C
C --- If so then there is no point doing any more
C
      IF (OLD)
     +  CALL CCPERR (1,' FATAL: out-of-date script, see documentation')
C
C --- Now examine other keywords:
C
C ****************
C --- Keyword DIFFMODE: OFF,IMOL,WATERS (default OFF)
C     Determines type of analysis:
C         OFF    - old AREAIMOL, all other k/w's as before (IDIFF=0)
C         IMOL   - needs two SYMMETRY k/w's  (IDIFF=1)
C         WATERS - doesn't require MODE  (IDIFF=2)
C         COMPARE- currently all k/w's as before (IDIFF=3)
C         CAVITY - similar to OFF I think (IDIFF=4)
C
C --- This must be the first keyword.
C
      IF (.NOT. DIFFMODE) THEN
C
        CALL ParseSubKey('DIFF','OFF',DIFFMODE)
        IF (DIFFMODE) IDIFF = 0
C
        IF (.NOT. DIFFMODE) THEN
          CALL ParseSubKey('DIFF','IMOL',DIFFMODE)
          IF (DIFFMODE) IDIFF = 1
        ENDIF
C
        IF (.NOT. DIFFMODE) THEN
          CALL ParseSubKey('DIFF','WATE',DIFFMODE)
          IF (DIFFMODE) IDIFF = 2
        ENDIF
C
        IF (.NOT. DIFFMODE) THEN
          CALL ParseSubKey('DIFF','COMP',DIFFMODE)
          IF (DIFFMODE) IDIFF = 3
        ENDIF
C
        IF (.NOT. DIFFMODE) THEN
          CALL ParseSubKey('DIFF','CAVI',DIFFMODE)
          IF (DIFFMODE) THEN
            IDIFF = 4
C           DIFFMODE CAVI not implemented for CCP4 4.2
            CALL CCPERR(1,'DIFFMODE CAVITY: option not implemented')
          ENDIF
        ENDIF
C
        DIFFMODE = .TRUE.
C
      ELSE
C
        REST = ' '
        CALL ParseKeyArg('DIFF',REST)
        IF (REST .NE. ' ') THEN
          CALL CCPERR 
     +       (2,' DIFFMODE can only appear once and must be first')
          FAIL = .TRUE.
        ENDIF       
C
      ENDIF
C
C ****************
C --- Keyword MODE: ALL,HOH,NOHOH,HOHALL (default NOHOH)
C     Specifies whether waters etc are included
C
C --- MODE is not required in DIFFMODE WATER (IDIFF=2)
C
      IF (IDIFF.NE.2) THEN
C
        CALL ParseSubKey('MODE','ALL',ALL)
        CALL ParseSubKey('MODE','HOH',HOH)
        CALL ParseSubKey('MODE','HOHA',HOHALL)
        CALL ParseSubKey('MODE','NOHO',NOHOH)
C
      ELSE
C
        REST = ' '
        CALL ParseKeyArg('MODE',REST)
        IF (REST .NE. ' ') THEN
          CALL CCPERR (3,' MODE redundant in DIFFMODE WATER')
        ENDIF
C
      ENDIF
C
C *****************
C --- Keyword SMODE: IMOL,OFF (default OFF)
C     Specifies whether symmetry-related atoms are generated
C
C --- SMODE is not required in DIFFMODE IMOL (IDIFF=1)
C
      IF (IDIFF.NE.1) THEN

        CALL ParseSubKey('SMOD','IMOL',SIMOL)
        CALL ParseSubKey('SMOD','OFF',SOFF)
C
      ELSE
C
        REST = ' '
        CALL ParseKeyArg('SMOD',REST)
        IF (REST .NE. ' ') THEN
          CALL CCPERR (3,' SMODE redundant in DIFFMODE IMOL')
        ENDIF
C
      ENDIF
C
C ********************
C --- Keyword SYMMETRY: (group name),(group no.),(operators)
C     Specifies space group
C
C --- In DIFFMODE IMOL (IDIFF=1) there can be 1 or 2 SYMM keywords
C
      IF (NSYM .EQ. 0) THEN
C
        CALL ParseSYMM(NAMSPG,NSPGRP,PGNAME,NSYM,NSYMP,TR)
C
      ELSE IF ((NSYM2 .EQ. 0) .AND. (IDIFF.EQ.1)) THEN
C
        CALL ParseSYMM(NAMSPG,NSPGRP,PGNAME,NSYM2,NSYMP,TR2)
C
      ELSE
C
        REST = ' '
        CALL ParseKeyArg('SYMM',REST)
        IF (REST .NE. ' ') THEN
          CALL CCPERR (2,' Multiple occurance of SYMMETRY')
          FAIL = .TRUE.
        ENDIF
C
      ENDIF
C
C *******************
C --- Keyword EXCLUDE: (res1 res2 ... res30)
C     Specifies residues to be excluded from the calculation
C
C --- Collect names of excluded residues
C
      IC = 1
 20   CONTINUE
      RTYPE = ' '
      CALL ParseSubChar('EXCL',' ',IC,.FALSE.,RTYPE)
      IF (RTYPE .NE. ' ') THEN
        NREJ = NREJ + 1
C
C ---   convert residue name to uppercase
C
        CALL CCPUPC (RTYPE)
C
C ---   Not allowed more than MAXRES residues ...
C
        IF (NREJ .GT. MAXRES) THEN
C
C ---     Send warning once only
C
          IF (.NOT. EXCEED) THEN
            CALL CCPERR (2,' More than MAXRES residues')
            CALL CCPERR (3,' Increase parameter MAXRES')
            EXCEED = .TRUE.
          ENDIF
          WRITE (6,FMT=7000) RTYPE
          NREJ = MAXRES
        ELSE
          REJTYP(NREJ) = RTYPE
        ENDIF
        IC = IC + 1
        GO TO 20
      ENDIF
C
 7000 FORMAT (/' Discarding residue type ',A3)
C
C *******************
C --- Keyword: PROBE (HOHRAD)
C     Sets the solvent radius of the probe used in the area calculation
C
      CALL ParseReal('PROB',HOHRAD)
C
C *******************
C --- Keyword: ATOM (name) (no.) (vdw radius, /AA)
C     Sets up user-defined atom type with name, atomic number and
C     associated van der Waal's radius
C
      ATOMNAME = ' '
      ATNUM = 0
      VDWR = 0
C
      CALL ParseSubChar('ATOM',' ',1,.TRUE.,ATOMNAME)
      CALL ParseSubInt('ATOM',' ',2,.TRUE.,ATNUM)
      CALL ParseSubReal('ATOM',' ',3,.TRUE.,VDWR)
C
      IF (ATOMNAME .NE. ' ') THEN
C
C --- Check if input is VAGUELY sensible ...
C
        IF (ATNUM.LE.0 .OR. VDWR.LE.0) THEN
          LATOM = .TRUE.
          CALL CCPERR (2,'ATOM no. or radius <= 0 : definition skipped')
          GO TO 50
        ENDIF
C
C --- Convert atom name to uppercase
C
        CALL CCPUPC(ATOMNAME)
C
C --- Right-justify atom name if necessary
C
        IF (LENSTR(ATOMNAME) .EQ. 1 .AND. ATOMNAME(1:1) .NE. ' ') THEN
          ATOMNAME(2:2) = ATOMNAME(1:1)
          ATOMNAME(1:1) = ' '
        ENDIF
C
C --- Check: is this atom type already defined?
C     (must match name AND no.) - if so then overwrite old radius
C
         DO 40 IC = NATOMTYP,1,-1
C
           IF (ATOMN(IC) .EQ. ATNUM) THEN
             IF (ATYPE(IC) .EQ. ATOMNAME) THEN
                RADIUS(IC) = VDWR
                ATMCHECK(IC:IC) = '*'
             ELSE
                LATOM = .TRUE.
                CALL CCPERR (2,
     $          'ATOM name/no. not consistent : definition skipped')
             ENDIF
             GO TO 50
           ENDIF
C
           IF (ATYPE(IC).EQ.ATOMNAME .AND. ATOMN(IC).NE.ATNUM) THEN
             LATOM = .TRUE.
             CALL CCPERR (2,
     $       'ATOM name/no. not consistent : definition skipped')
             GO TO 50
           ENDIF
C
 40      CONTINUE
C
C --- Add to the list of atoms
C
         CALL ADDATOMTYPE(MAXTYPE,ATYPE,RADIUS,ATOMN,ATMCHECK,
     +        ATOMNAME,VDWR,ATNUM,NATOMTYP)
C
      ENDIF
C
 50   CONTINUE
C
C *******************
C --- Keyword: PNTDEN (PNTDEN)
C     Sets the number of points/angstrom^2 used in the area calculation
C
      CALL ParseReal('PNTD',PNTDEN)
C
C *******************
C --- Keyword: VERBOSE
C     Turns on additional output
C
      CALL ParseKey('VERB',VERB)
C
C *******************
C --- Keyword: TRANS: NONE (default), 1, 2, BOTH
C     Turns on (or off) generation of translation vectors during
C     symmetry generation.
C
C     For DIFFMODE IMOL the options are:
C         NONE:   no translation vectors (ITRANS=0)
C         1   :   translations in 1st run only (ITRANS=1)
C         2   :   translations in 2nd run only (ITRANS=2)
C         BOTH:   translations in both 1st and 2nd runs (ITRANS=3)
C
C     For SMODE IMOL the TRANS keyword alone is required to switch on
C     the vectors (ITRANS=1); TRANS NONE (default) switches them off.
C
      LTRANS = .FALSE.
C
      CALL ParseSubKey('TRAN','NONE',LTRANS)
      IF (LTRANS) ITRANS = 0
C
      IF (IDIFF.EQ.1) THEN
        IF (.NOT. LTRANS) THEN
          CALL ParseSubKey('TRAN','1',LTRANS)
          IF (LTRANS) ITRANS = 1
        ENDIF
        IF (.NOT. LTRANS) THEN
          CALL ParseSubKey('TRAN','2',LTRANS)
          IF (LTRANS) ITRANS = 2
        ENDIF
        IF (.NOT. LTRANS) THEN
          CALL ParseSubKey('TRAN','BOTH',LTRANS)
          IF (LTRANS) ITRANS = 3
        ENDIF
      ELSE
        IF (.NOT. LTRANS) THEN
          CALL ParseKey('TRAN',LTRANS)
          IF (LTRANS) ITRANS = 1
        ENDIF
      END IF
C
C *****************
C --- Keyword CAV_PROBE <step> <max>
C     Specifies the probe sphere radius increment and maximum cutoff
C     when doing the search for cavities
C
C --- CAV_PROBE can only be used in DIFFMODE CAVITY (IDIFF=4)
C
      IF (IDIFF.EQ.4) THEN

        CALL ParseSubReal('CAV_',' ',1,.TRUE.,PROBE_STEP)
        CALL ParseSubReal('CAV_',' ',2,.TRUE.,PROBE_MAX)
C
      ENDIF
C
C *****************
C --- Keyword MATCHUP <mode>
C     Specifies the criteria for the matchup process in DIFFMODE COMPARE
C
C     Allowed modes are ALL (use all criteria, including coordinate data -
C     this is the default)
C     and NOCOORDS (don't use coordinate data)
C
C --- MATCHUP can only be used in DIFFMODE COMPARE (IDIFF=3)
C
      IF (IDIFF.EQ.3) THEN
C
        CALL ParseSubKey('MATC','ALL',USE_COORDS)
        CALL ParseSubKey('MATC','NOCO',NO_COORDS)
C
      END IF
C
C *******************
C --- Keyword: OUTPUT
C     Turns on output of psuedo-pdb output file containing
C     accessible areas.
C
      CALL ParseKey('OUTP',LOUTPUT)
C
C *******************
C --- CHECK: Is there any more input?
C
      CALL ParseDiagnose(CONT)
      IF (CONT) GO TO 10
C
C --- PERFORM CHECKING ... 
C
C --- Test to see if MODE is set
C
      IF (IDIFF .NE. 2) THEN
        MODE = 0
        IF (ALL) MODE = MODE + 1
        IF (HOH) MODE = MODE + 1
        IF (NOHOH) MODE = MODE + 1
        IF (HOHALL) MODE = MODE +1
C
        IF (MODE .EQ. 0) THEN
          CALL CCPERR (3, ' MODE not set: defaulting to NOHOH')
          NOHOH = .TRUE.
        ELSE IF (MODE .GE. 2) THEN
          CALL CCPERR (2, ' Multiple occurance of MODE')
          FAIL = .TRUE.
        ENDIF
      ENDIF
C
C --- Test to see if SMODE is set
C
      IF (IDIFF .NE. 1) THEN
        SMODE = 0
        IF (SIMOL) SMODE = SMODE + 1
        IF (SOFF) SMODE = SMODE + 1
C
        IF (SMODE .EQ. 0) THEN
          CALL CCPERR (3, ' SMODE not set: defaulting to OFF')
          SOFF = .TRUE.
        ELSE IF (SMODE .GE. 2) THEN
          CALL CCPERR (2, ' Multiple occurance of SMODE')
          FAIL = .TRUE.
        ENDIF
      ENDIF
C
C --- Test to see if symmetry operators have been supplied for IMOL
C     If not then default to SMODE = OFF
C
      IF (SIMOL .AND. (NSYM.EQ.0)) THEN
        CALL CCPERR (3,' Missing SYMMETRY operations: using P1')
      ENDIF
C
C --- If FAIL then something is wrong with input
C
      IF (FAIL) CALL CCPERR (1, 'ERROR(S) IN KEYWORDED INPUT')
C
C --- WRITE OUT MESSAGES TO SCREEN FOR DIFFERENT FUNCTIONS
C
      call ccp4h_summary_beg()
      call ccp4h_header('Mode','mode',3)
C --- First report which DIFFMODE we are running in
      IF (IDIFF.EQ.0) THEN
        WRITE (6,FMT=8010)
      ELSE IF (IDIFF.EQ.1) THEN
        SOFF  = .FALSE.
        SIMOL = .TRUE.
      ELSE IF (IDIFF.EQ.3) THEN
        WRITE (6,FMT=8013)
      ELSE IF (IDIFF.EQ.4) THEN
        WRITE (6,FMT=8014)
      ENDIF
C
 8010 FORMAT
     + (/,' DIFFMODE OFF: single area calculation')
 8013 FORMAT
     + (/,' DIFFMODE COMPARE: area differences between two sets of',
     +  /,'                   coordinates',/
     +  /,' Accessible areas will be calculated for both sets of atoms',
     +  /,' individually, and area differences will be calculated for',
     +  /,' those atoms which are common to both sets.')
 8014 FORMAT
     + (/,' DIFFMODE CAVITY: search for possible cavities')
C
C --- Values of MODE
      IF (IDIFF .EQ. 2) THEN
        WRITE (6,FMT=6010)
      ELSE IF (ALL) THEN
        WRITE (6,FMT=6002)
      ELSE IF (NOHOH) THEN
        WRITE (6,FMT=6004)
      ELSE IF (HOH) THEN
        WRITE (6,FMT=6006)
      ELSE IF (HOHALL) THEN
        WRITE (6,FMT=6008)
      ENDIF
C
 6002 FORMAT
     + (/' MODE ALL - accessible area will be calculated for all',
     +  /'            residues (including waters)',
     +  /'            Waters will be treated as solvent',
     + //' **WARNING: Use this mode with care - the presence of waters',
     +  /'            may lead to unrealistically inflated estimates',
     +  /'            of the accessible area (see documentation)')
 6004 FORMAT
     + (/' MODE NOHOH - Residue types HOH and WAT will be ignored')
 6006 FORMAT
     + (/' MODE HOH - Only residue types HOH and WAT will be considered'
     + ,/'            Waters will be treated as solvent')
 6008 FORMAT
     + (/' MODE HOHALL - Only residue types HOH and WAT will be',
     +  /'               considered',
     +  /'               Waters will be treated as protein atoms')
 6010 FORMAT
     + (/' DIFFMODE WATER - difference between MODE HOH and MODE',
     +  /'                  HOHALL',/,
     +  /' Only residue types HOH and WAT will be considered',/,
     +  /' The accessible area will be calculated twice, first',
     +  /' treating waters as solvent and then as protein')
C
C --- Values of SMODE
      IF (IDIFF .EQ. 1) THEN
        WRITE (6,FMT=6052)
      ELSE IF (SOFF) THEN
        WRITE (6,FMT=6050)
      ELSE IF (SIMOL) THEN
        WRITE (6,FMT=6051)
      ENDIF
      call ccp4h_summary_end()
C
 6050 FORMAT
     + (/' SMODE OFF - symmetry-related atoms will not be generated')
 6051 FORMAT
     + (/' SMODE IMOL - symmetry-related atoms will be generated from',
     +  /'              symmetry operations')
 6052 FORMAT
     + (/' DIFFMODE IMOL - differences due to intermolecular contacts',
     +  /' Separate accessible area calculations will be performed for',
     +  /' each set of symmetry operators.'
     +  /' The differences in accessible area will then be calculated')
C
C----  Apply translations option
C
      IF (IDIFF.EQ.1) THEN
C
        IF (ITRANS .EQ. 0) THEN
          WRITE (6,FMT=6060)
        ELSE IF (ITRANS .EQ. 1) THEN
          WRITE (6,FMT=6061)
        ELSE IF (ITRANS .EQ. 2) THEN
          WRITE (6,FMT=6062)
        ELSE IF (ITRANS .EQ. 3) THEN
          WRITE (6,FMT=6063)
        END IF
C
      ELSE IF (SIMOL) THEN
C
        IF (ITRANS .EQ. 0) THEN
          WRITE (6,FMT=6060)
        ELSE IF (ITRANS .EQ. 1) THEN
          WRITE (6,FMT=6064)
        END IF
C
      END IF
C
 6060 FORMAT
     + (/' TRANSlation mode NONE: lattice vector translations will not'
     +  /' be applied during the generation of symmetry-related'
     +  /' molecules')
 6061 FORMAT
     + (/' TRANSlation mode 1: lattice vector translations will be'
     +  /' applied during the generation of symmetry-related molecules'
     +  /' for the first area calculation only')
 6062 FORMAT
     + (/' TRANSlation mode 2: lattice vector translations will be'
     +  /' applied during the generation of symmetry-related molecules'
     +  /' for the second area calculation only')
 6063 FORMAT
     + (/' TRANSlation mode BOTH: lattice vector translations will be'
     +  /' applied during the generation of symmetry-related molecules'
     +  /' for both the first and second area calculations')
 6064 FORMAT
     + (/' TRANSlation mode ON: lattice vector translations will be'
     +  /' applied during the generation of symmetry-related molecules')
C
C----  Matchup criteria in DIFFMODE COMPARE
C
      IF (IDIFF.EQ.3) THEN
        IF (NO_COORDS) USE_COORDS = .FALSE.
        IF (USE_COORDS) THEN
          WRITE (6,FMT='(/,2A,/,2A)')
     +         ' All data will be used to perform matchup',
     +         ' of XYZIN and XYZIN2'
        ELSE 
          WRITE (6,FMT='(/,2A,/,2A)')
     +         ' Atomic coordinate data will not be used to perform',
     +         ' matchup of XYZIN and XYZIN2'
        END IF
      END IF
C
C----  Solvent radius
C
      IF (HOHRAD.GT.25.0) THEN
        CALL CCPERR (2,' PROBE radius > 25.0 A')
        HOHRAD = 25.0
      ELSEIF (HOHRAD.LE.0.0) THEN
        CALL CCPERR (2,' PROBE radius zero or less')
        HOHRAD = 1.4
      ENDIF
C
      WRITE (6,FMT=6025) HOHRAD
 6025 FORMAT (/,' Solvent/probe radius set to ',F6.2,' A')
C
C----  Point density
C
      IF (PNTDEN.LE.0.) THEN
        CALL CCPERR (2,' PNTDEN less than or equal to 0: set to 1.0')
        PNTDEN = 1.0
      ENDIF
C
      WRITE (6,FMT=6030) PNTDEN
C
 6030 FORMAT
     + (/,' Number of points/area used in AREA calculation set to '
     +  /,F6.2,' per angstrom squared',/)
C
C----  Output file
C
      IF (LOUTPUT) WRITE (6,FMT='(/,A,/,A,/)')
     +             ' Areas will be written to PDB-style output file',
     +             ' (area/area difference written to B-factor column)'
C
C----  Search for closed surfaces (only switched on for MODE NOHOH or ALL)
C
      IF (NOHOH .OR. ALL) THEN
        LCLOSED = .TRUE.
        WRITE (6,FMT='(/,2A,/,2A)')
     +             ' Search for isolated areas of surface will be',
     +             ' performed',' This will identify unconnected',
     +             ' surface areas, including cavities'
        IF (.NOT. SOFF .OR. IDIFF.EQ.1) THEN
          WRITE (6,FMT='(2A,/)')' inside the molecule or formed as a',
     +    ' result of intermolecular contacts'
        ELSE
          WRITE (6,FMT='(A,/)')' inside the molecule'
        END IF
      END IF
C
C---- List of defined atoms and radii
C
      IF (VERB) THEN
        WRITE (6,FMT=7050) NATOMTYP
        DO 30 IC = 1,NATOMTYP
          WRITE (6,FMT=7055) ATMCHECK(IC:IC), IC, ATYPE(IC), ATOMN(IC),
     +                       RADIUS(IC)
 30     CONTINUE
        IF (LATOM) WRITE (6,FMT=7060)
      END IF
C
 7050 FORMAT (/' Areaimol recognises the following ',I4,' atom types:',
     +        /' (* = user-defined atom/radius)',
     +        //'       No.   Name   Atmic no.  VdW rad. (A)',
     +         /'    ---------------------------------------')
 7055 FORMAT (5X,A1,I3,3X,A4,5X,I4,4X,F8.2)
 7060 FORMAT (/'NB: One or more ATOM keywords had inconsistent',
     +        /'arguments and so have been ignored',/)
C
C---- Write out residue types to be ignored (if any)
C
      IF (NREJ.NE.0) WRITE (6,FMT=6024) (REJTYP(IC),IC=1,NREJ)
 6024 FORMAT (//,
     +' The following residue types will be ignored in the area',/,
     +' calculations and will not be written to the output file:',//,
     + (1X,15 (A,2X)))
C
C
      RETURN
C
      END
C
C     =================================================================
      SUBROUTINE DIFFAREA(NREADIN,NPOINT,NPOINT2,PNTDEN,ATYPE,NATOMTYP,
     +                    ITYPE,RESTYP,IDIFF)
C     =================================================================
C     Compare solvent accessible areas after two areaimol runs
C     and output an array with differences in surface area points for
C     each atom.
C
C     Originally this subroutine was the unsupported program DIFFAREA.
C
C  INPUT
C
C    NREADIN (integer)       number of atoms originally read in
C
C    NPOINT (integer array)  no of surface points for each atom from
C                            first area calc.
C
C    NPOINT2 (integer array) no of surface points for each atom from
C                            second area calc.
C
C    PNTDEN (real)           point density (recipr. of area per point)
C
C    ATYPE (character*2 array)    list of recognised atomic names
C
C    NATOMTYP (integer)      number of recognised atom types
C
C    ITYPE (integer array)   identifies the type of each stored atom
C
C    RESTYP (logical*1 array) TRUE if atom is protein (not water)
C
C    IDIFF (integer)         flag identifying which DIFFMODE is in use
C
C  OUTPUT
C
C    NPOINT (integer array)  on output holds the difference of NPOINT and
C                            NPOINT2
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXTYPE
      PARAMETER (MAXTYPE=20)
C
C     .. Scalar arguments
      INTEGER NREADIN,NATOMTYP,IDIFF
      REAL PNTDEN
C     .. Array arguments
      INTEGER ITYPE(NREADIN),NPOINT(NREADIN),NPOINT2(NREADIN)
      CHARACTER ATYPE(0:MAXTYPE)*2
      LOGICAL*1 RESTYP(NREADIN)
      LOGICAL PROT
C     ..
C     .. Local Scalars ..
      REAL ACCESS1,ACCESS2,SUM1,SUM2,TDIFF
      INTEGER I,NOUT,IATYP,NATOMS
C     ..
C     .. Local Arrays ..
      REAL DIFF(MAXTYPE),SUMACCESS1(MAXTYPE),SUMACCESS2(MAXTYPE)
      INTEGER NDIFF(MAXTYPE),NTYPE(MAXTYPE)
C     ..
C
C --- WRITE INITIAL INFORMATION
C
      WRITE (6,FMT=6010)
 6010 FORMAT (//,11X,'================================================',
     +        //,14X,'*****  ANALYSIS OF AREA DIFFERENCES  *****',
     +         /,14X,'       ----------------------------')
C
      IF (IDIFF.EQ.1) THEN
        WRITE(6,FMT='(//,A,/,A)')
     +  ' DIFFMODE IMOL: area differences are due to the presence of',
     +  '                symmetry related atoms'
      ELSE IF (IDIFF.EQ.2) THEN
        WRITE(6,FMT='(//,A,/,A)')
     +  ' DIFFMODE WATERS: area differences are due to treating waters',
     +  '                  as either solvent or protein'
      ELSE IF (IDIFF.EQ.3) THEN
        WRITE(6,FMT='(//,A,/,A,/,A)')
     +  ' DIFFMODE COMPARE: area differences are due to the presence',
     +  '                   or absence of atoms between the two input',
     +  '                   files'
      END IF
C
C --- INITIALISE SUMS AND DIFFERENCES
C
      NOUT   = 0
      NATOMS = 0
C
      DO 1 I= 1,NATOMTYP
        DIFF(I) = 0.0
        SUMACCESS1(I) = 0.0
        SUMACCESS2(I) = 0.0
        NDIFF(I) = 0
        NTYPE(I) = 0
 1    CONTINUE
C
      DO 100 I = 1,NREADIN
C
C --- IN DIFFMODE WATER (IDIFF=2) IGNORE PROTEIN
C
        PROT = RESTYP(I)
        IF ((IDIFF .EQ. 2) .AND. (RESTYP(I) .EQV. .TRUE.)) GO TO 190
        NATOMS = NATOMS + 1
C
C----  Accessible area
C
        ACCESS1 = NPOINT(I)/PNTDEN
        ACCESS2 = NPOINT2(I)/PNTDEN
C
C----  Find atom type
C
      IATYP = ITYPE(I)
C
      SUMACCESS1(IATYP) = SUMACCESS1(IATYP) + ACCESS1
      SUMACCESS2(IATYP) = SUMACCESS2(IATYP) + ACCESS2
      NTYPE(IATYP) = NTYPE(IATYP) + 1
      NPOINT(I) = NPOINT(I) - NPOINT2(I)
C
      IF (NPOINT(I).NE.0) THEN
        NOUT = NOUT + 1
        NDIFF(IATYP) = NDIFF(IATYP) + 1
      END IF
C
 190  CONTINUE
C
 100  CONTINUE
C
      SUM1 = 0.0
      SUM2 = 0.0
      TDIFF = 0.0
      DO 60 I = 1,NATOMTYP
        SUM1 = SUMACCESS1(I) + SUM1
        SUM2 = SUMACCESS2(I) + SUM2
        DIFF(I) = SUMACCESS1(I) - SUMACCESS2(I)
        TDIFF = DIFF(I) + TDIFF
   60 CONTINUE
C
C----  Write out some stats
C
      call ccp4h_summary_beg()
      WRITE (6,FMT=6004)
 6004 FORMAT (//' Atom-type  Number   Area1    Area2  N-diff Area-diff',
     +       /)
      DO 1010 I = 1,NATOMTYP
        IF (NTYPE(I) .NE. 0) THEN
          WRITE (6,FMT=6006) ATYPE(I),NTYPE(I),SUMACCESS1(I),
     +    SUMACCESS2(I),NDIFF(I),DIFF(I)
        END IF
 1010 CONTINUE
 6006 FORMAT (4X,A,5X,I6,2F9.1,1X,I6,F9.1)
      WRITE (6,FMT=6008) NATOMS,SUM1,SUM2,NOUT,TDIFF
 6008 FORMAT (/' Totals    ',I6,2F9.1,1X,I6,F9.1,/)
      call ccp4h_summary_end()
C
      RETURN
C
      END
C
C     ==============================================================
      SUBROUTINE RESAREA(NATOM,ATOMSTORE,NPOINT,PNTDEN,RESTYP,NOHOH,
     +                   FLAG)
C     ==============================================================
C
C  Print out solvent accessible areas for each residue ,chain and
C  for the whole molecule.
C
C  Adapted from the program RESAREA.
C
C  ARGUMENTS:
C
C  NATOM (integer) - no. of atoms in list
C  ATOMSTORE (character*10 array) - part of ATOM line for each atom
C  NPOINT (integer array) - no. of surface points for each atom
C  PNTDEN (real) - no. of surface points per square Angstrom
C  RESTYP (logical*1 array) - flag set to TRUE if atom is protein
C  NOHOH (logical) - flag set to TRUE if MODE is NOHOH
C  FLAG (integer) - if 1 then raw areas, if zero then area differences
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER ILINE
      PARAMETER (ILINE=3)
C     ..
C     .. Local Scalars ..
      REAL ACCESS,CHAIN,SUMACCESS,TOTAL,PNTDEN
      INTEGER I,IPR,NATOM,ICOUNT,JCOUNT,FLAG
      CHARACTER CHNAME*1,RESNAM*3,SEQID*6,STRING*10
      LOGICAL NOHOH,WATER
C     ..
C     .. Local Arrays ..
      REAL SUMA(ILINE)
      INTEGER NPOINT(MAXATOM)
      CHARACTER RESNA(ILINE)*3,SEQA(ILINE)*6,ATOMSTORE(MAXATOM)*10
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C
      IPR = 0
C
C----  Introduce yourself
C
      IF (FLAG.EQ.1) THEN
        WRITE(6,FMT=7010)
 7010   FORMAT(//,' ANALYSIS OF ACCESSIBLE AREAS BY RESIDUE',/,
     +            ' ---------------------------------------',//)
      ELSE
        WRITE(6,FMT='(/,A,/,A,//)')' Differences for individual chains',
     +                             ' ---------------------------------'
      ENDIF    
C
C----  Read first residue
C
      ICOUNT = 1
 10   STRING = ATOMSTORE(ICOUNT)
      ACCESS = NPOINT(ICOUNT)/PNTDEN
C
      WATER = (.NOT. RESTYP(ICOUNT))
C
C---- check for water in MODE NOHOH
C
      IF (WATER .AND. NOHOH) THEN
C
C---- check whether the area difference is non-zero
C
        IF ((FLAG.EQ.0 .AND. NPOINT(ICOUNT).EQ.0) .OR. FLAG.EQ.1) THEN
C
          ICOUNT = ICOUNT + 1
          IF (ICOUNT.LE.NATOM) GO TO 10
          CALL CCPERR(1,'RESAREA: NO PROTEIN IN OUTPUT')
          RETURN
C
        ENDIF
C
      ENDIF
C
      RESNAM    = STRING(1:3)
      CHNAME    = STRING(5:5)
      SEQID     = STRING(5:10)
      SUMACCESS = ACCESS
      CHAIN     = ACCESS
      TOTAL     = ACCESS
C
C----  Main loop
C
      DO 30 JCOUNT = (ICOUNT+1),NATOM
C
C----  Check: ignore if this is water and MODE NOHOH
C
        IF (.NOT.(WATER .AND. NOHOH)) THEN
C
C---- then check if area difference is greater than zero
C
          IF (FLAG.EQ.1 .OR. (FLAG.EQ.0 .AND. NPOINT(JCOUNT).NE.0)) THEN
C
            STRING = ATOMSTORE(JCOUNT)
            ACCESS = NPOINT(JCOUNT)/PNTDEN
C
            TOTAL = TOTAL + ACCESS
C
            IF (STRING(5:10).EQ.SEQID) THEN
              SUMACCESS = SUMACCESS + ACCESS
            ELSE
C
C----  New residue
C
C  save info to print several residues per line
C
              IPR = IPR + 1
              RESNA(IPR) = RESNAM
              SEQA(IPR) = SEQID
              SUMA(IPR) = SUMACCESS
C
C
              IF (IPR.EQ.ILINE) THEN
                WRITE (6,FMT=6000) (RESNA(I),SEQA(I),SUMA(I),I=1,ILINE)
 6000           FORMAT (6 (5X,A,1X,A,2X,F8.2))
                IPR = 0
              END IF
C
              RESNAM = STRING(1:3)
              SEQID = STRING(5:10)
              SUMACCESS = ACCESS
            END IF
C
C
            IF (STRING(5:5).EQ.CHNAME) THEN
              CHAIN = CHAIN + ACCESS
            ELSE
C
C---- new chain
C
C  write save residues
C
              WRITE (6,FMT=6000) (RESNA(I),SEQA(I),SUMA(I),I=1,IPR)
              IPR = 0
C
              IF (FLAG.EQ.1) THEN
                call ccp4h_summary_beg()
                WRITE (6,FMT=6002) CHNAME,CHAIN
 6002           FORMAT (/' Total area of chain ',A,': ',F10.1,/)
                call ccp4h_summary_end()
              ELSE
                call ccp4h_summary_beg()
                WRITE (6,FMT=6003) CHNAME,CHAIN
 6003           FORMAT (/' Total area difference of chain ',A,
     +                   ': ',F10.1,/)
                call ccp4h_summary_end()
              ENDIF
C
              WRITE (6,FMT=7020)
 7020         FORMAT (' -------------------------------------------',//)
              CHNAME = STRING(5:5)
              CHAIN = ACCESS
            END IF
C
          ENDIF
C
        END IF
C
 30   CONTINUE
C
C----  End of file  - write out final residues
C
   40 WRITE (6,FMT=6888) (RESNA(I),SEQA(I),SUMA(I),I=1,IPR),RESNAM,
     +  SEQID,SUMACCESS
 6888     FORMAT (6 (5X,A,1X,A,2X,F8.2))
      IF (FLAG.EQ.1) THEN
        call ccp4h_summary_beg()
        WRITE (6,FMT=6002) CHNAME,CHAIN
        WRITE (6,FMT=7020)
        WRITE (6,FMT=6004) TOTAL
 6004   FORMAT (' TOTAL AREA: ',F10.1,/)
        call ccp4h_summary_end()
      ELSE
        call ccp4h_summary_beg()
        WRITE (6,FMT=6003) CHNAME,CHAIN
        WRITE (6,FMT=7020)
        WRITE (6,FMT=6005) TOTAL
 6005   FORMAT (' TOTAL AREA DIFFERENCE: ',F10.1,/)
        call ccp4h_summary_end()
      ENDIF
C
      RETURN
C
      END
C
C     ==============================================================
      SUBROUTINE CONTACTAREA(NATOM,ATOMSTORE,NPOINT,PNTDEN,RESTYP,
     +     NOHOH,HOHRAD,ITYPE,NATOMTYP,RADIUS,FLAG)
C     ==============================================================
C
C  Print out solvent contact areas for each residue, chain and
C  for the whole molecule.
C
C  The contact area is the area on the VdW surface of an atom which
C  can be contacted by the probe sphere
C  This turns out to be equal to the accessible area scaled by a
C  factor of (R^2)/(R+R_probe)^2 where R and R_probe are the radii
C  of the atom in question and the probe sphere respectively.
C
C  Adapted from the subroutine RESAREA.
C
C  ARGUMENTS:
C
C  NATOM (integer) - no. of atoms in list
C  ATOMSTORE (character*10 array) - part of ATOM line for each atom
C  NPOINT (integer array) - no. of surface points for each atom
C  PNTDEN (real) - no. of surface points per square Angstrom
C  RESTYP (logical*1 array) - flag set to TRUE if atom is protein
C  NOHOH (logical) - flag set to TRUE if MODE is NOHOH
C  HOHRAD (real) - probe radius
C  ITYPE (integer) - the atom type for each atom
C  NATOMTYP (integer) - number of atom types
C  RADIUS (real) - expanded VdW radius for each atom type
C  FLAG (integer) - if 1 then raw areas, if zero then area differences
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER ILINE
      PARAMETER (ILINE=3)
C     ..
C     .. Scalars ..
      REAL CONTACT,CHAIN,SUMCONTACT,TOTAL,PNTDEN,SCALE,HOHRAD
      INTEGER I,IPR,NATOM,ICOUNT,JCOUNT,ATYPE,NATOMTYP,FLAG
      CHARACTER CHNAME*1,RESNAM*3,SEQID*6,STRING*10
      LOGICAL NOHOH,WATER
C     ..
C     .. Arrays ..
      REAL SUMA(ILINE),RADIUS(0:NATOMTYP)
      INTEGER NPOINT(MAXATOM),ITYPE(MAXATOM)
      CHARACTER RESNA(ILINE)*3,SEQA(ILINE)*6,ATOMSTORE(MAXATOM)*10
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C
      IPR = 0
C
C----  Introduce yourself
C
      IF (FLAG.EQ.1) THEN
        WRITE(6,FMT=7010)
 7010   FORMAT(//,' ANALYSIS OF CONTACT AREAS BY RESIDUE',/,
     +            ' ------------------------------------',//)
      ELSE
        WRITE(6,FMT='(/,A,/,A,//)')
     +       ' Contact area differences for individual chains',
     +       ' ----------------------------------------------'
      ENDIF    
C
C----  Read first residue
C
      ICOUNT = 1
 10   STRING = ATOMSTORE(ICOUNT)
C
C---- Scale for accessible area to contact area for this atom
C
      ATYPE = ITYPE(ICOUNT)
      SCALE = ((RADIUS(ATYPE)-HOHRAD)/RADIUS(ATYPE))**2.0
      CONTACT = NPOINT(ICOUNT)/PNTDEN*SCALE
C
      WATER = (.NOT. RESTYP(ICOUNT))
C
C---- check for water in MODE NOHOH
C
      IF (WATER .AND. NOHOH) THEN
C
C---- check whether the area difference is non-zero
C
        IF ((FLAG.EQ.0 .AND. NPOINT(ICOUNT).EQ.0) .OR. FLAG.EQ.1) THEN
C
          ICOUNT = ICOUNT + 1
          IF (ICOUNT.LE.NATOM) GO TO 10
          CALL CCPERR(1,'CONTACTAREA: NO PROTEIN IN OUTPUT')
          RETURN
C
        ENDIF
C
      ENDIF
C
      RESNAM    = STRING(1:3)
      CHNAME    = STRING(5:5)
      SEQID     = STRING(5:10)
      SUMCONTACT = CONTACT
      CHAIN     = CONTACT
      TOTAL     = CONTACT
C
C----  Main loop
C
      DO 30 JCOUNT = (ICOUNT+1),NATOM
C
C----  Check: ignore if this is water and MODE NOHOH
C
        IF (.NOT.(WATER .AND. NOHOH)) THEN
C
C---- then check if area difference is greater than zero
C
          IF (FLAG.EQ.1 .OR. (FLAG.EQ.0 .AND. NPOINT(JCOUNT).NE.0)) THEN
C
            STRING = ATOMSTORE(JCOUNT)
C
C---- Scale for accessible area to contact area for this atom
            ATYPE = ITYPE(JCOUNT)
            SCALE = ((RADIUS(ATYPE)-HOHRAD)/RADIUS(ATYPE))**2.0
            CONTACT = NPOINT(JCOUNT)/PNTDEN*SCALE
C
            TOTAL = TOTAL + CONTACT
C
            IF (STRING(5:10).EQ.SEQID) THEN
              SUMCONTACT = SUMCONTACT + CONTACT
            ELSE
C
C----  New residue
C
C  save info to print several residues per line
C
              IPR = IPR + 1
              RESNA(IPR) = RESNAM
              SEQA(IPR) = SEQID
              SUMA(IPR) = SUMCONTACT
C
C
              IF (IPR.EQ.ILINE) THEN
                WRITE (6,FMT=6000) (RESNA(I),SEQA(I),SUMA(I),I=1,ILINE)
 6000           FORMAT (6 (5X,A,1X,A,2X,F8.2))
                IPR = 0
              END IF
C
              RESNAM = STRING(1:3)
              SEQID = STRING(5:10)
              SUMCONTACT = CONTACT
            END IF
C
C
            IF (STRING(5:5).EQ.CHNAME) THEN
              CHAIN = CHAIN + CONTACT
            ELSE
C
C---- new chain
C
C  write save residues
C
              WRITE (6,FMT=6000) (RESNA(I),SEQA(I),SUMA(I),I=1,IPR)
              IPR = 0
C
              IF (FLAG.EQ.1) THEN
                call ccp4h_summary_beg()
                WRITE (6,FMT=6002) CHNAME,CHAIN
 6002           FORMAT (/' Total contact area of chain ',A,': ',F10.1,/)
                call ccp4h_summary_end()
              ELSE
                call ccp4h_summary_beg()
                WRITE (6,FMT=6003) CHNAME,CHAIN
 6003           FORMAT (/' Total contact area difference of chain ',A,
     +                   ': ',F10.1,/)
                call ccp4h_summary_end()
              ENDIF
C
              WRITE (6,FMT=7020)
 7020         FORMAT (' -------------------------------------------',//)
              CHNAME = STRING(5:5)
              CHAIN = CONTACT
            END IF
C
          ENDIF
C
        END IF
C
 30   CONTINUE
C
C----  End of file  - write out final residues
C
   40 WRITE (6,FMT=6888) (RESNA(I),SEQA(I),SUMA(I),I=1,IPR),RESNAM,
     +  SEQID,SUMCONTACT
 6888     FORMAT (6 (5X,A,1X,A,2X,F8.2))
      IF (FLAG.EQ.1) THEN
        call ccp4h_summary_beg()
        WRITE (6,FMT=6002) CHNAME,CHAIN
        WRITE (6,FMT=7020)
        WRITE (6,FMT=6004) TOTAL
 6004   FORMAT (' TOTAL CONTACT AREA: ',F10.1,/)
        call ccp4h_summary_end()
      ELSE
        call ccp4h_summary_beg()
        WRITE (6,FMT=6003) CHNAME,CHAIN
        WRITE (6,FMT=7020)
        WRITE (6,FMT=6005) TOTAL
 6005   FORMAT (' TOTAL CONTACT AREA DIFFERENCE: ',F10.1,/)
        call ccp4h_summary_end()
      ENDIF
C
      RETURN
C
      END
C
C     ==============================================================
      SUBROUTINE WATERAREA(NATOM,ATOMSTORE,NPOINT,PNTDEN,HOH,HOHALL,
     +                     FLAG)
C     ==============================================================
C
C  9 AUG 88   Changed to print 3 residues per line  A.G.W.L.
C
C  Analyse solvent accessible areas for water molecules, dividing
C  them into clases with 0, lt 5, lt 10, gt 10 Ang**2 area
C
C  This program only considers residues type HOH or WAT, otherwise is
C  similar to RESAREA.
C  A.G.W.L. July 1988
C
C  Adapted from the program WATERAREA.
C
C  ARGUMENTS:
C
C  NATOM (integer) - no. of atoms in list
C  ATOMSTORE (character*10 array) - part of ATOM line for each atom
C  NPOINT (integer array) - no. of surface points for each atom
C  PNTDEN (real) - no. of surface points per square Angstrom
C  HOH (logical) - if true then we are in MODE HOH
C  FLAG (integer) - if 1 then raw areas, if zero then area differences
C
      IMPLICIT NONE
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
      INTEGER ILINE
      PARAMETER (ILINE=3)
C     ..
C     .. Local Scalars ..
      REAL ACCESS,CHAIN,TOTAL,PNTDEN
      INTEGER I,IPR,K,NATOM,ICOUNT,JCOUNT,FLAG
      CHARACTER CHNAME*1,STRING*10
      LOGICAL HOH,HOHALL
C     ..
C     .. Local Arrays ..
      REAL SUMA(ILINE)
      INTEGER NA(4),NPOINT(MAXATOM)
      CHARACTER RESA(ILINE)*3,SEQ(1000,4)*6,SEQA(ILINE)*6
      CHARACTER TITLE(4)*18,ATOMSTORE(MAXATOM)*10
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C     .. Data statements ..
      DATA TITLE/'of zero           ','less than 5 AngSq ',
     +     'less than 10 AngSq','more than 10 AngSq'/
C     ..
C
C
C----  Introduce yourself
C
      IF (FLAG.EQ.1) THEN
        WRITE(6,FMT='(//,A,/,A,//)')
     +                      ' ANALYSIS OF ACCESSIBLE AREAS OF WATERS',
     +                      ' --------------------------------------'
      ELSE
        WRITE(6,FMT='(/,A,/,A,//)')' Differences for individual waters',
     +                             ' ---------------------------------'
      ENDIF
C
C     Initialise variables + arrays to be zero
C    
      IPR = 0
C
      DO 60 I = 1,4
        NA(I) = 0
 60   CONTINUE
C
C----  Read first residue
C
      ICOUNT = 1
 10   STRING = ATOMSTORE(ICOUNT)
      ACCESS = NPOINT(ICOUNT)/PNTDEN
C
C----  check whether it's water
C
      IF (STRING(1:3).NE.'HOH' .AND. STRING(1:3).NE.'WAT') THEN
C
C----  check whether it has zero area difference
C
        IF ((FLAG.EQ.0 .AND. NPOINT(ICOUNT).EQ.0) .OR. FLAG.EQ.1) THEN
C
          ICOUNT = ICOUNT + 1
          IF (ICOUNT.LE.NATOM) GO TO 10
          CALL CCPERR(1,'WATERAREA: NO WATERS IN OUTPUT')
          RETURN
C
        ENDIF
C
      ENDIF
C
      IPR = IPR + 1
      RESA(IPR) = STRING(1:3)
      CHNAME = STRING(5:5)
      SEQA(IPR) = STRING(5:10)
      SUMA(IPR) = ACCESS
      I = 1
      IF (ACCESS.GT.0) I = 2
      IF (ACCESS.GT.5.0) I = 3
      IF (ACCESS.GT.10) I = 4
      NA(I) = NA(I) + 1
      SEQ(NA(I),I) = SEQA(IPR)
      TOTAL = ACCESS
      CHAIN = ACCESS
C
C
C---- Main loop
C
      DO 30 JCOUNT = (ICOUNT + 1),NATOM
C
        STRING = ATOMSTORE(JCOUNT)
        ACCESS = NPOINT(JCOUNT)/PNTDEN
C
C---- first check that it's water
C
        IF (STRING(1:3).EQ.'HOH' .OR. STRING(1:3).EQ.'WAT') THEN
C
C---- then check if area difference is greater than zero
C
          IF (FLAG.EQ.1 .OR. (FLAG.EQ.0 .AND. NPOINT(JCOUNT).NE.0)) THEN
C
            IF (STRING(5:5).EQ.CHNAME) THEN
              CHAIN = CHAIN + ACCESS
            ELSE
C
C---- new chain
C
            WRITE (6,FMT=6002) (RESA(I),SEQA(I),SUMA(I),I=1,IPR)
 6002       FORMAT (4 (5X,A,1X,A,2X,F8.2))
C
            IF (FLAG.EQ.1) THEN
              call ccp4h_summary_beg()
              WRITE (6,FMT=6000) CHNAME,CHAIN
 6000         FORMAT (/' Total area of chain ',A,': ',F10.1,//)
              call ccp4h_summary_end()
            ELSE
              call ccp4h_summary_beg()
              WRITE (6,FMT=6001) CHNAME,CHAIN
 6001         FORMAT (/' Total area difference of chain ',A,
     +                 ': ',F10.1,//)
              call ccp4h_summary_end()
            ENDIF
C
            CHNAME = STRING(5:5)
            CHAIN = ACCESS
            IPR = 0
          END IF
C
C
          TOTAL = TOTAL + ACCESS
C
          IPR = IPR + 1
          RESA(IPR) = STRING(1:3)
          CHNAME = STRING(5:5)
          SEQA(IPR) = STRING(5:10)
          SUMA(IPR) = ACCESS
C
          IF (IPR.EQ.ILINE) THEN
            WRITE (6,FMT=6002) (RESA(I),SEQA(I),SUMA(I),I=1,ILINE)
            IPR = 0
          END IF
C
C
          I = 1
          IF (ACCESS.GT.0) I = 2
          IF (ACCESS.GT.5.0) I = 3
          IF (ACCESS.GT.10) I = 4
          NA(I) = NA(I) + 1
          SEQ(NA(I),I) = STRING(5:10)
C
          ENDIF
C
        END IF

 30   CONTINUE
C
C----  End of file  - write out final residues
C
   40 WRITE (6,FMT=6002) (RESA(I),SEQA(I),SUMA(I),I=1,IPR)
C
      IF (FLAG.EQ.1) THEN
        call ccp4h_summary_beg()
        WRITE (6,FMT=6000) CHNAME,CHAIN
        WRITE (6,FMT=6004) TOTAL
 6004   FORMAT (/' TOTAL AREA: ',F10.1,/)
        call ccp4h_summary_end()
        WRITE (6,FMT=6006) (NA(I),I=1,4)
 6006   FORMAT (/1X,'Waters sorted by area',
     +          /1X,'---------------------',/,
     +          /1X,'Accessible area          0    .LT. 5   .LT. 10',
     +         '   .GT. 10',/1X,'Number of waters',4I10,/)
      ELSE
        call ccp4h_summary_beg()
        WRITE (6,FMT=6001) CHNAME,CHAIN
        WRITE (6,FMT=6005) TOTAL
 6005   FORMAT (/' TOTAL AREA DIFFERENCE: ',F10.1)
        call ccp4h_summary_end()
        WRITE (6,FMT=6007) (NA(I),I=1,4)
 6007   FORMAT (/1X,'Waters sorted by area differences',
     +          /1X,'---------------------------------',/,
     +          /1X,'Accessible area differences   0    .LT. 5    .LT.',
     +         '10   .GT. 10',/1X,'Number of waters     ',4I10,/)
      ENDIF
C
      IF (HOH .AND. FLAG.EQ.1 .AND. NA(1).GT.0) THEN
        WRITE (6,FMT='(/,2A)') ' Nb: Waters with zero accessible area',
     +  ' are completely enclosed by protein'
      ELSE IF (HOHALL .AND. FLAG.EQ.1 .AND. NA(1).GT.0) THEN
        WRITE (6,FMT='(/,2A,/,A)') ' Nb: Waters with zero accessible',
     +  ' area are completely enclosed by protein',' and/or waters'
      END IF
C
      DO 50 I = 1,4
        IF (FLAG.EQ.1) THEN
          WRITE (6,FMT=6008) TITLE(I), (SEQ(K,I),K=1,NA(I))
 6008     FORMAT (/1X,'Waters with accessible area ',A,/ (1X,12A))
        ELSE
          WRITE (6,FMT=6009) TITLE(I), (SEQ(K,I),K=1,NA(I))
 6009     FORMAT (/1X,'Waters with accessible area differences '
     +                ,A,/ (1X,12A))
        ENDIF
   50 CONTINUE
C
      WRITE(6,FMT='(//)')
C
      RETURN
C
      END
C
C     ================================================================
      SUBROUTINE BRICKSETUP(BINSIZE,HOHRAD,XYZLIM,MINX,MAXX,MINY,MAXY,
     +                      MINZ,MAXZ,MAXNET,MAXBIN,NET,VERB)
C     ================================================================
C
C  Sets up the brick array and checks its size prior to binning the
C  the atoms according to position.
C
      IMPLICIT NONE
C   ..
C   ..Scalar arguments
      REAL BINSIZE,HOHRAD
      INTEGER MINX,MAXX,MINY,MAXY,MINZ,MAXZ,MAXNET,MAXBIN
      LOGICAL VERB
C   ..
C   ..Array arguments
      REAL XYZLIM(2,3)
      INTEGER NET(MAXNET)
C   ..
C   ..Local scalars
      INTEGER ISIZE,ICOUNT
C   ..
C   ..External subroutines
      EXTERNAL CCPERR
C
C----  First zero the bricking array
C
      DO 10 ICOUNT=1,MAXNET
        NET(ICOUNT) = 0
 10   CONTINUE
C
C----  Set size of bricks for searching
C
      BINSIZE = (HOHRAD+2.0)*2.0
C
C      Find bin limits - and expand by two bins
C
      MINX = NINT(XYZLIM(1,1)/BINSIZE) - 2
      MAXX = NINT(XYZLIM(2,1)/BINSIZE) + 2
      MINY = NINT(XYZLIM(1,2)/BINSIZE) - 2
      MAXY = NINT(XYZLIM(2,2)/BINSIZE) + 2
      MINZ = NINT(XYZLIM(1,3)/BINSIZE) - 2
      MAXZ = NINT(XYZLIM(2,3)/BINSIZE) + 2
C
C----  Check size
C
      ISIZE = (MAXX-MINX+1)* (MAXBIN+1)* (MAXY-MINY+1)* (MAXZ-MINZ+1)
C
      IF (VERB) WRITE(6,'(//,A,/,21X,3(3X,2I5),/)') 
     +  ' Dimensions of the net:  MinX,MaxX    MinY,MaxY    MinZ,MaxZ',
     +  MINX, MAXX, MINY, MAXY, MINZ, MAXZ
C
C
      IF (ISIZE.GT.MAXNET) THEN
        WRITE (6,FMT=6028) ISIZE,MAXNET
 6028   FORMAT (' ERROR: Dimension of NET too small in program. Need: ',
     +         /,I10,'  Currently: ',I10)
C
C            *********************
        CALL CCPERR(1, 'Parameter MAXNET too small')
C            *********************
C
      END IF
C
      RETURN
C
      END
C
C     ======================================================
      SUBROUTINE SYMMAT(BRKMAT,IBRKMAT,NSYM,S,TR,TRVEC,VERB)
C     ======================================================
C
C   THIS SUBROUTINE PREPARES THE [S] MATRICES (MAXIMUM 12) THAT
C   WILL BE USED TO GENERATE SYMMETRY RELATED ATOMS .
C   IT COMBINES THE MATRIX FROM A BROOKHAVEN FILE WITH
C   SYMMETRIES SUPPLIED ON INPUT CARDS.
C
C   NB: also generates 27 translation vectors which contain
C   all linear combinations of shifts by one unit cell vector,
C   including the "identity translation" ie. zero vector.
C
C Arguments
C  
C INPUT
C   BRKMAT    Scale matrix from PDB input
C   IBRKMAT   inverse scale matrix
C   NSYM      No of symmetry operations
C   TR        Fractional transformation matrices
C   VERB      Verbose output flag
C
C OUTPUT
C   S         Orthogonal transformation matrices        
C   TRVEC     Translation vectors
C
      IMPLICIT NONE
C
C     .. Scalar arguments ..
      INTEGER NSYM
      LOGICAL VERB
C     ..
C     .. Array arguments ..
      REAL BRKMAT(4,4),IBRKMAT(4,4),S(4,4,12),TR(4,4,12),TRVEC(3,125)
C     ..
C     .. Local Scalars ..
      INTEGER I,I1,I2,J,K,NVEC
      LOGICAL IDENT
C     ..
C     .. Local Arrays ..
      REAL MAT1(4,4),MAT2(4,4),MAT3(4,4),TEMP(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL MATMUL4,CCPERR
C     ..
C     .. Intrinsic functions
      INTRINSIC REAL
C
C---- symmetry operations from SYMMETRY keyword(s)
C
C
      DO 10 I = 1,10
        S(4,4,I) = 1.0
        TR(4,4,I) = 1.0
   10 CONTINUE
C
      IF (NSYM.GT.12) THEN
        WRITE (6,FMT=6016) NSYM
 6016   FORMAT (' *** NSYM = ',I3,' MAXIMUM ALLOWED IS 12')
C
C              *********************
        CALL CCPERR(1, 'Too many symmetry operations')
C              *********************
C
      END IF
C
C --- Write fractional matrices to standard output
C
C --- (nb only if there are more than none ...)
C
      IF (NSYM.GT.0) THEN
        K = 1
C
        IF (VERB) WRITE (6,FMT=6020)
 6020   FORMAT (/,' SYMMETRY MATRICES FROM SYMMETRY KEYWORD(S)',/)
C
C --- Check for the identity matrix
C
 170    CONTINUE
C
        IDENT = .TRUE.
C
        DO 710 I = 1,4
C
          DO 720 J = 1,4
C
             IF (I.EQ.J .AND. TR(I,J,K).NE.1) IDENT = .FALSE.
             IF (I.NE.J .AND. TR(I,J,K).NE.0) IDENT = .FALSE.
C
 720      CONTINUE
C
 710    CONTINUE
C
C
        IF (IDENT) THEN
C --- Identity matrix found - overwrite from top of stack
C
          DO 740 I =1,4
            DO 730 J = 1,4
              TR(I,J,K) = TR(I,J,NSYM)
 730        CONTINUE
 740      CONTINUE
C
          NSYM = NSYM - 1
C
C --- Check if we've run out of symmetry operations ...
C
          IF (NSYM .EQ. 0) THEN
             CALL CCPERR (2,' POSSIBLE ERROR IN SYMMETRY OPERATIONS')
             WRITE (6,FMT=6100)
 6100        FORMAT (' There are no non-identity symmetry operations',/
     +               ' This is equivalent to using P1 symmetry',/)
          ELSE
C
C --- Repeat check for this matrix
C
             GO TO 170
          ENDIF
C
        ENDIF
C
C --- Write to standard output
C
        IF (VERB) THEN
C
        IF (K .LE. NSYM) THEN
          WRITE (6,FMT=6022) K, ((TR(I1,I2,K),I2=1,4),I1=1,4)
 6022     FORMAT (1X,'Fractional matrix ',I2,
     +    /,4 (1X,3F9.3,4X,F9.3,/))
C
          K = K + 1
          IF (K.LE.NSYM) GO TO 170
        END IF
C
        END IF
C
      END IF
C
C --- Generate the symmetry matrices in orthogonal coords
C
C --- (only if the number of matrices is non-zero)
C
      IF (NSYM.GT.0) THEN
C
C---    get the correct order for
C       matrix multiplications:  [S] = [IBRKMAT] * [TR] * [BRKMAT]
C
      IF (VERB) WRITE (6,FMT=6027)
 6027 FORMAT (/,
     + 'Final orthogonal matrices derived using matrix multiplication:',
     + //,
     + ' [Orth.Symm.] = [Orth.Scale]*[Frac.Symm.]*[Frac.Scale]',/)
C
C
      DO 90 I = 1,NSYM
        DO 60 J = 1,4
          DO 50 K = 1,4
            MAT2(K,J) = TR(K,J,I)
   50     CONTINUE
   60   CONTINUE
C
C
        CALL MATMUL4(MAT1,IBRKMAT,MAT2)
        CALL MATMUL4(MAT3,MAT1,BRKMAT)
C
C --- Write to standard output
C
        IF (VERB) WRITE (6,FMT=6026) I, ((MAT3(I1,I2),I2=1,4),I1=1,4)
 6026   FORMAT (1X,'Orthogonal matrix ',I2,
     +  /,4 (1X,3F9.3,4X,F9.3,/))
C
        DO 80 J = 1,4
          DO 70 K = 1,4
            S(K,J,I) = MAT3(K,J)
   70     CONTINUE
   80   CONTINUE
C
C
   90 CONTINUE
C
      ENDIF
C
C --- Now generate the translation vectors
C
C --- set up zero vector:
C
      TRVEC(1,1) = 0.0
      TRVEC(2,1) = 0.0
      TRVEC(3,1) = 0.0
C
C --- first fractional coordinates ...
      NVEC = 1
      DO 310 I = -2,2
        DO 320 J = -2,2
          DO 330 K = -2,2
C
C ---       don't regenerate the zero vector:
            IF ((I.EQ.0) .AND. ((J.EQ.0) .AND. (K.EQ.0))) GO TO 330
C
            NVEC = NVEC + 1
            TRVEC(1,NVEC) = REAL(I)
            TRVEC(2,NVEC) = REAL(J)
            TRVEC(3,NVEC) = REAL(K)
C
 330      CONTINUE
 320    CONTINUE
 310  CONTINUE
C
C --- V. IMPORTANT!!! Now convert to orthogonal coordinates:
C
      DO 340 I = 1,NVEC
C
        DO 350 J = 1,3
          TEMP(J) = IBRKMAT(J,1)*TRVEC(1,I) +
     +              IBRKMAT(J,2)*TRVEC(2,I) +
     +              IBRKMAT(J,3)*TRVEC(3,I)
 350    CONTINUE
C
        DO 360 J = 1,3
          TRVEC(J,I) = TEMP(J)
 360    CONTINUE
C
 340  CONTINUE
C
C
  100 RETURN
      END
C
C     ==================================================================
      SUBROUTINE SYMGEN(NATOM,XYZ,XYZLIM,DLIM,ITYPE,RESTYP,NSYM,S,TRVEC,
     +                  NOTRANS,VERB)
C     ==================================================================
C
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
C     ..
C     .. Scalar Arguments ..
      INTEGER NATOM,NSYM
      REAL DLIM
      LOGICAL NOTRANS,VERB
C     ..
C     .. Array arguments ..
      REAL S(4,4,12),XYZLIM(2,3),XYZ(3,MAXATOM),TRVEC(3,125)
      INTEGER ITYPE(MAXATOM)
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     .. Local Scalars ..
      INTEGER I,ICOUNT,J,K,L,NTR,NTRMAX
C     ..
C     .. Local Arrays ..
      REAL X(3),XB2(2,3),XT(3)
C     ..
C     ..
C
      ICOUNT = NATOM
C
C
      DO 10 K = 1,3
        XB2(1,K) = XYZLIM(1,K) - DLIM
        XB2(2,K) = XYZLIM(2,K) + DLIM
   10 CONTINUE
C
C --- LOOP over all atoms
      DO 50 I = 1,ICOUNT
C
C ---   LOOP over all symmetry ops. (treat op zero [J=0] as identity for
C                                    purpose of lattice translations)
        DO 40 J = 0,NSYM
C
          IF (J.EQ.0) THEN
            DO 25 K = 1,3
              X(K) = XYZ(K,I)
   25       CONTINUE
          ELSE
            DO 20 K = 1,3
              X(K) = S(K,1,J)*XYZ(1,I) + S(K,2,J)*XYZ(2,I) +
     +               S(K,3,J)*XYZ(3,I) + S(K,4,J)
   20       CONTINUE
          ENDIF
C
C ---     NOTRANS mode: switch off the all the lattice translations
C         except the zero vector (NTR=1).
C         Otherwise use all 27 vectors.
C         PJX 02/02/01: Updated to use -2 to +2 lattice translations
C         in each direction - this gives 125 vectors now
C
          IF (NOTRANS) THEN
            NTRMAX = 1
          ELSE
            NTRMAX = 125
          ENDIF
C
C ---     Do the translations and test each one
C
          DO 60 NTR = 1,NTRMAX
C
C ---       Ignore the zero vector if J=0
            IF (J.EQ.0 .AND. NTR.EQ.1) GO TO 60
C
C ---       Otherwise, use the zero vector too
            DO 70 L = 1,3
              XT(L) = X(L) + TRVEC(L,NTR)
              IF (XT(L).LT.XB2(1,L) .OR. XT(L).GT.XB2(2,L)) GO TO 60 
  70        CONTINUE
C
C ---       Symmetry-generated atom lies in the box: store attributes
C
            NATOM = NATOM + 1
            DO 30 L = 1,3
              XYZ(L,NATOM) = XT(L)
   30       CONTINUE
            ITYPE(NATOM) = ITYPE(I)
            RESTYP(NATOM) = RESTYP(I)
C
   60     CONTINUE
C
   40   CONTINUE
   50 CONTINUE
C
C
      IF (VERB )THEN
        WRITE (6,FMT=
     +  '(//,'' NO. OF REAL ATOMS :'',I6)') ICOUNT
        WRITE (6,FMT=
     +  '(/,'' no. of symmetry related atoms in all bricks :'',I6)')
     +   (NATOM - ICOUNT)
      END IF
C
      RETURN
      END
C
C     ===============================================================
      SUBROUTINE BINATOMS(X,NATOM,NET,MINX,MAXX,MINY,MAXY,MINZ,MAXZ,
     +                    MAXBIN,BINSIZE,RESTYP,HOHALL,VERB)
C     ===============================================================
C
C----  Put atoms in bins
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MAXATOM
      PARAMETER (MAXATOM=60000)
C     ..
C     .. Scalar Arguments ..
      REAL BINSIZE
      INTEGER MAXX,MAXY,MAXZ,MINX,MINY,MINZ,NATOM,MAXBIN
      LOGICAL HOHALL,VERB
C     ..
C     .. Array Arguments ..
      REAL X(3,MAXATOM)
      INTEGER NET(0:MAXBIN,MINX:MAXX,MINY:MAXY,MINZ:MAXZ)
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     .. Local Scalars ..
      INTEGER I,IN,INMAX,IX,IY,IZ
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC NINT
C     ..
C     ..
C
C pete DUNTEN: need to initialize inmax
      INMAX = 0
      DO 10 I = 1,NATOM
        IF (RESTYP(I) .OR. HOHALL) THEN
          IX = NINT(X(1,I)/BINSIZE)
          IY = NINT(X(2,I)/BINSIZE)
          IZ = NINT(X(3,I)/BINSIZE)
          NET(0,IX,IY,IZ) = NET(0,IX,IY,IZ) + 1
          IN = NET(0,IX,IY,IZ)
          IF (IN.LE.MAXBIN) NET(IN,IX,IY,IZ) = I
          IF (IN.GT.INMAX) INMAX = IN
        END IF
C
   10 CONTINUE
C
      IF (VERB) WRITE (6,FMT='(/'' Max atoms in bin: '',I5)') INMAX
C
      IF (INMAX.GT.MAXBIN) THEN
        WRITE (6,FMT=6000) INMAX,BINSIZE,MAXBIN,INMAX
 6000   FORMAT (//1X,'WARNING: there are ',I6,' atoms in a',
     +       ' bin of width ',F6.1,' A'/1X,
     +       'Currently a maximum of ',I6,' is allowed in each bin',
     +       //1X,'Check whether this number is reasonable - there',
     +       ' maybe an error',/1X,'in your input symmetry matrices.',
     +       //1X,'Otherwise, reset the MAXBIN parameter to ',I6,
     +       ' and recompile.',/1X)
        CALL CCPERR(1,'PARAMETER MAXBIN EXCEEDED')
      END IF
C
      RETURN
      END
C
C     =================================================================
      SUBROUTINE AREA(X,ITYPE,RADIUS,NATOMTYP,PNTDEN,NPOINT,NET,
     +                MINX,MAXX,MINY,MAXY,MINZ,MAXZ,MAXBIN,MAXATOM,
     +                BINSIZE,RESTYP,NOINTYP,NREADIN,HOH,HOHALL,
     +                SELECT_CALC)
C     =================================================================
C
C      This routine calculates the solvent accessible area of a molecule
C  using the Lee and Richards (Mike Connolly?) algorithm.
C
C  On entry:
C      X contains the atomic coords (in Angstrom orthog axes)
C      NREADIN - number of "real" atoms in X
C      ITYPE  - the atom type for each atom
C      RADIUS contains radius for each atom type.  This radius should be
C            the sum of the Van der Waals and solvent molecule.
C      NATOMTYP - number of atom types/radii
C      PNTDEN - Point density - dots per square angstrom
C      NOINTYP - number of atoms of each type
C      NET - grid with atoms listed in each spatial cell
C      MINX,MAXX,MINY,MAXY,MINZ,MAXZ - limits of the grid in NET
C      MAXATOM - maximum no. of atoms
C      BINSIZE - physical length of cubic cell edge in NET
C      RESTYP - whether atom is protein or water
C      NREADIN - no of "actual" ie, not symm-related atoms
C      HOH, HOHALL - flags for different MODES
C
C  On exit:
C      NPOINT - surface points for each atom - one point for each
C             1/PNTDEN Angstrom squared.
C
C
C    Max number neighbours       MXNNBR
C    Max number points on sphere MXPNT
C    Neighbours                  XYZNNB
C    coord. of points on sphere  XYZSPH
C    number of points sphere     NPTSPH
C
C
      IMPLICIT NONE
C
C     .. Parameters ..
      INTEGER MXNNBR
      PARAMETER (MXNNBR=2000)
      INTEGER MXPNT
      PARAMETER (MXPNT=15000)
C     ..
C     .. Scalar Arguments ..
      REAL BINSIZE,PNTDEN
      INTEGER NATOMTYP,MAXX,MAXY,MAXZ,MINX,MINY,MINZ,MAXBIN
      INTEGER NREADIN,MAXATOM
      LOGICAL HOH,HOHALL
C     ..
C     .. Array Arguments ..
      REAL RADIUS(0:NATOMTYP),X(3,MAXATOM)
      INTEGER ITYPE(MAXATOM),NPOINT(MAXATOM),NOINTYP(0:NATOMTYP)
      INTEGER NET(0:MAXBIN,MINX:MAXX,MINY:MAXY,MINZ:MAXZ)
      INTEGER SELECT_CALC(MAXATOM)
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     .. Local Scalars ..
      REAL DIJ2,DIJJ2,DX,DY,DZ,EQUAT,FI,FJ,HOR,PI,POINTS,RI,RJ,RJJ,
     +     SAREA,SUMR,TRUNC,VERT,XY,Z
      INTEGER I,II,IPNT,ITI,IX,IXX,IY,IYY,IZ,IZZ,J,JJ,KK,NHOR,NNBR,NPT,
     +        NVERT
C     ..
C     .. Local Arrays ..
      REAL CW(3),RNNB(MXNNBR),XYZNNB(3,MXNNBR),XYZSPH(3,MXPNT,0:20)
      INTEGER NPTSPH(0:20)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN,COS,NINT,REAL,SIN,SQRT
C     ..
C     ..
C
      PI = ATAN(1.0)*4.0
C
C ---  Build sphere around 0,0,0 for each of the atom types/radii
C
      DO 30 I = 0,NATOMTYP
C
C---- only build if atoms of this type have been found
C
       IF (NOINTYP(I).GT.0) THEN
C
        RI = RADIUS(I)
        SAREA = 4*PI* (RI*RI)
C
C---- number of points on sphere
C
        POINTS = SAREA*PNTDEN
        EQUAT = SQRT(POINTS*PI)
        VERT = EQUAT/2.0
        NVERT = NINT(VERT)
C
C----  Counter for points found on sphere
C
        NPT = 0
        TRUNC = 0.0
C
        DO 20 II = 0,NVERT
          FI = REAL(II)*PI/VERT
          HOR = SIN(FI)*EQUAT - TRUNC
          NHOR = NINT(HOR)
          IF (II.EQ.0) NHOR = 1
          TRUNC = REAL(NHOR) - HOR
          XY = SIN(FI)*RI
          Z = COS(FI)*RI
C
C
          DO 10 JJ = 0,NHOR - 1
            FJ = (2.0*PI*REAL(JJ))/REAL(NHOR)
            NPT = NPT + 1
            XYZSPH(1,NPT,I) = COS(FJ)*XY
            XYZSPH(2,NPT,I) = SIN(FJ)*XY
            XYZSPH(3,NPT,I) = Z
   10     CONTINUE
   20   CONTINUE
C
        IF (NPT.GT.MXPNT) THEN
          WRITE (6,FMT=1010) I,PNTDEN,INT(POINTS),MXPNT
C                         *******************************************
          CALL CCPERR(1,' Too many points: increase parameter MXPNT')
C                         *******************************************
 1010     FORMAT(/,'* Atom type ',I2,': with PNTDEN = ',F8.2,
     +           /,'  requires (approx.) ',I8,' points',
     +           /,'  but only has MXPNT = ',I8,' available',/)
        ENDIF
C
        NPTSPH(I) = NPT
C
C----  End of sphere building loop
C
       ENDIF
C
   30 CONTINUE
C
C
C  *********      MAIN LOOP OVER ALL ATOMS    **********
C
      DO 120 I = 1,NREADIN
C
C----  Initialise number of surface points to zero
C
        NPOINT(I) = 0
C
C----  Only calculate area for selected atoms
C
        IF (SELECT_CALC(I).EQ.1) THEN
C
C----  Reject protein atoms if only considering waters
C
          IF ((HOH.OR.HOHALL) .AND. RESTYP(I)) GO TO 120
          ITI = ITYPE(I)
          RI = RADIUS(ITI)
C     
C----  Find bin for this atom
C
          IXX = NINT(X(1,I)/BINSIZE)
          IYY = NINT(X(2,I)/BINSIZE)
          IZZ = NINT(X(3,I)/BINSIZE)
C
C---- FIND THE NEIGHBOURS
C      within the sum of Van der Waal + twice solvent radius
C  by looping over the 27 local bins

          NNBR = 0
C
          DO 80 IZ = IZZ - 1,IZZ + 1
            DO 70 IY = IYY - 1,IYY + 1
              DO 60 IX = IXX - 1,IXX + 1
C
C----  Now loop over atoms in this bin - there can be zero atoms in bin
C
                DO 50 IPNT = 1,NET(0,IX,IY,IZ)
                  J = NET(IPNT,IX,IY,IZ)
C                 Reject atom if same as source atom or if
C                 specifically unselected
                  IF (I.EQ.J .OR. SELECT_CALC(J).EQ.0) GO TO 50
                  RJ = RADIUS(ITYPE(J))
                  SUMR = RI + RJ
                  DX = X(1,I) - X(1,J)
                  IF (ABS(DX).GT.SUMR) GO TO 50
                  DY = X(2,I) - X(2,J)
                  IF (ABS(DY).GT.SUMR) GO TO 50
                  DZ = X(3,I) - X(3,J)
                  IF (ABS(DZ).GT.SUMR) GO TO 50
                  DIJ2 = DX*DX + DY*DY + DZ*DZ
                  IF (DIJ2.GE.SUMR*SUMR) GO TO 50
                  NNBR = NNBR + 1
C
                  IF (NNBR.GT.MXNNBR) THEN
                    WRITE (6,FMT=6000) MXNNBR,I,ITI,RI
 6000               FORMAT (//' ** More than ',I4,' neigbours for atom',
     +                   I5,' type',I4,' radius',F6.3,' - redimension ',
     +                   'parameter MXNNBR in s/r AREA')
                    CALL CCPERR(1,'AREA: PARAMETER MXNNBR TOO SMALL')
                  END IF
C
                  DO 40 KK = 1,3
                    XYZNNB(KK,NNBR) = X(KK,J)
 40               CONTINUE
C
                  RNNB(NNBR) = RJ
 50             CONTINUE
 60           CONTINUE
 70         CONTINUE
 80       CONTINUE
C
C --- Stick the sphere around the atom
C
C    Loop over points on sphere
C
          DO 110 II = 1,NPTSPH(ITI)
C
            DO 90 KK = 1,3
              CW(KK) = X(KK,I) + XYZSPH(KK,II,ITI)
 90         CONTINUE
C
C    Find the accessible area
C
C--- Test for neighbours that lie within RJJ of this point
C   Loop over neighbours
C
            DO 100 JJ = 1,NNBR
              RJJ = RNNB(JJ)
              DX = CW(1) - XYZNNB(1,JJ)
              IF (ABS(DX).GT.RJJ) GO TO 100
              DY = CW(2) - XYZNNB(2,JJ)
              IF (ABS(DY).GT.RJJ) GO TO 100
              DZ = CW(3) - XYZNNB(3,JJ)
              IF (ABS(DZ).GT.RJJ) GO TO 100
              DIJJ2 = DX*DX + DY*DY + DZ*DZ
C
C----  Skip to next point on sphere
C
              IF (DIJJ2.LT.RJJ*RJJ) GO TO 110
C
 100        CONTINUE
C
C  We have a surface point as no neighbour sufficiently close
C
            NPOINT(I) = NPOINT(I) + 1
C
C----  End of points on sphere loop
C
 110      CONTINUE
C
C---- End of if-endif for selected atoms
C
        ENDIF
C
C---- End of main atom loop
C
 120  CONTINUE
C
      RETURN
      END
C
C
C     ===================================================
      SUBROUTINE SETUP_SELECTION(MAXATOM,SELECT_CALC)
C     ===================================================
C
C---- Set up the selection array
C     Each element of the array corresponds to an atom
C     The element can be
C          1 (include in all calculations)
C          0 (exclude from all calculations)
C         -1 (include in calculations for other atoms
C             don't calculate for this one)
C
C     Currently simply selects all atoms for everything
C
      IMPLICIT NONE
C
      INTEGER MAXATOM,I,SELECT_CALC(MAXATOM)
C
      DO I = 1,MAXATOM
        SELECT_CALC(I) = 1
      END DO
      RETURN
      END
C
C     ===================================================
      SUBROUTINE SETRADII(MAXTYPE,NATOMTYP,HOHRAD,RADIUS)
C     ===================================================
C
C---- Increment radii by solvent radius
C
      IMPLICIT NONE
C
      INTEGER MAXTYPE,NATOMTYP,IC
      REAL HOHRAD,RADIUS(0:MAXTYPE)
C
      DO 100 IC = 0,NATOMTYP
        RADIUS(IC) = RADIUS(IC) + HOHRAD
  100 CONTINUE
      RETURN
      END
C
C     ==============================================================
      SUBROUTINE SEARCHISOLATED(SELECT_CALC,NATOMS,X,RADIUS,PROBE_RAD,
     +           NPT,ITYPE,NET,BINSIZE,MAXATOM,NATOMTYP,MINX,MAXX,MINY,
     +           MAXY,MINZ,MAXZ,MAXBIN,PNTDEN,NSURFACES,SURFACE,IMODE)
C     ==============================================================
C
C---- Search for closed/independent surfaces within the structure
C     Alternatively, search for clusters of buried atoms
C
C     IMODE tells the routine how to operate
C     IMODE = 1: cluster atoms with non-zero ASA
C                In this case the radii of the atoms must be the
C                expanded VdW radii
C     IMODE = 2: cluster buried atoms
C                In this case the radii of the atoms must be the
C                "bare" VdW radii
C
      IMPLICIT NONE
C
C     ..Parameters
C     MAXSURF is the maximum number of isolated surfaces that
C     can be found
      INTEGER MAXSURF
      PARAMETER (MAXSURF=80)
C
C     ..Scalar Arguments
      INTEGER NATOMS,MAXATOM,NATOMTYP,NSURFACES,IMODE
      INTEGER MAXBIN,MINX,MAXX,MINY,MAXY,MINZ,MAXZ
      REAL    BINSIZE,PNTDEN,PROBE_RAD
C
C     ..Array Arguments
      INTEGER NPT(MAXATOM),ITYPE(MAXATOM)
      INTEGER NET(0:MAXBIN,MINX:MAXX,MINY:MAXY,MINZ:MAXZ)
      INTEGER SURFACE(MAXATOM),SELECT_CALC(MAXATOM)
      REAL X(3,MAXATOM),RADIUS(0:NATOMTYP)
C
C     ..Local scalars
      INTEGER I,J,K,K1,IX,IY,IZ,IXX,IYY,IZZ,IPNT
      integer pjxcheck,pjxtotal,pjxburied
      INTEGER THISSURF
      REAL    RI,RJ,SUMR,DX,DY,DZ,DIJ2,RAD_FUZZ
      LOGICAL LCHECK
      logical diag
C
C     ..Local arrays
      INTEGER SURFLIST(0:MAXSURF),TAG(MAXSURF)
cpjx
cpjx  Diagnositics
      diag = .false.
C
C---- Check for valid IMODE
C
      IF (IMODE.NE.1 .AND. IMODE.NE.2) THEN
        CALL CCPERR(1,'SEARCHISOLATED: Bad IMODE')
      END IF
C
C---- Initialise surface array
C
      NSURFACES = 0
      DO 90 I=1,NATOMS
        SURFACE(I) = 0
 90   CONTINUE
C
C---- "Fuzz" factor for atom radii in IMODE 2
C
      RAD_FUZZ = 1.0
C
C---- Loop over atoms
C
      DO 10 I=1,NATOMS
C
        if (diag) then
          write(6,*)'pjx: ---------------------------------------------'
          write(6,*)'pjx: atom ',i,' has ',npt(i),' surface points'
        end if
C
C---- Use this atom? Depends on the mode
C     IMODE 1 checks atoms with non-zero ASA
C     IMODE 2 checks atoms with zero ASA (buried)
        LCHECK = .FALSE.
C
C     Only consider selected atoms
        IF (SELECT_CALC(I).EQ.1) THEN
          IF (IMODE.EQ.1 .AND. NPT(I).GT.0) LCHECK = .TRUE.
          IF (IMODE.EQ.2 .AND. NPT(I).EQ.0) LCHECK = .TRUE.
        END IF
C
        IF (LCHECK) THEN
C
C----  Find bin for this atom
C
          IXX = NINT(X(1,I)/BINSIZE)
          IYY = NINT(X(2,I)/BINSIZE)
          IZZ = NINT(X(3,I)/BINSIZE)
          if (diag) write(6,*)'pjx: binned in ',ixx,iyy,izz
C
C---- Assign radius
C     The values stored in RADIUS are the expanded VdW radii
          IF (IMODE.EQ.1) THEN
            RI = RADIUS(ITYPE(I))
          ELSE
            RI = RADIUS(ITYPE(I)) - PROBE_RAD + RAD_FUZZ
          END IF
          if (diag) write(6,*)'pjx: radius is ',ri
C
C---- Examine the neighbours
C     SURFLIST is a list of surfaces found on neighbouring atoms
C     SURLIST(0) holds the number of such surfaces
C
          SURFLIST(0) = 0
C
C---- Loop over neighbouring bins
          DO 80 IZ = IZZ - 1,IZZ + 1
            DO 70 IY = IYY - 1,IYY + 1
              DO 60 IX = IXX - 1,IXX + 1
C
C---- Now loop over atoms in the bin
                if (diag) then
                  write(6,*)'pjx: in neighbouring brick ',ix,iy,iz
                  write(6,*)'pjx: there are ',NET(0,IX,IY,IZ),' atoms'
                end if
                DO 50 IPNT = 1,NET(0,IX,IY,IZ)
                  J = NET(IPNT,IX,IY,IZ)
                  if (diag) write(6,*)'pjx: examining atom ',j
C
C---- Reject if it's the same atom
                  IF (I.EQ.J) GO TO 50
C
C---- Only consider other selected atoms with(out) surface points (depending
C     on IMODE)
C 
                  if (diag) write(6,*)'pjx: neighbour ',j,' has ',
     +                 npt(j),' points'
                  IF (SELECT_CALC(J).EQ.1 .AND.
     +                 ((IMODE.EQ.1 .AND. NPT(J).GT.0) .OR.
     +                 (IMODE.EQ.2 .AND. NPT(J).EQ.0))) THEN
C     The values stored in RADIUS are the expanded VdW radii
                    IF (IMODE.EQ.1) THEN
                      RJ = RADIUS(ITYPE(J))
                    ELSE
                      RJ = RADIUS(ITYPE(J)) - PROBE_RAD + RAD_FUZZ
                    END IF
                    if (diag) write(6,*)'pjx: and its radius is ',ri
                    SUMR = RI + RJ
                    if (diag) then
                      write(6,*)'pjx: total radius is ',sumr
                      write(6,*)'pjx: maximum square separation is ',
     +                     sumr*sumr
                    end if
                    DX = X(1,I) - X(1,J)
                    DY = X(2,I) - X(2,J)
                    DZ = X(3,I) - X(3,J)
                    DIJ2 = DX*DX + DY*DY + DZ*DZ
                    if (diag) write(6,*)'pjx: distance^2 = ',dij2
C
C---- Reject those not in direct contact
                    IF (DIJ2.GE.SUMR*SUMR) GO TO 50
C
C---- Check if a surface has already been assigned
C     Make a list of the neighbouring surfaces
                    IF (SURFACE(J).GT.0) THEN
                      if (diag) then
                        write(6,*)'pjx: neighbour already assigned to ',
     +                       surface(j)
                        write(6,*)'pjx: there are ',SURFLIST(0),
     +                       ' neighbouring surfaces'                  
                      end if
C
C---- Check if we've already seen this surface in the list
C     of neighbouring partial surfaces
                      DO 110 K=1,SURFLIST(0)
                        if (diag) write(6,*)'pjx: neighbouring surface '
     +                       ,k,' is ',surflist(k)
                        IF (SURFACE(J).EQ.SURFLIST(K)) GO TO 120
 110                  CONTINUE
                      if (diag)
     +                     write(6,*)'pjx: not already in current list'
                      SURFLIST(0) = SURFLIST(0) + 1
                      IF (SURFLIST(0).GT.MAXSURF) THEN
                        WRITE (6,FMT='(1X,A,A,I10,A)')
     +                       'In s/r SEARCHISOLATED: number of partial',
     +                       ' surfaces exceeds MAXSURF (',
     +                       MAXSURF,')'
                        CALL CCPERR
     +                  (1,' S/R SEARCHISOLATED: increase MAXSURF')
                      END IF
                      if (diag) write(6,*)'pjx: added surface ',
     +                     surface(j),' as ',surflist(0),
     +                     'th neighbouring surface'
                      SURFLIST(SURFLIST(0)) = SURFACE(J)
 120                  CONTINUE
                    ELSE
                      if (diag) write(6,*)'pjx: this neighbour is not',
     +                     ' yet assigned'
                    END IF
C
                  END IF
C     
 50             CONTINUE
 60           CONTINUE
 70         CONTINUE
 80       CONTINUE
C
C---- Are there any neighboring surfaces?
C
C---- No neighbouring surfaces: start new surface
          IF (SURFLIST(0).EQ.0) THEN
            NSURFACES  = NSURFACES + 1
            SURFACE(I) = NSURFACES
            if (diag) write(6,*)'pjx: new surface! assigned to surface '
     +           ,nsurfaces
C
C---- 1 neighbouring surface: join this atom on
          ELSE IF (SURFLIST(0).EQ.1) THEN
            SURFACE(I) = SURFLIST(1)
            if (diag) write(6,*)'pjx: joined to existing surface ',
     +           surflist(1)
C
C---- More than one neighbouring surface (gets complicated):
          ELSE IF (SURFLIST(0).GT.1) THEN
C---- Pick the lowest surface tag in the list
            if (diag) write(6,*)'pjx: looking for lowest value ',
     +           'neighbouring surface'
            THISSURF = SURFLIST(1)
            if (diag) write(6,*)'pjx: initially this is ',thissurf
            DO K = 2,SURFLIST(0)
              if (diag) write(6,*)'pjx: comparing with ',k,'th surface',
     +              '(',surflist(k),')'
              IF (SURFLIST(K).LT.THISSURF) THEN
                THISSURF = SURFLIST(K)
                if (diag) write(6,*)'pjx: reassigned this surface to ',
     +               thissurf
              END IF
            END DO
C---- Assign current atom to this surface
            SURFACE(I) = THISSURF
            if (diag) write(6,*)'pjx: all connected surfaces will be ',
     +           'reassigned to',thissurf
C---- Now you know that all the surfaces in this list are
C     connected - sweep through all atoms so far and reset
C     any surface assignments to join the surfaces
C     This could make the process quite slow!
            if (diag) write(6,*)'pjx: joining many surfaces'
            DO 150 K = 1,I-1
              IF (SURFACE(K).GT.0) THEN
                if (diag) then
                  write(6,*)'pjx: atom ',k,' already assigned to ',
     +                 'surface ',surface(k)
                  write(6,*)'pjx: list of neighbouring surfaces: ',
     +                 (surflist(k1),k1=1,surflist(0))
                end if
                DO 160 K1 = 2,SURFLIST(0)
                  IF (SURFACE(K).EQ.SURFLIST(K1)) THEN
                    SURFACE(K) = THISSURF
                    if (diag) write(6,*)'pjx: atom ',k,' reassigned to '
     +                   ,thissurf
                  END IF
 160            CONTINUE
              END IF
              if (diag) then
                if (surface(k).ne.thissurf) then
                  write(6,*)'pjx: atom ',k,' was not reassigned'
                end if
              end if
 150        CONTINUE
C
          END IF
C     
        END IF
C
 10   CONTINUE
      if (diag) write(6,*)'pjx: finished initial sweep over atoms'
C
C---- 2nd sweep to re-tag the surfaces from 1...Nsurfaces
C
      NSURFACES = 0
      DO I = 1,MAXSURF
        TAG(I) = -1
      END DO
C
      DO I = 1,NATOMS
C
        IF (SURFACE(I).GT.0) THEN
          DO J = 1,NSURFACES
            IF (SURFACE(I).EQ.TAG(J)) THEN
              SURFACE(I) = J
              GO TO 20
            END IF
          END DO
          IF (NSURFACES.GE.MAXSURF) THEN
            WRITE(6,FMT='(1X,A,I3,A)')
     +           ' S/R SEARCHISOLATED: maximum no of surfaces (',
     +           MAXSURF,') exceeded - increase parameter MAXSURF'
            NSURFACES = -NSURFACES
            RETURN
          END IF
          NSURFACES = NSURFACES + 1
          TAG(NSURFACES) = SURFACE(I)
          SURFACE(I) = NSURFACES
 20       CONTINUE
C
        END IF
C
      END DO
C
C---- Another check
      if (diag) then
        write(6,*)'pjx: checking for consistency'
        pjxburied = 0
        do j = 1,natoms
          if (npt(j).eq.0) then
            pjxburied = pjxburied + 1
            write(6,*)'pjx: atom ',j,' (surface ',surface(j),') ',
     +           'is buried (npt = ',npt(j),')'
          end if
        end do
        pjxtotal = 0
        do i = 0,nsurfaces
          pjxcheck = 0
          do j = 1,natoms
            if (surface(j).eq.i) then
              pjxcheck = pjxcheck + 1
              pjxtotal = pjxtotal + 1
            end if
          end do
          write(6,*)'pjx: CHECK! there are ',pjxcheck,
     +         ' atoms in surface',i
          if (pjxcheck.eq.0) write(6,*)'pjx: WARNING!!!! no atoms in ',
     +         'surface ',i
        end do
        write(6,*)'pjx: CHECK! there are ',pjxtotal,
     +       ' atoms in all surfaces - there should be ',natoms
        write(6,*)'pjx: of which ',pjxburied,' have no asa'
        write(6,*)'pjx: done consistency check'
      end if
C
C---- Finished
C
      RETURN
      END
C
C     ================================================================
      SUBROUTINE ADDATOMTYPE(MAXTYPE,ATYPE,RADIUS,ATOMN,ATMCHECK,
     +     ATOMNAME,VDWR,ATNUM,NATOMTYP)
C     ================================================================
C
C Add an atom type to the list of recognised atoms
C MAXTYPE is the maximum number of allowed types, NATOMTYP the current
C number actually stored. ATYPE, RADIUS, ATOMN are the storage arrays
C for the atom name, VdW radius and atomic number respectively, and
C ATOMNAME, VDWR and ATNUM are the values to be added for the new
C type.
C ATMCHECK 
      IMPLICIT NONE
C
C     ..Scalar arguments
      INTEGER MAXTYPE,NATOMTYP,ATNUM
      REAL    VDWR
      CHARACTER ATMCHECK*20,ATOMNAME*2
C
C     ..Array arguments
      INTEGER ATOMN(0:MAXTYPE)
      REAL    RADIUS(0:MAXTYPE)
      CHARACTER ATYPE(0:MAXTYPE)*2

      NATOMTYP = NATOMTYP + 1
      IF (NATOMTYP .GT. MAXTYPE) THEN
        CALL CCPERR (2,
     +       'S/R ADDATOMTYPE: too many user-defined atom types')
        NATOMTYP = NATOMTYP - 1
        RETURN
      ENDIF
      ATYPE(NATOMTYP) = ATOMNAME
      RADIUS(NATOMTYP) = VDWR
      ATOMN(NATOMTYP) = ATNUM
      ATMCHECK(NATOMTYP:NATOMTYP) = '*'
      RETURN
      END
C
C     ================================================================
      SUBROUTINE WRITEOUT(NAME,CELL,X,ATOMSTORE,ATYPE,ITYPE,NPOINT,
     +                    PNTDEN,MAXTYPE,MAXATOM,NREADIN,IDIFF,
     +                    HOH,HOHALL,RESTYP)
C     ================================================================
C
C Writes out a psuedo-pdb file which contains the accessible area for
C each atom (if DIFFMODE OFF) or the area difference (any other DIFFMODE).
C
C This is an attempt to restore part of the function of the old areaimol
C program, although note that the data written by this new version is not
C as complete. Uses routines in RWBROOK and writes the areas (or
C differences) to the B-factor column...
C
C Arguments:
C
C INPUT:
C
C   NAME (character)               logical name of file to write to, i.e.
C                                  'XYZOUT'
C   CELL (real array)              cell parameters taken from input file
C   X    (real array)              x,y,z positions of atoms
C   ATOMSTORE (character*10 array) chunk of ATOM record containing seq.id.
C   ATYPE (character*2 array)      list of recognised atomic names
C   ITYPE (integer array)          no. of type of each atom (matched to ATYPE)
C   NPOINT (integer array)         no. of surface points of each atom
C   PNTDEN (real)                  no if points per sq. angstrom
C   MAXTYPE (integer)              maximum number of atom types
C   MAXATOM (integer)              maximum number of atoms
C   NREADIN (integer)              actual number of atoms stored
C   IDIFF (integer)                sets DIFFMODE
C   HOH (logical)                  set if in MODE HOH
C   HOHALL (logical)               set if in MODE HOHALL
C   RESTYP (logical*1 array)       set TRUE if atom is not water
C
      IMPLICIT NONE
C
C     ..Scalar arguments
      INTEGER MAXTYPE,MAXATOM,NREADIN,IDIFF
      REAL PNTDEN
      CHARACTER NAME*(*)
      LOGICAL HOH,HOHALL
C     ..
C     ..Array arguments
      INTEGER NPOINT(MAXATOM),ITYPE(MAXATOM)
      REAL CELL(6),X(3,MAXATOM)
      CHARACTER ATOMSTORE(MAXATOM)*10,ATYPE(0:MAXTYPE)*2
      LOGICAL*1 RESTYP(MAXATOM)
C     ..
C     ..Local Scalars
      INTEGER IFAIL,NOUT,I,IZ,IXYZOUT,IRESN
      REAL BISO,OCC
      CHARACTER LINE*80,RESNAM*3,CHNAM*1,RESNO*4,ATNAM*4
      CHARACTER INSCOD*1,ALTCOD*1,ID*4
C     ..
C     ..Local Arrays
      REAL BOUT(6)
C     ..
C     ..External Subroutines
      EXTERNAL XYZINIT,XYZOPEN,WREMARK,WBCELL,XYZATOM,XYZCOORD,
     +         XYZADVANCE,XYZCLOSE
C
C---- Initialise dummy variables
      NOUT = 0
      ALTCOD = ' '
      IZ     = 0
      BISO   = 0.0
      OCC    = 1.0
      BOUT(1)= 0.0
      BOUT(2)= 0.0
      BOUT(3)= 0.0
      BOUT(4)= 0.0
      BOUT(5)= 0.0
      BOUT(6)= 0.0
C
C---- Initialise for writing new file
C
      CALL XYZINIT
      IFAIL   = 0
      IXYZOUT = 0
      CALL XYZOPEN(NAME,'OUTPUT',' ',IXYZOUT,IFAIL)
C
C---- Write comments to header
C
      CALL WREMARK (IXYZOUT,
     + 'REMARK Output from AREAIMOL: this is not a standard PDB file')
C
      IF (IDIFF.EQ.0) THEN
        LINE = 
     + 'REMARK Accessible area for each atom is in B-factor column'
      ELSE
        LINE =
     + 'REMARK Area difference for each atom is in B-factor column'
      END IF
C
      CALL WREMARK(IXYZOUT,LINE)
C
      LINE = ' '
C
C---- Write the cell information
C
      CALL WBCELL(IXYZOUT,CELL,1)
C
C---- Write out the atoms with accessible area in "occ" column
C
      DO 200 I = 1,NREADIN
C
C---- Reject unwanted atoms:
C
C---- 1: If in HOH or HOHALL mode, do not write protein to output file
C
        IF (HOH.OR.HOHALL .AND. RESTYP(I)) GO TO 200
C
C---- 2: If in DIFFMODE WATER do not write protein to output file
C
        IF (IDIFF.EQ.2 .AND. RESTYP(I)) GO TO 200
C
C---- 3: If in DIFFMODE other than OFF then do not write zero areas out
C
        IF (IDIFF.NE.0 .AND. NPOINT(I).EQ.0) GO TO 200
C
C----   Increment counter
C
        NOUT = NOUT + 1
C
        LINE   = ATOMSTORE(I)
        ATNAM  = ATYPE(ITYPE(I))
        RESNAM = LINE(1:3)
        CHNAM  = LINE(5:5)
        RESNO  = LINE(6:9)
        READ(LINE(6:9),FMT='(I4)')IRESN
        INSCOD = LINE(10:10)
        ID     = ATNAM
        IF(ID(1:1).EQ.' ')ATNAM = ID(2:4)
C
        CALL XYZATOM(IXYZOUT,NOUT,ATNAM,RESNAM,CHNAM,IRESN,RESNO,
     +               INSCOD,ALTCOD,' ',IZ,ID)
C
        BOUT(1) = NPOINT(I)/PNTDEN
C
        CALL XYZCOORD(IXYZOUT,'O','U',X(1,I),X(2,I),X(3,I),
     +                OCC,BISO,BOUT)
C
        CALL XYZADVANCE(IXYZOUT,0,0,*200,*200)
C
C---- Next atom
C
 200  CONTINUE
C
C---- Close file
C
      CALL XYZCLOSE(IXYZOUT)
C
C---- Write some information to the log
C
      WRITE (6,FMT=6030) NOUT
 6030 FORMAT (/1X,I6,' ATOMS WRITTEN TO OUTPUT BROOKHAVEN FILE',/)
C
      RETURN
C
      END
C
C---- END OF WRITING OUTPUT FILE
C
C
C --- END OF AREAIMOL ---