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      ================
      PROGRAM CONTACTM
C      ================
C
C  March 1998 additional features
C
C  **********************************
C  Metal coordination geometry option
C  **********************************
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  METAL ZN
C  eof
C
C   or if there is no spacegroup name on the CRYST1 line
C   in the input PDB file use (or you can still select individual
C   symmetry operations with SYMMETRY keyword as in the original
C   contact.f )
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  METAL ZN
C  SPACEGROUP R3
C  eof
C
C This gives the output:
C   source atoms     target atoms       distance trans      symmetry operation
C  zn     1  ZN  His    10B CD2  ...  3.08    [      ]   1: -Y,  X-Y,  Z
C                His    10B CE1  ...  3.02    [      ]   9:   X, Y, Z
C                His    10B CE1  ...  3.02    [      ]   1: -Y,  X-Y,  Z
C                His    10B NE2  ...  2.11 ***[      ]   9:   X, Y, Z
C                His    10B NE2  ...  2.10 ***[      ]   1: -Y,  X-Y,  Z
C                His    10B NE2  ...  2.11 ***[      ]   2: -X+Y,  -X,  Z
C                Wat   201  O    ...  2.19 ***[      ]   1: -Y,  X-Y,  Z
C                Wat   201  O    ...  2.19 ***[      ]   2: -X+Y,  -X,  Z
C                Wat   251  O    ...  3.32    [      ]   9:   X, Y, Z  Target Spec Posn
C                Wat   251  O    ...  3.32    [      ]   1: -Y,  X-Y,  Z  Target Spec Posn
C                Wat   251  O    ...  3.32    [      ]   2: -X+Y,  -X,  Z  Target Spec Posn
C                His    10B CD2  ...  3.09    [      ]   9:   X, Y, Z
C                His    10B CD2  ...  3.09    [      ]   2: -X+Y,  -X,  Z
C                His    10B CE1  ...  3.02    [      ]   2: -X+Y,  -X,  Z
C                Wat   201  O    ...  2.19 ***[      ]   9:   X, Y, Z
C  zn     2  ZN  His    10D CD2  ...  3.14    [      ]   1: -Y,  X-Y,  Z
C                His    10D CD2  ...  3.14    [      ]   2: -X+Y,  -X,  Z
C                His    10D CE1  ...  2.96    [      ]   1: -Y,  X-Y,  Z
C                His    10D CE1  ...  2.96    [      ]   2: -X+Y,  -X,  Z
C                His    10D NE2  ...  2.08 ***[      ]   1: -Y,  X-Y,  Z
C                His    10D NE2  ...  2.08 ***[      ]   2: -X+Y,  -X,  Z
C                Wat   513  O    ...  2.32 ***[    -C]   9:   X, Y, Z
C                Wat   513  O    ...  2.32 ***[    -C]   2: -X+Y,  -X,  Z
C                His    10D CD2  ...  3.14    [      ]   9:   X, Y, Z
C                His    10D CE1  ...  2.96    [      ]   9:   X, Y, Z
C                His    10D NE2  ...  2.08 ***[      ]   9:   X, Y, Z
C                Wat   513  O    ...  2.32 ***[    -C]   1: -Y,  X-Y,  Z
C
C for Metal atom  ZN    ZN      1
C=============================================
C  1 His    10B NE2       
C  2 His    10B NE2           98.9
C  3 His    10B NE2           98.7    98.8
C  4 Wat   201  O             93.9    90.4   163.0
C  5 Wat   201  O            163.1    93.9    90.2    74.9
C  6 Wat   201  O             90.2   163.2    93.7    74.8    74.8
C                                1       2       3       4       5
C
C for Metal atom  ZN    ZN      2
C=============================================
C  1 His    10D NE2       
C  2 His    10D NE2          103.4
C  3 Wat   513  O            154.0    96.6
C  4 Wat   513  O             96.6    87.7    67.5
C  5 His    10D NE2          103.4   103.4    87.7   154.0
C  6 Wat   513  O             87.7   154.0    67.5    67.5    96.6
C                                1       2       3       4       5
C
C EXAMPLE WHERE METAL AND LIGANDS ARE IN SAME RESIDUE AND
C  NO SYMMETRY BONDS FOUND
C
C./Ncontact xyzin /homes/henrick/pdb/pdb1kdn.ent <<eof
CMETAL AL 2.25
Ceof
C
C   source atoms     target atoms       distance trans      symmetry operation
C Af3   157A AL  Adp   158A O1B  ...  3.42    [      ]   6:   X, Y, Z
C                Adp   158A PB   ...  3.39    [      ]   6:   X, Y, Z
C                Adp   158A O3B  ...  2.41    [      ]   6:   X, Y, Z
C                 mg   156A MG   ...  3.54    [      ]   6:   X, Y, Z
C                His   122A CB   ...  3.54    [      ]   6:   X, Y, Z
C                His   122A CG   ...  3.36    [      ]   6:   X, Y, Z
C                His   122A ND1  ...  2.43    [      ]   6:   X, Y, Z
C                His   122A CE1  ...  3.39    [      ]   6:   X, Y, Z
C                Af3   157A F2   ...  1.86 ***[      ]   6:   X, Y, Z
C                Af3   157A F3   ...  1.94 ***[      ]   6:   X, Y, Z
C                Af3   157A F1   ...  1.82 ***[      ]   6:   X, Y, Z
C Af3   157B AL  Adp   158B PB   ...  3.36    [      ]   6:   X, Y, Z
C                Adp   158B O3B  ...  2.45    [      ]   6:   X, Y, Z
C                His   122B CE1  ...  3.50    [      ]   6:   X, Y, Z
C                 mg   156B MG   ...  3.54    [      ]   6:   X, Y, Z
C                Af3   157B F3   ...  1.87 ***[      ]   6:   X, Y, Z
C                Adp   158B O1B  ...  3.29    [      ]   6:   X, Y, Z
C                His   122B CG   ...  3.51    [      ]   6:   X, Y, Z
C                His   122B ND1  ...  2.58    [      ]   6:   X, Y, Z
C                Af3   157B F1   ...  1.81 ***[      ]   6:   X, Y, Z
C                Af3   157B F2   ...  1.89 ***[      ]   6:   X, Y, Z
C Af3   157C AL  His   122C CB   ...  3.56    [      ]   6:   X, Y, Z
C                His   122C CG   ...  3.37    [      ]   6:   X, Y, Z
C                His   122C ND1  ...  2.48    [      ]   6:   X, Y, Z
C                His   122C CE1  ...  3.39    [      ]   6:   X, Y, Z
C                Af3   157C F2   ...  1.81 ***[      ]   6:   X, Y, Z
C                Af3   157C F3   ...  1.87 ***[      ]   6:   X, Y, Z
C                Adp   158C O3B  ...  2.39    [      ]   6:   X, Y, Z
C                Af3   157C F1   ...  1.91 ***[      ]   6:   X, Y, Z
C
C for Metal atom AF3  A AL    157
C=============================================
C  1 Af3   157A F2        
C  2 Af3   157A F3           108.6
C  3 Af3   157A F1           121.9   129.1
C for Metal atom AF3  B AL    157
C=============================================
C  1 Af3   157B F3        
C  2 Af3   157B F1           116.6
C  3 Af3   157B F2           117.0   126.4
C for Metal atom AF3  C AL    157
C=============================================
C  1 Af3   157C F2        
C  2 Af3   157C F3           115.9
C  3 Af3   157C F1           125.4   118.5
C
C *******************************
C Additional output for MODE AUTO
C *******************************
C
C
C  (as in the original contact.f the default is do SOURCE = all input
C   atoms to TARGET = all input atoms)
C   using
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  MODE AUTO
C  ATYPE ALL
C  LIMITS   2 3.66
C  eof
C
C  gives for LIMITS   2 6.0  Angstrom for PDB file 4ins
C
C  Sorted summation for Number of Contacts for
C  atoms between symmetry related molecules
C  Excluding water molecules
C
C Num contacts TransSymm       Symmetry
C=================================================
C     886        555 002      2: -X+Y,  -X,  Z
C     574        555 001      1: -Y,  X-Y,  Z
C     364        554 007      7: -Y+1/3,  X-Y+2/3,  Z+2/3
C     252        455 005      5: -X+Y+2/3,  -X+1/3,  Z+1/3
C     129        556 009      9:   X, Y, Z
C     118        556 001      1: -Y,  X-Y,  Z
C      90        554 002      2: -X+Y,  -X,  Z
C      63        554 009      9:   X, Y, Z
C       9        554 008      8: -X+Y+1/3,  -X+2/3,  Z+2/3
C       9        555 004      4: -Y+2/3,  X-Y+1/3,  Z+1/3
C       4        555 008      8: -X+Y+1/3,  -X+2/3,  Z+2/3
C       4        554 004      4: -Y+2/3,  X-Y+1/3,  Z+1/3
C       3        555 007      7: -Y+1/3,  X-Y+2/3,  Z+2/3
C
C Note: The 2 Zinc insulin dodecamer requires the symmetry operations
C  (in addition to identity)  -X+Y,  -X,  Z and  -Y,  X-Y,  Z, hence
C  in the above table of contacts the cutoff at 6Ang for a significant
C  oligomer is clear.
C
C  For LIMITS 2 3.66 Angstrom one gets
C    Sorted summation for Number of Contacts for
C    atoms between symmetry related molecules
C     Excluding water molecules
C
C Num contacts TransSymm       Symmetry
C=================================================
C      48        555 002      2: -X+Y,  -X,  Z
C      31        555 001      1: -Y,  X-Y,  Z
C      10        554 007      7: -Y+1/3,  X-Y+2/3,  Z+2/3
C       7        455 005      5: -X+Y+2/3,  -X+1/3,  Z+1/3
C       7        556 009      9:   X, Y, Z
C       7        554 009      9:   X, Y, Z
C       5        556 001      1: -Y,  X-Y,  Z
C       4        554 002      2: -X+Y,  -X,  Z
C
C Again the cutoff is clear.
C
C Note: The summation of all atom contacts between symmetry related
C  molecules excludes water molecules, but includes any other HETATM
C  labelled atoms.
C
C  Again if there is no spacegroup name on the CRYST1 record of the input
C  PDB file then the keyword SPACEGROUP can be used.
C
C
C **********************************
C Additional Atom Selection Keywords
C **********************************
C
C  Ncontact.f uses the same code as distang.f to allow
C  for atom selection. The original SOURCE and TARGET is still
C  functional. The new options are keywords FROM and TO.
C  e.g.
C
C  allowed inputs are
C
C  (i) based on atom numbers
C
C  FROM ATOM 1 TO 561
C  TO ATOM 1 TO 561
C
C  (ii) based on residue numbers but now allowing for chain names
C
C  FROM RESIDUE ALL CHAIN A 1 to 125
C  TO   RESIDUE ALL CHAIN W 1 to 256
C
C  The general input expression is:
C
C  FROM [ [ ATOM <inat0> [TO] <inat1> ] | 
C          [ RESIDUE <ires0> [TO] <ires1> ]] ...
C             [ CHAIN ALL | ONS | CA ]
C TO [ [ ATOM <jnat0> [TO] <jnat1> ] | 
C          [ RESIDUE <jres0> [TO] <jres1> ]] ... 
C             [ CHAIN ALL | ONS | CA ]
C  
C If ATOM is specified it is followed by <inat0> and <inat1>,
C respectively the first and last target atoms checked. I.e. FROM atoms
C <inat0> to <inat1> are checked against TO atoms <jnat0> to <jnat1>.
C (Atoms are numbered 1 to NAT, in the order read, but the residue order
C can be varied without restriction. Beware: atoms with occ=0.0 are not
C counted.)
C   
C If RESIDUE is specified it is followed by <jres0> and <jres1>,
C respectively the first and last target residues checked, and
C optionally subsidiary keywords:
C   
C CHAIN ALL | ONS | CA
C      ALL (default) will select all atoms in the requested residues.
C      ONs will select just the oxygen and nitrogen atoms. 
C      CA will select just the CA atoms.
C
C i.e. FROM residues <ires0> to <ires1> (CHAIN <type>) are
C     checked against TO residues <jres0> to <jres1> for the
C     appropriate class of atoms.
C          
C
C Example:
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  LIMITS   2 3.66
C  FROM RESIDUE ALL CHAIN A 1 to 21
C  TO   RESIDUE ALL CHAIN B 1 to 30
C  eof
C
C  will give just source CHAIN A to TARGET CHAIN B
C
C The ANGLE operation should be the same (except perhaps as I am
C  not certain about the original all waters with occupancy less that 1.0
C  are excluded from water analysis in this version)
C
C  In addition both residue names for waters of HOH or WAT are allowed.
C
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  LIMITS   2 3.66
C  ANGLE
C  eof
C
C  *****************************************
C  For MODE AUTO with TO and FROM keywords
C  *****************************************
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  MODE AUTO
C  LIMITS   2 3.66
C  FROM RESIDUE ALL CHAIN A 1 to 21
C  TO   RESIDUE ALL CHAIN D 1 to 30
C  eof
C
C will give all symmetry related chain D atoms to identity chain A atoms, i.e.
C
C Thr     8A O   His     5D CD2  ...  3.00         +C1: -Y,  X-Y,  Z
C Ser     9A C   His     5D NE2  ...  3.45         +C1: -Y,  X-Y,  Z
C Ser     9A O   His     5D CD2  ...  3.39         +C1: -Y,  X-Y,  Z
C                His     5D CE1  ...  3.62         +C1: -Y,  X-Y,  Z
C                His     5D NE2  ...  2.60 ***     +C1: -Y,  X-Y,  Z
C Gln    15A CG  Phe    25D CE2  ...  3.57     -A    5: -X+Y+2/3,-X+1/3,Z+1/3
C Asn    18A CG  Thr    27D OG1  ...  3.64     -A    5: -X+Y+2/3,-X+1/3,Z+1/3
C Asn    18A OD1 Thr    27D OG1  ...  2.75 *** -A    5: -X+Y+2/3,-X+1/3,Z+1/3
C
C
C  **************************************************
C  Combining selection TO/FROM with AUTO and ANGLE
C  **************************************************
C
C  ./Ncontact xyzin   /homes/henrick/pdb/pdb4ins.ent <<eof
C  MODE AUTO
C  LIMITS   2 3.66
C  FROM RESIDUE ALL CHAIN A 1 to 21
C  TO   RESIDUE ALL CHAIN D 1 to 30
C  ANGLE
C  eof
C
C
C
C  29 NOV 1988  The keyworded version. (t.s.)
C
C  25 NOV 1988  Automatic version which in AUTO mode
C               generates additional (primitive) lattice translation
C               and combines them with space-group symmetry supplied
C               on input cards in search for intermolecular contacts.
C               Since only single translations are used (e.g. +A ,-A,
C               +A-B ...etc) it works well only if the molecule is
C               reasonably positioned within the cell (not outside).
C               Limitation of DLIM restored again (it is necessary because
C               of the box size used). Some other changes i.e. increased
C               size of arrays for water contacts calculations. (T.S.)
C
C  8 SEPT 1988      Correction to IMOL mode so that contacts between
C            main-chain to main-chain or side-chain to side-chain
C            atoms on the same or adjacent Residues
C            are ONLY suppressed if the target symmetry
C            operation is the identifier. Otherwise true intermolecular
C            contacts (eg across a crystallographic 2-fold axis) will
C            be suppressed.
C            Remove limitation of DLIM to a maximum of 6.0A     A.G.W.L.
C
C *  This program computes various types of contacts for protein
C    structures. It is very fast and reads a standard Brookhaven data
C    bank file.
C
C  If MODE = IMOL or AUTO then the coordinate file must contain SCALE cards.
C
C    There are five modes of operation:
C
C    'ALL '  - for intramolecular distances (all atomic contacts for
C              chosen Residues)
C    'IRES'  - (default) for interResidual contacts (distances between
C              main-chain atoms from adjacent Residues suppressed)
C    'ISUB'  - for intersubunit contacts
C    'IMOL'  - for intermolecular contacts (this mode requires symmetry
C              cards). ** main-chain to main-chain and side-chain to
C              side-chain contacts are suppressed if the atoms are on
C              the same or adjacent Residues AND the target symmetry
C              operation is the identity.
C    'AUTO'  - as IMOL, but the search for intermolecular contacts
C              is done automatically. Only space group symmetry cards
C              should be supplied on input (excluding identity if you want
C              to find intermolecular contacts only)
C
C  The computations can be made using all
C  atom types (ATYPE  ALL) (default)
C  or non-carbon atoms only (ATYPE NONCARBON).
C
C  Command ANGLE (ANGLE angh ango dmin maxnb bmax) may be used to
C  compute hydrogen-bond angles and water hydrogen-bond statistics.
C
C  Maximum Residue number is 6000. Maximum number of atoms is 48000.
C  You may use up to 48 symmetry operations in IMOL and AUTO mode.
C
C      PARAMETER (MXATOM=48000,MAXRES=9000,MXSYMM=96)
C  MXSYMM is: no. of symmetry operations + 125
C
C  When using the ANGLE option, the program calculates the hydrogen
C  position for those target nitrogen atoms where the hydrogen position
C  is unambigous (ie excluding NZ on Lys and N terminus). The
C  angle O...H...N is calculated and printed. For source...oxygen
C  hydrogen bonds, the angle source...O__Bonded carbon is
C  calculated. Limits on both these angles (ANGH and ANGO) must be
C  supplied, and bonds with angles less than these limits are
C  rejected.
C  Suitable values are 120 and 90 degrees.
C  The ANGLE option can be used to search for hydrogen bonds within
C  the protein, and the bond angles will be calculated as described
C  above. Note that mainchain-mainchain and sidechain-sidechain
C  contacts within one Residue or to an immediatley adjacent Residue
C  are suppressed in this mode.
C The minimum bondlength (DMIN) for bonds to be included and the maximum
C  number of bonds (MAXNB) are read from the same card, in free format,
C  as is the MAXIMUM temperature factor for a source atom to be
C  included in contact search and water analysis (BMAX).
C
C  (See examples below)
C
C    Tadeusz Skarzynski,Imperial College,London   (February 1986)
C
C  Modified to provide analysis of H bonds involving waters, including
C  calculation of hydrogen positions July 1988, A.G.W. Leslie
C
C  Examples of  command files (note use of default values):
C
C 1)Intermolecular contacts:
C
C       contact xyzin holo9ca.brk << eof
C       TITLE INTER-MOLECULAR CONTACTS
C       ATYPE ALL
C       MODE  IMOL
C       limits  0.0  3.9
C       symm    1+X,Y,Z
C       symm    X,Y,1+Z
C       symm    -X,1/2+Y,-Z
C       symtit  +A
C       symtit  +C
C       symtit  21
C       eof
C
C 2) Intermolecular contacts (AUTO mode)
C
C       contact xyzin holo1-c127.brk << eof
C       mode   auto
C       symm 1/2+X,-Y,-Z
C       symm -X,1/2+Y,-Z
C       symm 1/2-X,1/2-Y,Z
C       eof
C
C 3)Contacts between specified parts of the molecule:
C
C       contact xyzin holo9c4.brk << eof
C       TITLE  NAD-WATER CONTACTS
C       limits 0.0  3.9
C       SOURCE 336
C       TARGET 337 1000
C       eof
C
C 4)Bonds lenghts in NAD molecule (Residue number 336):
C
C       contact xyzin holo9c4.brk << eof
C       TITLE NAD BONDS
C       mode    ALL
C       limits  0.0  1.9
C       SOURCE  336
C       TARGET  336
C       eof
C
C 5)Intersubunit contacts .Subunits must have different chain names.
c
C       contact xyzin holo9c4.brk<< eof
C       title  INTERSUBUNIT CONTACTS
C       mode   ISUB
C       limits 0.0  3.9
C       source 180 210
C       eof
C
C 6) Analysis of water hydrogen bonding (all water contacts)
C      contact xyzin holo9c4.brk << eof
C      TITLE WATER CONTACTS
C       angle
C      SOURCE  340 1000            !water molecules
C      eof
C
C 7) Hydrogen bonds and angles for a piece of chain
C
C       contact xyzin holo9c4.brk << eof
C       ANGLE
C       LIMITS   0.0 3.4
C       SOURCE   33 88
C       TARGET  1 334
C       eof
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalars in Common ..
      REAL AngleHY,AngleOx,Bmax,Dclose,DistMetalLigands,DistMin,Dlimit,
     +     FRdist,Xmax,Xmin,Ymax,Ymin,Zmax,Zmin
      INTEGER Icount,Icounts,Idents,maxNB,MetalCounter,NumSymm
      LOGICAL AutoFlag,LookForMetal,BigSearch
      CHARACTER ElementMetal*2,ModeAtomFlag*2,ModeFlag*2,Title*80
C     ..
C     .. Arrays in Common ..
      REAL Batom,FRclose,RFin,Rsymm,S,TR,XYZ,XYZS
      INTEGER NumMetalAngles,IsitAtom,IsitResidue,IsitChains
      INTEGER*2 InumAtmSymm,IresidNums,Nsymop,NUM,NUS,TranslationCode
      INTEGER*2 FromAtmSelected,ToAtomSelected
      LOGICAL SourceResFlag,TargetResFlag,DoChains,Hexclude
      CHARACTER ChainNames*1,AtomNames*4,ResidNames*4,
     +          TranslationTitles*9,Metals*80,SymmTitles*80
      CHARACTER ChainSelected*1
      INTEGER KtoAtmBegin,KtoAtmEnd,KFromAtmBegin,KFromAtmEnd,
     +        KtoResidSelect,KFromResidSelect,KSelectAtmType
C     ..
C     .. Local Scalars ..
      REAL Biso,Occupancy,TooCloseLimit,VOL,X,Y,Z
      INTEGER IcurrResIDnum,Ifail,Iserial,IXYZIN,Lcode,MaxBox,NumBigB,
     +        NX,NY,NZ,ZatomicNumber,Mto,Mfrom,Ldo,IRW,Ktofrom,Jdo
      CHARACTER AltCode*1,CurrChainName*1,InsertCode*1,CurrAtomName*4,
     +          CurrResIDname*4,CurrResno*4,pdbID*4,Rwork*4,segmentID*4
C     ..
C     .. Local Arrays ..
      REAL Bfactors(6),cell(6),Recipcell(6),ROin(4,4)
      INTEGER*2 NET(MaxSize)
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR,CCPFYP,CCPRCS,CONTACT,FRCLIM,GETRecipcell,
     +         KEYWORDS,RBcell,RBRORF2,SETNET,SYMGEN,SYMMAT,XYZADVANCE,
     +         XYZATOM,XYZCLOSE,XYZCOORD,XYZINIT,XYZOPEN
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC INT
C     ..
C     .. Common blocks ..
      COMMON /ATMCHR/ChainSelected(MaxSelectedChains,2)
      COMMON /SelectControl/IsitAtom,IsitResidue,IsitChains,DoChains,
     +       KtoAtmBegin(MaxSelectedAtoms),
     +       KtoAtmEnd(MaxSelectedAtoms),
     +       KFromAtmBegin(MaxSelectedAtoms),
     +       KFromAtmEnd(MaxSelectedAtoms),
     +       KtoResidSelect(MaxSelectedResidues),
     +       KFromResidSelect(MaxSelectedResidues),
     +       KSelectAtmType(MaxSelectedResidues)
      COMMON /AutoCom/AutoFlag,TranslationCode(2*MaxAtoms),
     +       TR(4,4,MaxSymmetry),BigSearch
      COMMON /BoundsCom/Xmin,Xmax,Ymin,Ymax,Zmin,Zmax
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharMetal/ElementMetal,Metals(MaxMetals)
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /CharTT/TranslationTitles(125)
      COMMON /DatMetals/MetalCounter,DistMetalLigands,LookForMetal,
     +       NumMetalAngles(MaxMetals)
      COMMON /InputNames/AtomNames(MaxAtoms),ResidNames(MaxAtoms),
     +       ChainNames(MaxAtoms)
      COMMON /Params/NUM(MaxAtoms),XYZ(3,MaxAtoms),Icount,
     +       Batom(MaxAtoms)
      COMMON /ResNumCom/IresidNums(MaxAtoms),
     +                  FromAtmSelected(MaxAtoms),
     +                  ToAtomSelected(MaxAtoms)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /WatCom/AngleHY,AngleOx,DistMin,maxNB,Bmax
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
       CHARACTER HYD(52)*4
C     ..
C     .. Data statements ..
       DATA HYD/'1HA','1HB','1HD','1HD1','1HD2','1HE','1HG','1HG1',
     +          '1HG2','1HZ','2HA','2HB','2HD','2HD1','2HD2','2HE',
     +          '2HG','2HG1','2HG2','2HZ','3HD1','3HD2','3HE','3HG1',
     +          '3HG2','H','HA','HB','HD1','HD2','HE1','HE2','HE3',
     +          'HG','HH2','HZ','HZ2','HZ3',
     +'1DH1','2DH1','1DH2','2DH2',' D  ','1DZ ','2DZ ','3DZ ',
     +'1DD2','2DD2',' DH ',' DG1','1DE2','2DE2'/
      DATA TranslationTitles/
     +      '-2A-2B-2C','-2A-2B-1C','-2A-2B   ','-2A-2B+1C','-2A-2B+2C',
     +      '-2A-1B-2C','-2A-1B-1C','-2A-1B   ','-2A-1B+1C','-2A-1B+2C',
     +      '-2A   -2C','-2A   -1C','-2A      ','-2A   +1C','-2A   +2C',
     +      '-2A+1B-2C','-2A+1B-1C','-2A+1B   ','-2A+1B+1C','-2A+1B+2C',
     +      '-2A+2B-2C','-2A+2B-1C','-2A+2B   ','-2A+2B+1C','-2A+2B+2C',
c
     +      '-1A-2B-2C','-1A-2B-1C','-1A-2B   ','-1A-2B+1C','-1A-2B+2C',
     +      '-1A-1B-2C','-1A-1B-1C','-1A-1B   ','-1A-1B+1C','-1A-1B+2C',
     +      '-1A   -2C','-1A   -1C','-1A      ','-1A   +1C','-1A   +2C',
     +      '-1A+1B-2C','-1A+1B-1C','-1A+1B   ','-1A+1B+1C','-1A+1B+2C',
     +      '-1A+2B-2C','-1A+2B-1C','-1A+2B   ','-1A+2B+1C','-1A+2B+2C',
c
     +      '   -2B-2C','   -2B-1C','   -2B   ','   -2B+1C','   -2B+2C',
     +      '   -1B-2C','   -1B-1C','   -1B   ','   -1B+1C','   -1B+2C',
     +      '      -2C','      -1C','         ','      +1C','      +2C',
     +      '   +1B-2C','   +1B-1C','   +1B   ','   +1B+1C','   +1B+2C',
     +      '   +2B-2C','   +2B-1C','   +2B   ','   +2B+1C','   +2B+2C',
c
     +      '+1A-2B-2C','+1A-2B-1C','+1A-2B   ','+1A-2B+1C','+1A-2B+2C',
     +      '+1A-1B-2C','+1A-1B-1C','+1A-1B   ','+1A-1B+1C','+1A-1B+2C',
     +      '+1A   -2C','+1A   -1C','+1A      ','+1A   +1C','+1A   +2C',
     +      '+1A+1B-2C','+1A+1B-1C','+1A+1B   ','+1A+1B+1C','+1A+1B+2C',
     +      '+1A+2B-2C','+1A+2B-1C','+1A+2B   ','+1A+2B+1C','+1A+2B+2C',
c
     +      '+2A-2B-2C','+2A-2B-1C','+2A-2B   ','+2A-2B+1C','+2A-2B+2C',
     +      '+2A-1B-2C','+2A-1B-1C','+2A-1B   ','+2A-1B+1C','+2A-1B+2C',
     +      '+2A   -2C','+2A   -1C','+2A      ','+2A   +1C','+2A   +2C',
     +      '+2A+1B-2C','+2A+1B-1C','+2A+1B   ','+2A+1B+1C','+2A+1B+2C',
     +      '+2A+2B-2C','+2A+2B-1C','+2A+2B   ','+2A+2B+1C','+2A+2B+2C'/
C     ..
C
C---- DEFAULT acceptable Metal---N  or Metal----O distance to be a bond
C
      DistMetalLigands = 2.35
      LookForMetal = .false.
      MetalCounter = 0
C
C          *********************************************
      CALL CCPFYP
      CALL CCPRCS(6,'CONTACT','$Date: 2002/07/31 13:22:53 $')
      CALL XYZINIT
C          **********************************************
C
      write (6,fmt=6000)
 6000 FORMAT (' ************** PROGRAM CONTACT *********')
C
C
      CALL KEYWORDS
C
C---- too close is angstrom
C
      TooCloseLimit = 0.5
C
C
      Xmin = 1000.0
      Xmax = -1000.0
      Ymin = 1000.0
      Ymax = -1000.0
      Zmin = 1000.0
      Zmax = -1000.0
      Icount = 0
      Ifail = 0
      IXYZIN = 0
      Ktofrom = 0
      Mfrom = 0
      Mto = 0
C
C          **************************************
      CALL XYZOPEN('XYZIN','INPUT',' ',IXYZIN,Ifail)
C          **************************************
C
      NumBigB = 0
      ROin(1,1) = -0.001
      RFin(1,1) = -0.001
C
C     Only get cell info in MODEs IMOL or AUTO
      IF (ModeFlag.EQ.'AU' .OR. ModeFlag.EQ.'IM') THEN
        CALL RBcell(cell,VOL)
        CALL RBRORF2(ROin,RFin,Lcode)
        CALL GETRecipcell(cell,Recipcell)
        CALL FRCLIM(Recipcell,TooCloseLimit,FRclose,FRdist)
      END IF
C
   10 CONTINUE
      CALL XYZADVANCE(IXYZIN,0,0,*20,*20)
      CALL XYZATOM(IXYZIN,Iserial,CurrAtomName,CurrResIDname,
     +             CurrChainName,IcurrResIDnum,CurrResno,InsertCode,
     +             AltCode,segmentID,ZatomicNumber,pdbID)
c
c--- track altcode for same residue
c    if so take 1st 
c
      CALL XYZCOORD(IXYZIN,'O','U',X,Y,Z,Occupancy,Biso,Bfactors)
C
C--- skip ANISOU lines
        IF (Bfactors(2) .ne. 0.0 .and. 
     +      Bfactors(3) .ne. 0.0) THEN
          IF (X .eq. 0.0 .and. 
     +        Y .eq. 0.0 .and. 
     +        Z .eq. 0.0) GO TO 10
        END IF
C
C--- exclude hydrogens ?
C
      IF (Hexclude) THEN
        DO 5555 Jdo=1,52
          IF (CurrAtomName .eq. HYD(Jdo)) GO TO 10
 5555   CONTINUE
      END IF
C
C--- is it Metal
C
      IF (LookForMetal) THEN
C
C--- will only work for connections to metal within
C    same residue if metal is first atom in list
C
C--- this will work
cHETATM 1150 AL   AF3 A 157      69.404  29.227   2.379  1.00 22.45
cHETATM 1151  F1  AF3 A 157      69.611  30.838   3.207  1.00 23.32
cHETATM 1152  F2  AF3 A 157      68.410  27.856   3.138  1.00 23.97
cHETATM 1153  F3  AF3 A 157      69.936  28.691   0.589  1.00 19.15
C
C--- this will fail
C
cHETATM 1151  F1  AF3 A 157      69.611  30.838   3.207  1.00 23.32
cHETATM 1152  F2  AF3 A 157      68.410  27.856   3.138  1.00 23.97
cHETATM 1153  F3  AF3 A 157      69.936  28.691   0.589  1.00 19.15
cHETATM 1150 AL   AF3 A 157      69.404  29.227   2.379  1.00 22.45
C
C---- Check for anisotropic Bs, skip ANISOU
C
CPJX: I'm moved this here because I think you need to skip
C     all the atoms (not just non-metals, otherwise metals get
C     counted twice)
C
        IF (CurrAtomName(1:2) .eq. ElementMetal) THEN
          MetalCounter = MetalCounter + 1
          IF (MetalCounter .gt. MaxMetals) THEN
            write (*,fmt='(/,x,a,i4,/,x,a,//,x,a,/)')
     +      'Warning: Maximum number of metal atoms is MaxMetals = ',
     +      MaxMetals,
     +      'and contact has found more than this in the input file',
     +      'Increase value of parameter MaxMetals'
            CALL CCPERR(1,' Too many Metal atoms found')
          END IF
          write (Rwork,fmt=6002) IcurrResIDnum
 6002     FORMAT (i4)
          Metals(MetalCounter) = CurrResIDname//' '//CurrChainName//
     +                           ' '//CurrAtomName//' '//Rwork
          SourceResFlag(IcurrResIDnum) = .true.
        END IF
      END IF
C
C---- accepting criteria here
C
      IF (SourceResFlag(IcurrResIDnum) .or.
     +    TargetResFlag(IcurrResIDnum)) THEN

        IF (ModeAtomFlag .eq. 'HO' .and. 
     +      ( CurrResIDname .eq. 'WAT ' .or. 
     +        CurrResIDname .eq. 'HOH ' )  ) THEN
            IF (Occupancy .lt. 0.01) GO TO 10
            IF (Biso .gt. Bmax) NumBigB = NumBigB + 1
        ELSE
C
C----  store all atoms of specified Residue ranges
C
          IF (ModeAtomFlag .ne. 'AL' .and. 
     +        ModeAtomFlag .ne. 'HO' .and.
     +        CurrAtomName(1:1) .eq. 'C') GO TO 10
          GO TO 2000
          END IF
          END IF
C
C---- add the extra selection criteria
C
      IF (IsitAtom .ne. 0 .or. IsitResidue .ne. 0) THEN
      IF (Occupancy .lt. 0.01) GO TO 10
C
C---- Find if this atom is used and set iw(..) ir(..)
C
      Ktofrom = Ktofrom + 1
      IRW = 0
      Mfrom = 0
      Mto = 0
C
C--- testing selection on atom numbers
C
      IF (IsitAtom .ne. 0) THEN
        DO 140 Ldo = 1,IsitAtom
          IF ((Ktofrom .eq. KFromAtmEnd(Ldo) .or. 
     +         Ktofrom .eq. KFromAtmBegin(Ldo))  .or. 
     +        (Ktofrom .le. KFromAtmEnd(Ldo) .and. 
     +         Ktofrom .ge. KFromAtmBegin(Ldo))) 
     +           MFrom = 1
          IF ((Ktofrom .eq. KtoAtmEnd(Ldo) .or. 
     +         Ktofrom .eq. KtoAtmBegin(Ldo))  .or. 
     +        (Ktofrom .le. KtoAtmEnd(Ldo) .and. 
     +         Ktofrom .ge. KtoAtmBegin(Ldo))) 
     +           MTo = 1
C
C----  No point in checking further if atom is in
C
          IF (Mfrom .eq. 1  .and.  Mto .eq. 1) GO TO 160
  140   CONTINUE
      END IF
C
C
      IF (IsitResidue .ne. 0) THEN 
        DO 150 Jdo = 1,IsItResidue
          IRW = 0
C
C---- ne/and GOTO 150  <OR> eq/or GOTO 157   which works ????
C
          IF (IcurrResIDnum .ne. KFromResidSelect(Jdo)  .and.  
     +        IcurrResIDnum .ne. KtoResidSelect(Jdo) )  GO TO 150
C
C
          IF (DoChains) THEN
          DO 155 Ldo = 1,IsitChains
          IF (CurrChainName .eq. ChainSelected(Ldo,1)   .or.
     +        CurrChainName .eq. ChainSelected(Ldo,2) ) GO TO 157
 155      CONTINUE
          GO TO 150
          END IF
 157      CONTINUE
C
C---- KSelectAtmType 1- ca only;  2 - O or N;  3- all atoms
C      Choose all atoms
C
          IF (KSelectAtmType(Jdo) .eq. 3) IRW = 1
C
C----  OR Choose CAs
C
          IF (CurrAtomName .eq. 'CA  '  .and.  
     +        KSelectAtmType(Jdo) .eq. 1) IRW = 1
C
C----  OR Choose Os and Ns
C
          IF (KSelectAtmType(Jdo) .eq. 2) THEN
            IF (CurrAtomName(1:1) .eq. 'O'  .or.  
     +          CurrAtomName(1:1) .eq. 'N') IRW = 1
          END IF
C
C
          IF (KFromResidSelect(Jdo) .ne. 0) 
     +           Mfrom = IRW
          IF (KtoResidSelect(Jdo) .ne. 0) 
     +           Mto = IRW
         IF (Mfrom .eq. 1  .and. Mto .eq. 1) GO TO 160
  150   CONTINUE
      END IF
C
C
  160 CONTINUE
        IF (Mfrom .eq. 0  .and. Mto .eq. 0) GO TO 10
      END IF
C
C
          Icount = Icount + 1
          FromAtmSelected(Icount) = Mfrom
          ToAtomSelected(Icount) = Mto
          GO TO 2010
C
C
 2000     CONTINUE
          Icount = Icount + 1

 2010     CONTINUE
          IF (Icount .gt. MaxAtoms) CALL CCPERR(1,
     +        ' Max atoms exceeded program terminating')
          NUM(Icount) = Icount
          IF (X .lt. Xmin) Xmin = X
          IF (X .gt. Xmax) Xmax = X
          IF (Y .lt. Ymin) Ymin = Y
          IF (Y .gt. Ymax) Ymax = Y
          IF (Z .lt. Zmin) Zmin = Z
          IF (Z .gt. Zmax) Zmax = Z
C
C
          AtomNames(Icount) = CurrAtomName
          ResidNames(Icount) = CurrResIDname
          ChainNames(Icount) = CurrChainName
          IresidNums(Icount) = IcurrResIDnum
          Batom(Icount) = Biso
          XYZ(1,Icount) = X
          XYZ(2,Icount) = Y
          XYZ(3,Icount) = Z
      GO TO 10
C     *****
C     at this point, the following arrays are loaded by selected atoms
C     XYZ([dimension <=3],[atom number <MaxAtoms=48000]) -> coordinates
C     [AtomNames,ResidNames,ChainNames,IresidNums,Batom]([atom number <MaxAtoms=48000])
C     NUM([atom number <MaxAtoms=48000]) -> pointer to atom number initialized with
C       atom number
C     *****
C     bricking:
C---- boxes 6 x 6 x 6 Angstr.
C

   20 CONTINUE
      IF (Xmin .lt. 0.0) THEN
          NX = NINT((Xmax-Xmin+6.1)/6.0) + 3
      ELSE
          NX = INT((Xmax-Xmin)/6.0) + 3
      END IF
      IF (Ymin .lt. 0.0) THEN
          NY = NINT((Ymax-Ymin+6.1)/6.0) + 3
      ELSE
          NY = INT((Ymax-Ymin)/6.0) + 3
      END IF
      IF (Zmin .lt. 0.0) THEN
          NZ = NINT((Zmax-Zmin+6.1)/6.0) + 3
      ELSE
          NZ = INT((Zmax-Zmin)/6.0) + 3
      END IF
      IF (NumBigB .ne. 0) write (6,fmt=6004) NumBigB,Bmax
 6004 FORMAT (/' ',I5,' Atoms (WATERS) Rejected because B greater than',
     +       F6.1)
      MaxBox = MaxSize/51 - 1
C
C
      IF (NX*NY*NZ .gt. MaxBox) THEN
        CALL CCPERR(1,'Too many boxes - change parameter MaxSize')
      ELSE
C
C----   symmetry related atoms generated here
C
        IF (ModeFlag .eq. 'IM') THEN
          CALL SYMMAT
          CALL SYMGEN
        END IF
C
C----   the array of boxes created here
C
C     NET(1,X,Y,Z) -> number of atoms in brick (X,Y,Z)
C     NET(I,X,Y,Z) -> pointer to general atom number NUM([index]) for atom I
C       in brick (X,Y,Z)
        IF (ModeFlag .ne. 'IM') 
     +     CALL SETNET(Xmin,Ymin,Zmin,NX,NY,NZ,NUM,
     +                             XYZ,Icount,NET)
C     NUS([atom number <MaxAtoms=48000]): numbers of symmetry related atoms 
        IF (ModeFlag .eq. 'IM') 
     +     CALL SETNET(Xmin,Ymin,Zmin,NX,NY,NZ,NUS,
     +                             XYZS,Icounts,NET)
C
C
      write (6,fmt=6800) Xmin,Xmax,Ymin,Ymax,Zmin,Zmax,
     +  NX,NY,NZ
 6800 FORMAT (/,' Xmin Xmax Ymin Ymax Zmin Zmax :',6F8.2,/,
     +          ' NX NY NZ boxes  ',3i8)
C
C---- computing contacts
C
C            ************************************
        CALL CONTACT(Xmin,Ymin,Zmin,NX,NY,NZ,NET)
C            ************************************
C
      END IF
      CALL XYZCLOSE(IXYZIN)
      CALL CCPERR(0,'Normal termination')
C
C
      END
C
C
C
C     ====================
      SUBROUTINE KEYWORDS
C     ====================
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains,MaxParams,Mkeys
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30,
     +           MaxParams=200,Mkeys=19)
C     ..
C     .. Scalars in Common ..
      REAL AngleHY,AngleOx,Bmax,Dclose,DistMetalLigands,DistMin,Dlimit,
     +     FRdist
      INTEGER Icounts,Idents,maxNB,MetalCounter,NumSymm
      LOGICAL AutoFlag,LookForMetal,BigSearch
      CHARACTER ElementMetal*2,ModeAtomFlag*2,ModeFlag*2,Title*80
C     ..
C     .. Arrays in Common ..
      REAL FRclose,RFin,Rsymm,S,TR,XYZS
      INTEGER NumMetalAngles,IsitAtom,IsitResidue,IsitChains,
     +        Kend,IRSC,J1,J2
      INTEGER*2 InumAtmSymm,Nsymop,NUS,TranslationCode
      LOGICAL SourceResFlag,TargetResFlag,DoChains
      CHARACTER TranslationTitles*9,Metals*80,SymmTitles*80
C     ..
C     .. Local Scalars ..
      REAL VOL
      INTEGER I,Ifail,IFL,Iflag,IFLT,Ikey,IXYZIN,J,Jdo,KE,L,LL,LLx,Lspg,
     +        LST,Mbeg,Mend,Mresidue,Mrest,Msource,Mtarget,Mwork,NSPG,
     +        NSYM,Nsymp,Ntok,NumSpg,KsymmTit,Nures,INAT0,INAT1,
     +        IRES0,IRES1,IMODE

      LOGICAL AngleFlag,GotSymmetry,Lend,NeedSymmetry,Hexclude,
     +        OldStyle
      CHARACTER TModeFlag*2,KEY*4,Cwork*6,PGname*10,SpgName*11,
     +          SpgWork*11,line*80,Sline*80,SpgNameCif*15,errlin*80
      CHARACTER CHNAM*1
      CHARACTER ChainSelected*1
C     ..
C     .. Local Arrays ..
      INTEGER JJ(6)
      REAL cell(6),Fvalue(MaxParams),T2(4,4,MaxSymmetry)
      INTEGER Ibeg(MaxParams),Idec(MaxParams),
     +        Iend(MaxParams),Ityp(MaxParams)
      INTEGER*2 SourceWork(MaxResidues),TargetWork(MaxResidues)
      CHARACTER Cvalue(MaxParams)*4,KEYS(Mkeys)*8,New2Old(3,2)*11,
     +          Rhomb(2,2)*11
      INTEGER KtoAtmBegin,KtoAtmEnd,KFromAtmBegin,KFromAtmEnd,
     +        KtoResidSelect,KFromResidSelect,KSelectAtmType
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR,CCPUPC,CHKKEY,GTPINT,GTPREA,MSYMLB3,
     +         PARSER,RBcell,RBSPGRP,RDSYMM,SIDENT,SYMFR3,SYMTRN2,
     +         TRsymm,XYZCLOSE,XYZOPEN
C     ..
C     .. Common blocks ..
      COMMON /ATMCHR/ChainSelected(MaxSelectedChains,2)
      COMMON /SelectControl/IsitAtom,IsitResidue,IsitChains,DoChains,
     +       KtoAtmBegin(MaxSelectedAtoms),
     +       KtoAtmEnd(MaxSelectedAtoms),
     +       KFromAtmBegin(MaxSelectedAtoms),
     +       KFromAtmEnd(MaxSelectedAtoms),
     +       KtoResidSelect(MaxSelectedResidues),
     +       KFromResidSelect(MaxSelectedResidues),
     +       KSelectAtmType(MaxSelectedResidues)
      COMMON /AutoCom/AutoFlag,TranslationCode(2*MaxAtoms),
     +       TR(4,4,MaxSymmetry),BigSearch
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharMetal/ElementMetal,Metals(MaxMetals)
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /CharTT/TranslationTitles(125)
      COMMON /DatMetals/MetalCounter,DistMetalLigands,LookForMetal,
     +       NumMetalAngles(MaxMetals)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /WatCom/AngleHY,AngleOx,DistMin,maxNB,Bmax
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C     .. External Functions ..
      INTEGER lenstr
      EXTERNAL lenstr
C     ..
C     .. Data statements ..
C
C
      DATA KEYS/'ATYPE','AMODE','MODE','ANGLE','LIMITS','SOURCE',
     +     'TARGET','SYMM','SYMTIT','TITLE','HOH','WAT','METAL',
     +     'SPACEGRO','FROM','TO','HEXCLUDE','BIGSEARC','END'/
      DATA New2Old(1,1),New2Old(1,2)/'P 21       ','P 1 21 1   '/
      DATA New2Old(2,1),New2Old(2,2)/'P 2        ','P 1 2 1    '/
      DATA New2Old(3,1),New2Old(3,2)/'C 2        ','C 1 2 1    '/
      DATA Rhomb(1,1),Rhomb(1,2)/'R 3        ','H 3        '/
      DATA Rhomb(2,1),Rhomb(2,2)/'R 32       ','H 32       '/
C     ..
C
C---- defaults
C
      AutoFlag = .false.
      AngleFlag = .false.
      DoChains = .false.
      NeedSymmetry = .false.
      GotSymmetry = .false.
      Hexclude = .false.
      BigSearch = .false.
      Idents = 0
      Dclose = 0.0
      Dlimit = 3.6
      Mresidue = 0
      Mrest = 0
      Msource = 0
      Mtarget = 0
      NumSymm = 0
      KsymmTit = 0
      Title = ' '
      ModeFlag = '**'
      ModeAtomFlag = 'AL'
      AngleHY = 120.0
      AngleOx = 90.0
      DistMin = 2.3
      maxNB = 4
      Bmax = 50
      IsitAtom = 0
      IsitResidue = 0
      IsitChains = 0
      KFromAtmBegin(1) = 0
      KtoAtmBegin(1) = 0
      KFromAtmEnd(1) = 0
      KtoAtmEnd(1) = 0
C
C
      DO 10 Jdo = 1,MaxResidues
        TargetResFlag(Jdo) = .true.
        SourceResFlag(Jdo) = .true.
   10 CONTINUE
      DO 15 Jdo = 1,MaxSelectedResidues
        KFromResidSelect(Jdo) = 0
        KtoResidSelect(Jdo) = 0
        KSelectAtmType(Jdo) = 0
 15   CONTINUE
C
C----  Loop to read control input
C
   20 KEY = ' '
      line = ' '
      Ntok = MaxParams
C
C          **********************************************************
      CALL PARSER(KEY,line,Ibeg,Iend,Ityp,Fvalue,Cvalue,Idec,Ntok,Lend,
     +            .true.)
C          **********************************************************
C
      IF (Lend) GO TO 270
      IF (KEY .eq. 'END') GO TO 270
C
C          ***************************
      CALL CHKKEY(KEY,KEYS,Mkeys,Ikey)
C          ***************************
C
C PJX: Reinstated the check below (missing from Kim's version)
      IF (IKEY.LE.0) THEN
        WRITE(ERRLIN,'(A,A4)') 
     +         'Error, unrecognised or ambiguous keyword: ',KEY
        CALL CCPERR(1,ERRLIN)
      ENDIF
C     IF (Ikey .le. 0) GO TO 20
      GO TO (30,30,40,50,80,90,130,170,250,
     +       260,50,50,60,200,600,700,1000,1010,270) Ikey
C
C---- Hexclude
C
 1000 Hexclude = .true.
      GO TO 20
C
C---- BIGSEARCH
C     =========
C
 1010 BigSearch = .true.
      GO TO 20
C
C---- ATYPE or ModeAtomFlag
C     ==============
C
   30 ModeAtomFlag = line(Ibeg(2) : (Ibeg(2)+1))
      CALL CCPUPC(ModeAtomFlag)
      GO TO 20
C
C---- ModeFlag
C     ====
C
   40 ModeFlag = line(Ibeg(2) : (Ibeg(2)+1))
      CALL CCPUPC(ModeFlag)
      GO TO 20
C
C---- ANGLE or HOH or WAT
C     ===================
C
   50 AngleFlag = .true.
C
C          *******************************
      CALL GTPREA(2,AngleHY,Ntok,Ityp,Fvalue)
      CALL GTPREA(3,AngleOx,Ntok,Ityp,Fvalue)
      CALL GTPREA(4,DistMin,Ntok,Ityp,Fvalue)
      CALL GTPINT(5,maxNB,Ntok,Ityp,Fvalue)
      CALL GTPREA(6,Bmax,Ntok,Ityp,Fvalue)
C          ********************************
C
      GO TO 20
C
C---- METAL
C     ===========
C
   60 CONTINUE
      LookForMetal = .true.
      ElementMetal = Cvalue(2) (1:2)
      CALL CCPUPC(ElementMetal)
      DistMetalLigands = 2.35
      IF (Ntok .ge. 3) 
     +     CALL GTPREA(3,DistMetalLigands,Ntok,Ityp,Fvalue)
      GO TO 20
C
C---- LIMITS
C     ======
C
C            *******************************
   80 CALL GTPREA(2,Dclose,Ntok,Ityp,Fvalue)
      CALL GTPREA(3,Dlimit,Ntok,Ityp,Fvalue)
C            *******************************
C
      GO TO 20
C
C---- SOURCE (source atoms)
C     =====================
C
   90 IF (Msource .eq. 0) THEN
C
        DO 100 I = 1,MaxResidues
          SourceResFlag(I) = .false.
  100   CONTINUE
C
        Msource = 1
      END IF
C
      IF (Ntok .eq. 2) THEN
C
C            ********************************
        CALL GTPINT(2,Mbeg,Ntok,Ityp,Fvalue)
C            ********************************
C
        Mend = Mbeg
      ELSE
C
C            ********************************
        CALL GTPINT(2,Mbeg,Ntok,Ityp,Fvalue)
        CALL GTPINT(3,Mend,Ntok,Ityp,Fvalue)
C            ********************************
C
      END IF
C
      IF (Mend .lt. Mbeg) THEN
        Mwork = Mend
        Mend = Mbeg
        Mbeg = Mwork
      END IF
C
      Nures = Mend - Mbeg + 1
      Mbeg = Mbeg - 1
C
      DO 120 I = 1,Nures
        Mresidue = Mresidue + 1
        IF (Mresidue .le. MaxResidues) GO TO 110
        write (line,fmt=6000) MaxResidues
 6000   FORMAT (
     +' *** Number of source Residues in one chain greater ',
     +         'than ',I5,' !!!')
        CALL CCPERR(1,line)
  110   SourceWork(Mresidue) = Mbeg + I
  120 CONTINUE
C
      GO TO 20
C
C---- TARGET (target atoms)
C     =====================
C
  130 IF (Mtarget .eq. 0) THEN
C
        DO 140 I = 1,MaxResidues
          TargetResFlag(I) = .false.
  140   CONTINUE
C
        Mtarget = 1
      END IF
C
      IF (Ntok .eq. 2) THEN
C
C            ********************************
        CALL GTPINT(2,Mbeg,Ntok,Ityp,Fvalue)
C            ********************************
C
        Mend = Mbeg
      ELSE
C
C            ********************************
        CALL GTPINT(2,Mbeg,Ntok,Ityp,Fvalue)
        CALL GTPINT(3,Mend,Ntok,Ityp,Fvalue)
C            ********************************
C
      END IF
      IF (Mend .lt. Mbeg) THEN
        Mwork = Mend
        Mend = Mbeg
        Mbeg = Mwork
      END IF
      Nures = Mend - Mbeg + 1
      Mbeg = Mbeg - 1
      DO 160 I = 1,Nures
        Mrest = Mrest + 1
        IF (Mrest .le. MaxResidues) GO TO 150
        write (line,fmt=6002) MaxResidues
 6002   FORMAT (
     +' *** Number of target Residues in one chain greater ',
     +         'than ',I5,' !!!')
        CALL CCPERR(1,line)
  150   TargetWork(Mrest) = Mbeg + I
  160 CONTINUE
      GO TO 20
C
C---- SYMM
C     =====
C
  170 CONTINUE
C
      Sline = ' '
      OldStyle = .true.
C
C---- Convert to upper case
C
C          ************
      CALL CCPUPC(line)
C          ************
C
      Sline = line(Ibeg(2) :)
      NumSymm = NumSymm + 1
      CALL SYMFR3(Sline,1,NumSymm,T2,KE)
      IF (KE .ne. 0) CALL CCPERR(1,'Bad symmetry operation')
C
C     elements are expected transposed from symfr3 convention
C
      DO 190 L = 1,4
        DO 180 LL = 1,4
          TR(LL,L,NumSymm) = T2(L,LL,NumSymm)
  180   CONTINUE
  190 CONTINUE
C
      write (6,fmt=6004) ((TR(I,J,NumSymm),I=1,4),J=1,3)
 6004 FORMAT (1X,'Matrix and vector',/3 (1X,3F10.4,5X,F10.4,/))
C
C---- test if this is unit matrix, no translation
C     this is needed to suppress contacts between identical atoms
C     or main-chain to main-chain contacts in same or adjacent Residue
C     when using IMOL ModeFlag but including the identity operation.
C     A.G.W.L. 8/9/88
C
C          **************************
      CALL SIDENT(TR(1,1,NumSymm),Iflag)
C          **************************
C
      IF (Iflag .eq. 1) THEN
C
C----  first test that identity has not been put in already
C
        IF (Idents .ne. 0) CALL CCPERR(1,
     +       'Identity matrix has been entered already')
C
C
        Idents = NumSymm
        write (6,fmt=6006) NumSymm
 6006   FORMAT (' Identity is symm op number ',I6)
      END IF
      GotSymmetry = .true.
      GO TO 20
C
C---- SPACEGROUP
C     ==========
C
  200 CONTINUE
      NSYM = 0
      OldStyle = .false.
      CALL RDSYMM(2,line,Ibeg,Iend,Ityp,Fvalue,Ntok,SpgName,NumSpg,
     +            PGname,NSYM,Nsymp,Rsymm)
C
C .... mistake here for P1 where nsym=1
c

      IF (NSYM .gt. 1) then
      NSYM = NSYM - 1
C
C---- IDENTS deliberately skipped
C
      DO 230 Jdo = 2,NSYM + 1
        LLx = Jdo - 1
        DO 220 L = 1,4
          DO 210 LL = 1,4
            TR(LL,L,LLx) = Rsymm(L,LL,Jdo)
  210     CONTINUE
  220   CONTINUE
        CALL SIDENT(TR(1,1,LLx),Iflag)
C
C--- shouldnt ever get here with Iflag = 1
C
        IF (Iflag .eq. 1) THEN
          IF (Idents .ne. 0) CALL CCPERR(1,' identity error ')
          Idents = LLx
        END IF
  230 CONTINUE
      KsymmTit = NSYM
      CALL SYMTRN2(NSYM+1,Rsymm,SymmTitles)
      DO 240 Jdo = 2,NSYM + 1
        write (Cwork,fmt=6008) Jdo - 1
 6008   FORMAT (i4,': ')
        SymmTitles(Jdo-1) = Cwork//SymmTitles(Jdo) (1:
     +                      lenstr(SymmTitles(Jdo)))
  240 CONTINUE
      GotSymmetry = .true.
      NumSymm = NSYM
         else
        DO 221 L = 1,4
          DO 211 LL = 1,4
            TR(LL,L,1) = Rsymm(L,LL,1)
 211     CONTINUE
 221  CONTINUE
        CALL SIDENT(TR(1,1,1),Iflag)
        NSYM = 1
      KsymmTit = 1
      CALL SYMTRN2(1,Rsymm,SymmTitles)
        write (Cwork,fmt=6008) 1
        SymmTitles(Jdo-1) = Cwork//SymmTitles(1) (1:
     +                      lenstr(SymmTitles(1)))
      GotSymmetry = .true.
      NumSymm = -1
        end if
C
C
      GO TO 20
C
C---- SYMTIT
C     ======
C
  250 KsymmTit = KsymmTit + 1
      SymmTitles(KsymmTit) = line(Ibeg(2) :Ibeg(2)+14)
      write (Cwork,fmt=6008) KsymmTit
      SymmTitles(KsymmTit) = Cwork//SymmTitles(KsymmTit) 
     +                     (1:lenstr(SymmTitles(KsymmTit)))

      GO TO 20
C
C---- Title
C     =====
C
  260 Title = line(Ibeg(2) :80)
      GO TO 20
C
C
C  FROM ATOM 1 TO 1000
C  TO ATOM 1 to 200000
C
C  FROM RES ALL CHAIN A 1 to 125
C  TO   RES ALL CHAIN W 1 to 256
C
C  FROM ATOM 1 TO 561
C  TO ATOM 1 TO 561
C
C  FROM ATOM 1320 TO 1350
C  TO ATOM 1320 to 1350
C
C  FROM ATOM 1 TO 10000
C  TO ATOM 1 TO 10000
C
C  FROM [ [ ATOM <inat0> [TO] <inat1> ] | [ RESIDUE <ires0> [TO] <ires1> ]] ...
C  [ CHAIN ALL | ONS | CA ]
C  
C  TO [ [ ATOM <jnat0> [TO] <jnat1> ] | [ RESIDUE <jres0> [TO] <jres1> ]] ... [
C  CHAIN ALL | ONS | CA ]
C  
C   Sets atom limits. (Default all atoms to all.)
C   
C   If ATOM is specified it is followed by <inat0> and <inat1>,
C   respectively the first and last target atoms checked. I.e. FROM atoms
C   <inat0> to <inat1> are checked against TO atoms <jnat0> to <jnat1>.
C   (Atoms are numbered 1 to NAT, in the order read, but the residue order
C   can be varied without restriction. Beware: atoms with occ=0.0 are not
C   counted.)
C   
C   If RESIDUE is specified it is followed by <jres0> and <jres1>,
C   respectively the first and last target residues checked, and
C   optionally subsidiary keywords:
C   
C   CHAIN ALL | ONS | CA
C          ALL (default) will select all atoms in the requested residues.
C          ONs will select just the oxygen and nitrogen atoms. CA will
C          select just the CA atoms.
C          
C   
C          i.e. FROM residues <ires0> to <ires1> (CHAIN <type>) are
C          checked against TO residues <jres0> to <jres1> for the
C          appropriate class of atoms.
C          
C
C---- FROM - must be followed by word atom or residu
C
C
C---- start FROM line
C  FROM [ [ ATOM <inat0> [TO] <inat1> ] | 
C       [ RESIDUE <ires0> [TO] <ires1> ]] 
C       [ CHAIN ALL | ONS | CA ]
C
  600 CONTINUE
C
C---- RDATOMSELECT gives the correct syntax automatically
C
          IFAIL = -1
          IMODE = 0
          CALL RDATOMSELECT(2,INAT0,INAT1,IRES0,IRES1,CHNAM,
     +                        IMODE,Ntok,line,Ibeg,Iend,Ityp,Idec,
     +                        Fvalue,IFAIL)
          IF (IFAIL.NE.0)
     +      CALL CCPERR (1,' Error reading FROM atom/residue selection')
C
C         ATOM
C         ====
          IF (INAT0.GT.0) THEN
            IsitAtom = IsitAtom + 1
            KFromAtmBegin(IsitAtom) = INAT0
            KFromAtmEnd(IsitAtom)   = INAT1
C
C         RESIDUE
C         =======
          ELSE IF (IRES0.GT.0) THEN
C
C           ALL/ONS/CA
C           ==========
            IF(IMODE.EQ.1) THEN
              IRSC = 3
              ModeAtomFlag = 'AL'
            ELSE IF (IMODE.EQ.2) THEN
              IRSC = 2
              ModeAtomFlag = 'NO'
            ELSE IF (IMODE.EQ.3) THEN
              IRSC = 1
            END IF
C
C           Get residue range
C           =================
            J1 = IRES0
            J2 = IRES1
            DO 603 Jdo = J1,J2
              IsitResidue = IsitResidue + 1
              KFromResidSelect(IsitResidue) = Jdo
              KSelectAtmType(IsitResidue) = IRSC
  603       CONTINUE
C
C           CHAIN
C           =====
            IF (CHNAM.NE.' ') THEN
              DoChains = .true.
              IsitChains = IsitChains + 1
              ChainSelected(IsitChains,1) = CHNAM
              ChainSelected(IsitChains,2) = ' '
            END IF
            GO TO 20
          END IF
C
C----  TO [ [ ATOM <jnat0> [TO] <jnat1> ] | 
C         [ RESIDUE <jres0> [TO] <jres1> ]]  
C         [  CHAIN ALL | ONS | CA ]
C
  700 CONTINUE
C
C---- RDATOMSELECT gives the correct syntax automatically
C
          IFAIL = -1
          IMODE = 0
          CALL RDATOMSELECT(2,INAT0,INAT1,IRES0,IRES1,CHNAM,
     +                        IMODE,NTOK,line,Ibeg,Iend,Ityp,Idec,
     +                        Fvalue,IFAIL)
          IF (IFAIL.NE.0)
     +      CALL CCPERR (1,' Error reading TO atom/residue selection')
C
C         ATOM
C         ====
          IF (INAT0.GT.0) THEN
C
C-------- Is this a second "TO" keyword for no "FROM"?
C         Fill in the "FROM" info
C
            IF (KtoAtmBegin(IsitAtom) .gt. 0) THEN
              IsitAtom = IsitAtom + 1
              KFromAtmBegin(IsitAtom) = KFromAtmBegin(IsitAtom-1)
              KFromAtmEnd(IsitAtom)   = KFromAtmEnd(IsitAtom-1)
            END IF
            KtoAtmBegin(IsitAtom) = INAT0
            KtoAtmEnd(IsitAtom)   = INAT1
C
C         RESIDUE
C         =======
          ELSE IF (IRES0.GT.0) THEN
C
C           ALL/ONS/CA
C           ==========
            IF(IMODE.EQ.1) THEN
              IRSC = 3
              ModeAtomFlag = 'AL'
            ELSE IF (IMODE.EQ.2) THEN
              IRSC = 2
              ModeAtomFlag = 'NO'
            ELSE IF (IMODE.EQ.3) THEN
              IRSC = 1
            END IF
C
C           Get residue range
C           =================
C
            J1 = IRES0
            J2 = IRES1
            DO 710 Jdo = J1,J2
              IsitResidue = IsitResidue + 1
              KtoResidSelect(IsitResidue) = Jdo
              KSelectAtmType(IsitResidue) = IRSC
  710       CONTINUE
C
C           CHAIN
C           =====
            IF (CHNAM.NE.' ') THEN
              DoChains = .true.
              IsitChains = IsitChains + 1
              ChainSelected(IsitChains,2) = CHNAM
              ChainSelected(IsitChains,1) = ' '
            END IF
          END IF
          GO TO 20
C
  270 CONTINUE
C
C---- checking input
C
C SET up for METAL keyword...
C
      IF (LookForMetal) THEN
C
        write (6,fmt='(/,x,a,/)')
     +  'Setting up options for METAL keyword:'
        IF (ModeFlag.eq.'**') THEN
          ModeFlag = 'AU'
          write (6,fmt='(x,a,/)') '   MODE not set: defaulting to AUTO'
        END IF
        IF (Dclose.lt.0.25) THEN
          Dclose = 0.25
          write (6,fmt='(x,2a,f5.2,/)')
     +    '   Minimum distance for contacts is too small',
     +    ' - resetting to ', Dclose
        END IF
        IF (Msource .eq. 0) THEN
          DO 70 I = 1,MaxResidues
            SourceResFlag(I) = .false.
   70     CONTINUE
          Msource = 1
        END IF
        write (6,fmt='(x,a,f5.2,a,/)')
     +  '   Atoms less than ',DistMetalLigands,
     +  ' A apart will be marked with a *** symbol' 
      END IF
C
C SET up the default Mode...
C
      IF (ModeFlag.eq.'**')THEN
        ModeFlag = 'IR'
        write (6,fmt='(/,x,a,/)') 'MODE not set: defaulting to IRES'
      END IF
C
      IF (IsitAtom .ne. 0  .or. IsitResidue .ne. 0) THEN
        DO 275 Jdo = 1,MaxResidues
          TargetResFlag(Jdo) = .false.
          SourceResFlag(Jdo) = .false.
 275    CONTINUE
      END IF
C
C
      IF (ModeFlag .eq. 'IM' .or. 
     +    ModeFlag .eq. 'AU') NeedSymmetry = .true.
C
C
      IF (ModeFlag .ne. 'AU') BigSearch = .false.
C
C---- For normal search need to reload the TranslationTitles array
      IF (.NOT. BigSearch) CALL RELOAD(TranslationTitles)
C
C if 1001
      IF (NeedSymmetry) THEN
        IF ( .not. GotSymmetry) THEN
C
C---- try getting from pdb file
C
          NSYM = 0
          SpgWork = ' '
          Ifail = 0
          IXYZIN = 0
          CALL XYZOPEN('XYZIN','INPUT',' ',IXYZIN,Ifail)
          CALL RBSPGRP(SpgWork)
          CALL RBcell(cell,VOL)
          CALL XYZCLOSE(IXYZIN)
          IF (SpgWork(1:1) .ne. ' ') THEN
C
C---- patch fudge etc for pdb file spacegroup names
C     to match symop.lib
C
C---- possible spacegroups??
C
cCRYST1   82.500   82.500   34.000  90.00  90.00 120.00 R 3          18
c A2      B2     C121     C1211   C2     C21    C222    C2221   C4212
c F222    F23    F422     F432    H3     H32    I121    I2      I212121
c I213    I222   I23      I4      I41    I4122  I4132   I422    I432
c P1      P1121  P1211    P2      P21    P21212 P212121 P21221  P213
c P22121  P2221  P3       P31     P3112  P3121  P32     P321    P3212
c P3221   P4     P41      P41212  P4122  P4132  P42     P4212   P422
c P42212  P4222  P43      P432    P43212 P4322  P4332   P6      P61
c P6122   P622   P6222    P63     P6322  P64    P6422   P65     P6522
c R3      R32
C
            DO 280 Jdo = 1,3
              IF (SpgWork .eq. New2Old(Jdo,2)) 
     +              SpgWork = New2Old(Jdo,1)
  280       CONTINUE
            DO 290 Jdo = 1,2
                IF (spgwork.eq.Rhomb(Jdo,1)) THEN
c
c---- check real compares
c
            IF (cell(4).lt.cell(5)+0.002     .and. 
     +          cell(4).gt.cell(5)-0.002     .and.
     +          cell(4).lt.   90.0 +0.002    .and.
     +          cell(4).gt.   90.0 -0.002    .and.
     +          cell(6).lt.  120.0 +0.002    .and.
     +          cell(6).gt.  120.0 -0.002) THEN
              spgwork = Rhomb(Jdo,2)
              GO TO 81
            END IF
          END IF
c
c
          IF (spgwork.eq.Rhomb(Jdo,2)) THEN
            IF (cell(1).lt.cell(2)+0.002 .and.
     +          cell(1).gt.cell(2)-0.002 .and.
     +          cell(1).lt.cell(3)+0.002 .and.
     +          cell(1).gt.cell(3)-0.002 ) THEN
              spgwork = Rhomb(Jdo,1)
             END IF
          GO TO 81
          END IF
  290       CONTINUE
C
C
 81         Lspg = lenstr(SpgWork)
            NSPG = 1
            SpgName = ' '
            DO 300 Jdo = 1,Lspg
              IF (SpgWork(Jdo:Jdo) .eq. ' ') GO TO 300
              SpgName(NSPG:NSPG) = SpgWork(Jdo:Jdo)
              NSPG = NSPG + 1
  300       CONTINUE
            LST = 22
            NumSpg = -99
            CALL MSYMLB3(LST,NumSpg,SpgName,SpgNameCif,PGname,symp,NSYM,
     1      Rsymm)
            IF (NumSpg .eq. -99) CALL CCPERR(1,' error in symmetry')
c
c
            write(6,8888) NSYM
 8888       format(' NSYM    ',i10)
c
c
            if (nsym.gt.1) then
            NSYM = NSYM - 1
C
C---- identity deliberately skipped
C
            DO 330 Jdo = 2,NSYM + 1
              LLx = Jdo - 1
              DO 320 L = 1,4
                DO 310 LL = 1,4
                  TR(LL,L,LLx) = Rsymm(L,LL,Jdo)
  310           CONTINUE
  320         CONTINUE
              Iflag = -99
              CALL SIDENT(TR(1,1,LLx),Iflag)
C
C--- shouldnt ever get here with Iflag = 1
C
              IF (Iflag .eq. 1) THEN
                IF (Idents .ne. 0) CALL CCPERR(1,' identity error ')
                Idents = LLx
              END IF
  330       CONTINUE

            KsymmTit = NSYM
            CALL SYMTRN2(NSYM+1,Rsymm,SymmTitles)
            DO 340 Jdo = 2,NSYM + 1
              write (Cwork,fmt=6008) Jdo - 1
              SymmTitles(Jdo-1) = Cwork//SymmTitles(Jdo) (1:
     +                            lenstr(SymmTitles(Jdo)))
  340       CONTINUE

            GotSymmetry = .true.
            NumSymm = NSYM
             else
            GotSymmetry = .true.
            NumSymm = -1
            KsymmTit = 1
              go to 350
            end if
c
c
            GO TO 350
C! spacegroup name from CRYST1 record
          END IF   
C
C---- need symmetry cant find it
C
      CALL CCPERR(1,'Contact options requested - require symm')
C
C
C  ! end not Gotsymmetry
        END IF     
C ! end if needsymmetry
      END IF       
C
C
  350 CONTINUE
C
C
      IF (AngleFlag) ModeAtomFlag = 'HO'
      write (6,fmt=6010) Title
 6010 FORMAT (1X,A)
C
C
      IF (ModeAtomFlag .eq. 'AL') THEN
        ModeAtomFlag = 'AL'
        write (6,fmt=6012)
 6012   FORMAT (' Atoms of all types will be used. ',/)
      END IF
C
C
      IF (ModeAtomFlag .eq. 'NO') THEN
        ModeAtomFlag = 'NO'
        write (6,fmt=6014)
 6014   FORMAT (/' Only non-carbon atoms will be used.',/)
      END IF
C
C
      IF (ModeAtomFlag .eq. 'HO') THEN
        ModeAtomFlag = 'HO'
        write (6,fmt=6016)
 6016   FORMAT (/,
     +' Hydrogen-bond angles will be computed, - ',
     +     ' only noncarbon contacts will be printed',/)
      END IF
C
      IF (ModeFlag .eq. 'IS') THEN
        TModeFlag = 'IS'
        write (6,fmt=6018)
 6018   FORMAT (/' Intersubunit contacts will be printed')
      END IF
C
      IF (ModeFlag .eq. 'AU') THEN
        write (6,fmt=6020)
 6020   FORMAT (/,
     +' *AMODE* ModeFlag - Automatic search for ',
     +     'intermolecular contacts')
        AutoFlag = .true.
        TModeFlag = 'IM'
        IF (BigSearch) write (6,fmt=6090)
 6090   FORMAT (/,
     +' *BigSearch* - automatic search will be extended to +/-2'
     +     ,/,'        lattice vectors in all directions')
      END IF
C
C
      IF (ModeFlag .eq. 'IM') THEN
        IF (NumSymm .eq. 0) CALL CCPERR(1,
     +      'You must supply SYMMETRY information')
        TModeFlag = 'IM'
        write (6,fmt=6022)
 6022   FORMAT (/,
     +' Intermolecular contacts will be printed',//,
     +'Contacts between main-chain atoms or between side-chain atoms',
     +/,' of the same or adjacent Residues suppressed if the target',
     +/,' symmetry operation is the identity')
      END IF
C
C
      IF (ModeFlag .eq. 'AL' .or. 
     +    ModeFlag .eq. 'IR') THEN
        TModeFlag = 'AL'
        write (6,fmt=6024)
 6024   FORMAT (
     +' Contacts between Residues from selected ranges will ',
     + 'be printed ',/,
     +' -  distances for main-chain atoms of ',
     +'adjacent Residues suppressed')
      END IF
C
C
      IF (ModeFlag .eq. 'IR') THEN
        TModeFlag = 'IR'
        write (6,fmt=6026)
 6026   FORMAT (' -  no intraResidual distances')
      END IF
C
C
      ModeFlag = TModeFlag
C
C
      IF (ModeAtomFlag .eq. 'HO') THEN
        write (6,fmt=6028) AngleHY,AngleOx,DistMin,maxNB,Bmax
 6028   FORMAT (/,
     +' Angle at hydrogen greater than    ',F6.1,/,
     +' Angle at oxygen greater than      ',F6.1,/,
     +' Minimum allowed separation        ',F5.2,' Angstroms',/,
     +' Maximum number of neighbours      ',I3,/,
     +' Only those atoms with B less than ',F6.1,' will be considered')
      END IF
C
C---- NumSymm is a flag
C
      IF (NumSymm .eq. -1) NumSymm = 0
C
C
      write (6,fmt=6030) Dclose,Dlimit
 6030 FORMAT (/,
     +' Atom-atom search will be done within the limits of ',
     +       F5.2,' to ',F5.2,' angstroms ')
C
C
      IF (ModeFlag .eq. 'IM') THEN
        IF (KsymmTit .ne. NumSymm .and. 
     +      KsymmTit .ne. 0) THEN
          write (6,fmt=6032)
 6032     FORMAT (/,
     +' *** WARNING - no. of SYMM and SYMTIT cards different')
          KsymmTit = 0
        END IF
C
C
        IF (NumSymm .gt. MaxSymmetry-125) THEN
          write (line,fmt=6034) NumSymm,MaxSymmetry - 125
 6034     FORMAT (' *** NumSymm = ',I3,' MAXIMUM ALLOWED IS ',I3)
          CALL CCPERR(1,line)
        END IF
C
C
        IF (KsymmTit .eq. 0) THEN
          KsymmTit = NumSymm
            IF (OLdStyle) THEN
ccx
c   SYMTRN2 expects matrix in different order to in array TR()
ccx
      DO 1360 Jdo = 1,KsymmTit
        DO 2360 L = 1,4
          DO 3360 LL = 1,4
            Rsymm(LL,L,Jdo) = TR(L,LL,Jdo)
 3360   CONTINUE
 2360   CONTINUE
 1360   CONTINUE
            CALL SYMTRN2(KsymmTit,Rsymm,SymmTitles)
            DO 360 Jdo = 1,KsymmTit
              write (Cwork,fmt=6008) Jdo
              SymmTitles(Jdo) = Cwork//SymmTitles(Jdo) 
     +                         (1:lenstr(SymmTitles(Jdo)))
  360       CONTINUE
          ELSE
          DO 365 Jdo = 1,NumSymm
            write (Cwork,fmt=6036) Jdo
 6036       FORMAT (I4,': ')
            SymmTitles(Jdo) = Cwork//'SYMMETRY '
 365     CONTINUE
          END IF
        END IF
C
C----  AutoFlag ModeFlag
C
        IF (AutoFlag) THEN
C
C--- if OldStyle should identity be added ---- NO 
C     as only required if requested - old style only added
C     extra TITLE 
C
          IF (OldStyle) THEN
            write (Cwork,fmt=6008) NumSymm + 1
            SymmTitles(NumSymm+1) = Cwork//'  X, Y, Z'
          ELSE IF (Idents .eq. 0) THEN
            Idents = NumSymm + 1
            write (Cwork,fmt=6008) NumSymm + 1
            SymmTitles(NumSymm+1) = Cwork//'  X, Y, Z'
          END IF
          if (NumSymm .eq.0) then
            Idents = 14
            if (BigSearch) Idents=63
          end if
          CALL TRsymm(NumSymm,TR,BigSearch)
        END IF
      END IF
C
C
      IF (Numsymm.le.1) then
          LLw= 1
        else
          llw=Numsymm + 1
        end if
      DO 370 Jdo = 1,LLw
        write (6,fmt=6038) SymmTitles(Jdo) (1:lenstr(SymmTitles(Jdo)))
        SymmTitles(Jdo) = SymmTitles(Jdo) (1:38)
        SymmTitles(Jdo) = SymmTitles(Jdo) (1:lenstr(SymmTitles(Jdo)))
 6038   FORMAT (a)
  370 CONTINUE
C
C
      IF ( .not. LookForMetal) THEN
        DO 380 I = 1,Mresidue
          IFL = SourceWork(I)
          SourceResFlag(IFL) = .true.
 380    CONTINUE
        DO 390 I = 1,Mrest
          IFLT = TargetWork(I)
          TargetResFlag(IFLT) = .true.
 390    CONTINUE
      ELSE
        DO 395 I = 1,MaxResidues
          TargetResFlag(I) = .true.
 395    CONTINUE
       END IF
C
C
      RETURN
 6040 FORMAT (/' Symmetry operation----  ',/' ',80A1,/)
C
      END
C
C
C
C     ===============================================
      SUBROUTINE CONTACT(Xmin,Ymin,Zmin,NX,NY,NZ,NET)
C     ===============================================
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains,MaxAngles
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30,
     +           MaxAngles=20)
C     ..
C     .. Scalar Arguments ..
      REAL Xmin,Ymin,Zmin
      INTEGER NX,NY,NZ
C     ..
C     .. Array Arguments ..
      INTEGER*2 NET(51,NX,NY,NZ)
C     ..
C     .. Scalars in Common ..
      REAL AngleHY,AngleOx,Bmax,Dclose,DistMetalLigands,DistMin,Dlimit,
     +     FRdist
      INTEGER Icount,Icounts,Idents,maxNB,MetalCounter,NumSymm
      LOGICAL AutoFlag,LookForMetal,BigSearch
      CHARACTER ElementMetal*2,ModeAtomFlag*2,ModeFlag*2,Title*80
C     ..
C     .. Arrays in Common ..
      REAL Batom,FRclose,RFin,Rsymm,S,TR,XYZ,XYZS
      INTEGER NumMetalAngles,IsitAtom,IsitResidue,IsitChains
      INTEGER*2 InumAtmSymm,IresidNums,Nsymop,NUM,NUS,TranslationCode
      INTEGER*2 FromAtmSelected,ToAtomSelected
      LOGICAL SourceResFlag,TargetResFlag,DoChains
      CHARACTER ChainNames*1,AtomNames*4,ResidNames*4,
     +          TranslationTitles*9,Metals*80,SymmTitles*80
C     ..
C     .. Local Scalars ..
      REAL ANGL,Anglim,Bwater,cutoff,D,Dclose2,Dlimit2,Dmark,X,Y,Z,ZANG
      INTEGER I,Idel,Iend,IPR,IR,Isymmetry,Itranslate,IX,IY,IZ,Jdo,K,
     +        MsourceAtom,KatominBrick,KAtransLoop,KBtransLoop,
     +        KCtransLoop,Kdo,KresidTarget,KsymmTarget,Ktarget,
     +        KtargetAtmNum,KtransTarget,Ldo,Mdo,MsourceResid,
     +        MtargetResid,Ncons,Ndo,Nhit,Nmet,Nnocon,NumAtmBrick
      LOGICAL Contct,IsSourceMainChainAtm,IsSourceSideChainAtm,
     +        IsTargetMainChainAtm,IsTargetSideChainAtm,Lclose,
     +        Lspecdone,Lspecial,Hexclude
      CHARACTER H1*1,H2*1,Ework*3,H3*3,H4*3,mark*3,RESI*4,Swork*4,
     +          Twork*4,ResSource*5,ResTarget*5,TargetSpecial*20
      CHARACTER ChainSelected*1
      INTEGER KtoAtmBegin,KtoAtmEnd,KFromAtmBegin,KFromAtmEnd,
     +        KtoResidSelect,KFromResidSelect,KSelectAtmType,NTRANS
C     ..
C     .. Local Arrays ..
      REAL delX(3),SaveAngles(3,MaxAngles,MaxMetals),temp1(3),temp2(3),
     +     temp3(3),WS(3),Zangles(MaxAngles),Zcoord(3,MaxMetals)
      INTEGER Inocon(MaxResidues),NAhits(MaxSymmetry),
     +        NBhits(MaxSymmetry),NChits(MaxSymmetry),
C<<     +        NsymmCounter(MaxSymmetry,27)
     +        NsymmCounter(MaxSymmetry,125)
C>>   modified HP990407 for 5 translations
      CHARACTER ortsym(125)*3,SaveNames(MaxAngles,MaxMetals)*25
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR,CCPLWC,TRANSFRM,GETANGL,GETH,GETO,getspecial,
     +         gettooclose,sorthits,WATER
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,INT,sqrt
C     ..
C     .. Common blocks ..
      COMMON /ATMCHR/ChainSelected(MaxSelectedChains,2)
      COMMON /AutoCom/AutoFlag,TranslationCode(2*MaxAtoms),
     +       TR(4,4,MaxSymmetry),BigSearch
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharMetal/ElementMetal,Metals(MaxMetals)
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /CharTT/TranslationTitles(125)
      COMMON /DatMetals/MetalCounter,DistMetalLigands,LookForMetal,
     +       NumMetalAngles(MaxMetals)
      COMMON /InputNames/AtomNames(MaxAtoms),ResidNames(MaxAtoms),
     +       ChainNames(MaxAtoms)
      COMMON /Params/NUM(MaxAtoms),XYZ(3,MaxAtoms),Icount,
     +       Batom(MaxAtoms)
      COMMON /ResNumCom/IresidNums(MaxAtoms),
     +                  FromAtmSelected(MaxAtoms),
     +                  ToAtomSelected(MaxAtoms)
      COMMON /SelectControl/IsitAtom,IsitResidue,IsitChains,DoChains,
     +       KtoAtmBegin(MaxSelectedAtoms),
     +       KtoAtmEnd(MaxSelectedAtoms),
     +       KFromAtmBegin(MaxSelectedAtoms),
     +       KFromAtmEnd(MaxSelectedAtoms),
     +       KtoResidSelect(MaxSelectedResidues),
     +       KFromResidSelect(MaxSelectedResidues),
     +       KSelectAtmType(MaxSelectedResidues)

      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /WatCom/AngleHY,AngleOx,DistMin,maxNB,Bmax
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C     .. External Functions ..
      INTEGER lenstr
      EXTERNAL lenstr
C     ..
C     .. Data statements ..
      DATA ortsym/
     +     '333','334','335','336','337',
     +     '343','344','345','346','347',
     +     '353','354','355','356','357',
     +     '363','364','365','366','367',
     +     '373','374','375','376','377',
c
     +     '433','434','435','436','437',
     +     '443','444','445','446','447',
     +     '453','454','455','456','457',
     +     '463','464','465','466','467',
     +     '473','474','475','476','477',
c
     +     '533','534','535','536','537',
     +     '543','544','545','546','547',
     +     '553','554','555','556','557',
     +     '563','564','565','566','567',
     +     '573','574','575','576','577',
c
     +     '633','634','635','636','637',
     +     '643','644','645','646','647',
     +     '653','654','655','656','657',
     +     '663','664','665','666','667',
     +     '673','674','675','676','677',
c
     +     '733','734','735','736','737',
     +     '743','744','745','746','747',
     +     '753','754','755','756','757',
     +     '763','764','765','766','767',
     +     '773','774','775','776','777'/
C     ..
C
C---- For normal search need to reload the ortsym array
C
      IF (.NOT. BigSearch) CALL RELOAD(ortsym)
C
C---- special cutoff is fractional
C
      cutoff = 0.02
          KMadd = 1
          IF (LookForMetal) KMadd = 3
C
C
C<<      DO 20 Jdo = 1,27
          IF (BigSearch) THEN
            NTRANS = 125
          ELSE
            NTRANS = 27
          END IF
C
          DO 20 Jdo = 1,NTRANS
C>> HP 990407 modified for 5 translations -2,-1,0,1,2
        DO 10 Kdo = 1,MaxSymmetry
          NsymmCounter(Kdo,Jdo) = 0
   10   CONTINUE
   20 CONTINUE
C
C
      Nmet = 0
      Nnocon = 0
      Iend = 0
      Dmark = 3.3
      IF (LookForMetal) Dmark = DistMetalLigands
      Dlimit2 = Dlimit*Dlimit
      Dclose2 = Dclose*Dclose
      mark = '   '
      write (6,fmt=6000)
 6000 FORMAT (/' LIST OF CONTACTS :',/' ==================',/)
      IF (ModeAtomFlag .eq. 'HO') write (6,fmt=6002)
 6002 FORMAT (
     +' If the target atom is a Nitrogen whose hydrogen position',/,
     +' is unambiguous, the angle given is that subtended at the ',/,
     +' calculated hydrogen position. If the target is an Oxygen, ',/,
     +' the angle is: source...oxygen___carbon atom bonded to the ',/,
     +' target oxygen.')
C
cpjx      IF (LookForMetal) THEN
cpjx        write (6,fmt=6001)
cpjx 6001   FORMAT (/,
cpjx     +'   source atoms     target atoms       distance trans',
cpjx     +        '      symmetry operation',/)
cpjx
cpjx      ELSE IF (ModeFlag .eq. 'IM') THEN
      IF (ModeFlag .eq. 'IM') THEN
        write (6,fmt=6004)
 6004   FORMAT (/,
     +'   source atoms     target atoms       distance angle',
     +        '      symmetry operation',/)
      ELSE
        write (6,fmt=6006)
 6006   FORMAT (/,
     +'   source atoms        target atoms       distance angle',/)
      END IF
C
C----  source atoms here
C      source atoms loop
C
      DO 160 MsourceAtom = 1,Icount
        IR = IresidNums(MsourceAtom)
C
C
       IF (DoChains ) THEN
          DO 700 Ldo = 1,IsitChains
            IF ( ChainNames(MsourceAtom) .eq. 
     +           ChainSelected(Ldo,1)) GO TO 710
 700      CONTINUE
          GO TO 160
 710      CONTINUE
       END IF
C
C  if 1000
        IF (SourceResFlag(IR) .or.
     +       FromAtmSelected(MsourceAtom) .ne. 0 ) THEN
C
C
          IF (LookForMetal) THEN
            IF (AtomNames(MsourceAtom) (1:2) .ne. 
     +          ElementMetal) GO TO 160
            Nmet = Nmet + 1
            NumMetalAngles(Nmet) = 0
            Zcoord(1,Nmet) = XYZ(1,MsourceAtom)
            Zcoord(2,Nmet) = XYZ(2,MsourceAtom)
            Zcoord(3,Nmet) = XYZ(3,MsourceAtom)
          END IF
C
C  if 1001
C----  suppress carbon source atoms in HOH ModeFlag
          IF (ModeAtomFlag .ne. 'HO' .or.
     +        AtomNames(MsourceAtom) (1:1) .ne. 'C' ) THEN
            H3 = AtomNames(MsourceAtom) (1:3)
            RESI = ResidNames(MsourceAtom)
            IF (RESI .eq. 'HOH ') RESI = 'WAT '
            Bwater = Batom(MsourceAtom)
            IsSourceMainChainAtm = (( H3 .eq. 'N  '.or.
     +                                H3 .eq. 'O  '.or.
     +                                H3 .eq. 'C  '.or.
     +                                H3 .eq. 'CA ') .and.
     +                              (RESI .ne. 'WAT '))
            IsSourceSideChainAtm = ( .not. IsSourceMainChainAtm .and.
     +                              (RESI .ne. 'WAT '))
            IPR = 1
            Contct = .false.
            X = XYZ(1,MsourceAtom) - Xmin
            Y = XYZ(2,MsourceAtom) - Ymin
            Z = XYZ(3,MsourceAtom) - Zmin
            IF (Xmin .lt. 0.0) THEN
               IX = NINT((X+6.1)/6.0) + 1
            ELSE
               IX = INT(X/6.0) + 2
            END IF
            IF (Ymin .lt. 0.0) THEN
               IY = NINT((Y+6.1)/6.0) + 1
            ELSE
               IY = INT(Y/6.0) + 2
            END IF
            IF (Zmin .lt. 0.0) THEN
               IZ = NINT((Z+6.1)/6.0) + 1
            ELSE
               IZ = INT(Z/6.0) + 2
            END IF
            write (ResSource,fmt='(I4,a)') IresidNums(MsourceAtom),
     +        ChainNames(MsourceAtom)
            Ework = ResidNames(MsourceAtom) (2:4)
            CALL CCPLWC(Ework)
            Swork = ResidNames(MsourceAtom) (1:1)//Ework
            IF (Swork .eq. 'Hoh') Swork = 'Wat '
C
C----  (27 bricks searched) now 125 bricks HP990407
C
            DO 150 KAtransLoop = IX - KMadd,IX + KMadd
               IF (KAtransLoop .le. 0 .or. 
     +                 KAtransLoop .gt. NX) GO TO 150
              DO 140 KBtransLoop = IY - KMadd,IY + KMadd
                 IF (KBtransLoop .le. 0 .or. 
     +                 KBtransLoop .gt. NY) GO TO 140
                DO 130 KCtransLoop = IZ - KMadd,IZ + KMadd
                   IF (KCtransLoop .le. 0 .or. 
     +                 KCtransLoop .gt. NZ) GO TO 130
C
C---- number of atoms in the brick
C
                  NumAtmBrick = NET(1,KAtransLoop,KBtransLoop,
     +                              KCtransLoop)
C
C  if  1002
                  IF (NumAtmBrick .ne. 0) THEN
C
C---- "IM" loop
C  if 1003
                    IF (ModeFlag .eq. 'IM') THEN
C
C----  loop for intermolecular contacts
C
C---- atoms-in-one-brick loop
C
                      DO 70 KatominBrick = 2,NumAtmBrick + 1
C
C----Target atom no.
c original  IA = NET(II,J,K,L)  KtargetAtmNum
c           IAA = NUMS(IA)      Ktarget
C
                        KtargetAtmNum = NET(KatominBrick,KAtransLoop,
     +                                   KBtransLoop,KCtransLoop)
                        KsymmTarget = Nsymop(KtargetAtmNum)
                        KtransTarget = TranslationCode(KtargetAtmNum)
                        Ktarget = InumAtmSymm(KtargetAtmNum)
                        KresidTarget = IresidNums(Ktarget)
C
C
       IF (DoChains ) THEN
          DO 800 Ldo = 1,IsitChains
            IF ( ChainNames(Ktarget) .eq. 
     +           ChainSelected(Ldo,2)) GO TO 810
 800      CONTINUE
          GO TO 70
 810      CONTINUE
       END IF
C
C  if 1004
       IF (TargetResFlag(KresidTarget)   .or.
     +      ToAtomSelected(Ktarget) .ne. 0) THEN
C
C  if 1005
           IF (ModeAtomFlag .ne. 'HO' .or.
     +         AtomNames(Ktarget) (1:1) .ne. 'C'  .or. 
     +      ToAtomSelected(Ktarget) .ne. 0) THEN
C
C----  suppress identical atoms contacts if target sym op is identity
C     is this .or. or is .and. ????
C  if 1006
        IF (AtomNames(MsourceAtom) .ne. AtomNames(Ktarget) .or.
     +      KsymmTarget .ne. Idents  ) THEN

          IF (.not. LookForMetal .and.
     +      (IresidNums(MsourceAtom) .eq. IresidNums(Ktarget)) )
     +      GO TO 333
C
C----  suppress main chain-main chain and side chain-side chain
C      contacts on same or adjacent Residues, but only if symmetry
C      operation for target is the identity.
C
          H4 = AtomNames(Ktarget) (1:3)
          RESI = ResidNames(KtargetAtmNum)
          IF (RESI .eq. 'HOH ') RESI = 'WAT '
          IsTargetMainChainAtm = ((H4 .eq. 'N  ' .or.
     +                             H4 .eq. 'O  ' .or.
     +                             H4 .eq. 'C  ' .or.
     +                             H4 .eq. 'CA ') .and. 
     +                             (RESI .ne. 'WAT '))
          IsTargetSideChainAtm = ( .not. IsTargetMainChainAtm .and.
     +                           (RESI .ne. 'WAT '))
C
C
         IF ((IsSourceMainChainAtm .and. IsTargetMainChainAtm .or.
     +         IsSourceSideChainAtm .and.
     +         IsTargetSideChainAtm) .and.
     +         KsymmTarget .eq. Idents) THEN
             Idel = ABS(IresidNums(MsourceAtom)-
     +                  IresidNums(Ktarget))
          IF (.not. LookForMetal .and. Idel .lt. 2) GO TO 60
         END IF
C
C
         DO 30 Jdo = 1,3
            WS(Jdo) = XYZS(Jdo,KtargetAtmNum)
            delX(Jdo) = ABS(XYZS(Jdo,KtargetAtmNum) -
     +                      XYZ(Jdo,MsourceAtom))
ccx            IF (delX(Jdo) .gt. Dlimit ) GO TO 60
   30    CONTINUE
C
C
         D = delX(1)*delX(1) + 
     +        delX(2)*delX(2) +
     +         delX(3)*delX(3)
         Lspecial = .false.
         Lspecdone = .false.
         CALL transfrm(WS,RFin)
ccx         CALL getspecial(WS,Lspecial,Lspecdone,cutoff)
ccx         IF (Lspecdone) GO TO 60
C
C
         IF (Lspecial) THEN
             TargetSpecial = '  Target Spec Posn'
         ELSE
             Lclose = .false.
             CALL gettooclose(WS,Lclose)
ccx             IF (Lclose) GO TO 60
             TargetSpecial = ' '
         END IF
C
C  if 1007
         IF (D .le. Dlimit2 .and. D .ge. Dclose2) THEN
           D = sqrt(D)
           H1 = AtomNames(MsourceAtom) (1:1)
           H2 = AtomNames(Ktarget) (1:1)
           IF (H1 .eq. 'O' .or. 
     +         H1 .eq. 'N' ) THEN
C
C---- this didnt mark Al---F bonds
C 
               IF (H2 .eq. 'O' .or. 
     +             H2 .eq. 'N') THEN
                   mark = '*  '
                   IF (D .lt. Dmark) mark = '***'
               END IF
           END IF
C
C--- so do metal Dmark here
C  so what bonds can a metal have ???? O,N,S,Cl,Br,I,F
C howabout NOT carbon -- no good for some case ?
C
           IF (LookForMetal .and. 
     +         H1 .eq. ElementMetal(1:1)   .and.
     +         D .lt. Dmark) mark = '***'
C
C
           ANGL = 0.0
           KsymmTarget = Nsymop(KtargetAtmNum)
           KtransTarget = TranslationCode(KtargetAtmNum)
           MsourceResid = IresidNums(Ktarget)
C
C
           IF (ResidNames(MsourceAtom) .ne. 'WAT ' .and. 
     +         ResidNames(MsourceAtom) .ne. 'HOH ' .and.
     +         ResidNames(Ktarget) .ne. 'WAT ' .and.
     +         ResidNames(Ktarget) .ne. 'HOH ') 
     +       NsymmCounter(KsymmTarget,KtransTarget) = 
     +           NsymmCounter(KsymmTarget,KtransTarget) + 1
C
C
          IF (ModeAtomFlag .eq. 'HO') THEN
            IF (ResidNames(Ktarget) .ne. 'WAT ' .and. 
     +          ResidNames(Ktarget) .ne. 'HOH ') THEN
C
C----  If target is an oxygen, calculate the source...target
C      O---bonded carbon angle. If it is a nitrogen and the
C      hydrogen position is fixed, calculate
C      the source...hydrogen---nitrogen angle.
C
         IF (AtomNames(Ktarget) (1:1) .eq. 'O') THEN
             CALL GETO(AtomNames(Ktarget),ResidNames(Ktarget),
     +                 MsourceResid,KtargetAtmNum,KsymmTarget,
     +                 XYZ(1,MsourceAtom),XYZS(1,KtargetAtmNum),
     +                 ANGL)
             Anglim = AngleOx
         ELSE IF (AtomNames(Ktarget) (1:1) .eq. 'N') THEN
              CALL GETH(AtomNames(Ktarget),ResidNames(Ktarget),
     +                  MsourceResid,KtargetAtmNum,KsymmTarget,
     +                  XYZ(1,MsourceAtom),XYZS(1,KtargetAtmNum),
     +                  ANGL)
               Anglim = AngleHY
         END IF
        END IF
       END IF
C
C
       KsymmTarget = Nsymop(KtargetAtmNum)
       KtransTarget = TranslationCode(KtargetAtmNum)
       MsourceResid = IresidNums(MsourceAtom)
       MtargetResid = IresidNums(Ktarget)
C
C----  compile statistics on water H bonds
C
      IF (ResidNames(MsourceAtom) .eq. 'WAT ' .and. 
     +    ModeAtomFlag .eq. 'HO') CALL WATER(MsourceResid,
     +    MtargetResid,AtomNames(Ktarget),ResidNames(Ktarget),D,
     +    ANGL,Iend,MsourceAtom,KtargetAtmNum,Bwater)
      IF (ResidNames(MsourceAtom) .eq. 'HOH ' .and. 
     +    ModeAtomFlag .eq. 'HO') CALL WATER(MsourceResid,
     +    MtargetResid,AtomNames(Ktarget),ResidNames(Ktarget),D,
     +    ANGL,Iend,MsourceAtom,KtargetAtmNum,Bwater)
C
C----  test angle
C  if 1008
      IF (ANGL .eq. 0.0 .or.
     +    ANGL .ge. Anglim) THEN
          Contct = .true.
          write (ResTarget,fmt='(I4,a)')
     +           IresidNums(Ktarget),ChainNames(Ktarget)
          Ework = ResidNames(Ktarget) (2:4)
          CALL CCPLWC(Ework)
          Twork = ResidNames(Ktarget) (1:1)// Ework
          IF (Twork .eq. 'Hoh') Twork = 'Wat '
C
C  if 1009
          IF (IPR .eq. 1) THEN
            IF (ModeAtomFlag .eq. 'HO') THEN
              write (6,fmt=6008) Swork,ResSource,
     +    AtomNames(MsourceAtom),Twork,ResTarget,
     +    AtomNames(Ktarget),D,mark,ANGL,
     +    TranslationTitles(KtransTarget),
     +    SymmTitles(KsymmTarget) (1:lenstr(SymmTitles(KsymmTarget))),
     +    TargetSpecial(1:lenstr(TargetSpecial))
          ELSE
              write (6,fmt=6010) Swork,ResSource,
     +    AtomNames(MsourceAtom),Twork,ResTarget,
     +    AtomNames(Ktarget),D,mark,
     +    TranslationTitles(KtransTarget),
     +    SymmTitles(KsymmTarget) (1:lenstr(SymmTitles(KsymmTarget))),
     +    TargetSpecial(1:lenstr(TargetSpecial))
 6008     FORMAT (1X,A,1X,A,1X,A,A,1X,A,1X,A,' ... ',f5.2,1X,A,F6.1,
     +            '[',A,']',A,A)
 6010     FORMAT (1X,A,1X,A,1X,A,A,1X,A,1X,A,' ... ',f5.2,1X,A,'[',A,
     +            ']',A,A)
C
C---- this made angles with Dmark
C
          IF (mark .eq. '***' .and. LookForMetal) THEN
              NumMetalAngles(Nmet) = NumMetalAngles(Nmet) + 1
              IF (NumMetalAngles(Nmet) .gt. MaxAngles) CALL CCPERR(1,
     +        ' too many angles')
              SaveNames(NumMetalAngles(Nmet),Nmet) = Twork//' '//
     +             ResTarget//' '// AtomNames(Ktarget)
              DO 40 Jdo = 1,3
                SaveAngles(Jdo,NumMetalAngles(Nmet),Nmet) = 
     +               XYZS(Jdo,KtargetAtmNum)
   40         CONTINUE
          END IF
        END IF
C
C
       IPR = 0
       ELSE IF (ModeAtomFlag .eq. 'HO') THEN
           write (6,fmt=6012) Twork,ResTarget,AtomNames(Ktarget),D,
     +     mark,ANGL,TranslationTitles(KtransTarget),
     +     SymmTitles(KsymmTarget) (1:lenstr(SymmTitles(KsymmTarget))),
     +     TargetSpecial(1:lenstr(TargetSpecial))
       ELSE
           write (6,fmt=6014) Twork,ResTarget,AtomNames(Ktarget),D,
     +     mark,TranslationTitles(KtransTarget),
     +     SymmTitles(KsymmTarget) (1:lenstr(SymmTitles(KsymmTarget))),
     +     TargetSpecial(1:lenstr(TargetSpecial))
 6012  FORMAT (16X,A,1X,A,1X,A,' ... ',f5.2,1X,1A,F6.1,'[',A,']',A,
     +         A)
 6014  FORMAT (16X,A,1X,A,1X,A,' ... ',f5.2,1X,1A,'[',A,']',A,A)
       IF (mark .eq. '***' .and.
     +     LookForMetal) THEN
           NumMetalAngles(Nmet) = NumMetalAngles(Nmet) + 1
           IF (NumMetalAngles(Nmet) .gt. MaxAngles) CALL CCPERR(1,
     +         ' too many angles')
           SaveNames(NumMetalAngles(Nmet),Nmet) = Twork//' '//
     +             ResTarget//' '//AtomNames(Ktarget)
           DO 50 Jdo = 1,3
             SaveAngles(Jdo,NumMetalAngles(Nmet),Nmet) = 
     +          XYZS(Jdo,KtargetAtmNum)
   50      CONTINUE
       END IF
C
C  end if 1009
                                   END IF
C
C end if 1008
                                END IF
C
C end if 1007
                              END IF
C
C end if 1006
                            END IF
 333                          CONTINUE
C
C end if 1005
                          END IF
C
C end if 1004
                        END IF
C
C
   60 mark = '   '
   70 CONTINUE
C
C else if 1003
                    ELSE
C
C  ModeFlag  .ne.  'IM'
C
C----  loop for internal contacts
C       atoms-in-one-brick loop
C
      DO 120 KatominBrick = 2,NumAtmBrick + 1
C
C---- target atom no.
C
       KtargetAtmNum = NET(KatominBrick,KAtransLoop,
     +                                   KBtransLoop,KCtransLoop)
       IF (DoChains ) THEN
cpjx        IF (ModeFlag .eq. 'IR') GO TO 910
          DO 900 Ldo = 1,IsitChains
            IF ( ChainNames(KtargetAtmNum) .eq. 
     +           ChainSelected(Ldo,2)) GO TO 910
 900      CONTINUE
          GO TO 120
 910      CONTINUE
       END IF
C
C  if 1101
       IF (KtargetAtmNum .ne. MsourceAtom) THEN
           KresidTarget = IresidNums(KtargetAtmNum)
C
C  if 1102
           IF (TargetResFlag(KresidTarget)  .or. 
     +      ToAtomSelected(KtargetAtmNum) .ne. 0) THEN
C
C----  suppress contacts involving carbons if ModeFlag is HOH
C  if 1103
      IF (ModeAtomFlag .ne. 'HO' .or.
     +   AtomNames(KtargetAtmNum) (1:1) .ne. 'C') THEN
C
C--- suppress 
C  if 1104 
C
ccx         IF (ModeFlag .ne. 'IR' ) THEN
         IF (ModeFlag .eq. 'AL' ) GO TO 5555

          IF (ModeFlag .eq. 'IR' .and. 
     +     IresidNums(MsourceAtom) .eq. IresidNums(KtargetAtmNum) )
     +       GO TO 111

       IF (ModeFlag .ne. 'IR' ) THEN
          IF (.not. LookForMetal .and.
     +     IresidNums(MsourceAtom) .eq. IresidNums(KtargetAtmNum) )
     +       GO TO 111

          IF (.not. LookForMetal .and.
     +    ChainNames(MsourceAtom) .eq. ChainNames(KtargetAtmNum) )
     +       GO TO 111
       END IF
 5555     CONTINUE
C
C  if 1105
C<<       IF (ModeFlag .ne. 'IS' ) THEN
          IF (ModeFlag .ne. 'IS' .or. chainNames(MsourceAtom) 
     +         .ne. ChainNames(KtargetAtmNum) ) THEN
C>> correction from CCP4, Kim
       IF (ModeFlag .ne. 'IR' ) THEN
           IF (ChainNames(MsourceAtom) .eq.
     +         ChainNames(KtargetAtmNum) .and. ModeFlag .ne. 'AL')
     +          GO TO 111
       END IF
C
C----  suppress contacts between main-chain atoms
C      on adjacent Residues
C
           IF (ModeFlag .ne. 'AL') THEN
              H4 = AtomNames(KtargetAtmNum) (1:3)
              RESI = ResidNames(KtargetAtmNum)
              IF (RESI .eq. 'HOH ') RESI = 'WAT '
              IsTargetMainChainAtm = ((H4 .eq. 'N  ' .or.
     +                                 H4 .eq. 'O  ' .or.
     +                                 H4 .eq. 'C  ' .or.
     +                                 H4 .eq. 'CA ') .and.
     +                                (RESI .ne. 'WAT '))
              IsTargetSideChainAtm = ( .not. 
     +             IsTargetMainChainAtm .and.
     +             (RESI .ne. 'WAT '))
              IF (IsSourceMainChainAtm .and.
     +            IsTargetMainChainAtm) THEN
                  Idel = ABS(IresidNums(MsourceAtom) - 
     +                    IresidNums(KtargetAtmNum))
          IF (.not. LookForMetal .and. Idel .lt. 2) GO TO 110
                                       END IF
           END IF
C
C
           DO 80 Jdo = 1,3
            delX(Jdo) = ABS(XYZ(Jdo,KtargetAtmNum) -
     +                      XYZ(Jdo,MsourceAtom))
            IF (delX(Jdo) .gt. Dlimit) GO TO 110
   80      CONTINUE
           D = delX(1)*delX(1) +
     +          delX(2)*delX(2) + 
     +           delX(3)*delX(3)
C
C  if 1106
           IF (D .le. Dlimit2 .and.
     +         D .ge. Dclose2) THEN
               D = sqrt(D)
               H1 = AtomNames(MsourceAtom) (1:1)
               H2 = AtomNames(KtargetAtmNum) (1:1)
               IF (H1 .eq. 'O' .or. 
     +             H1 .eq. 'N' .or.
     +             H1 .eq. ElementMetal(1:1)) THEN
                   IF (H2 .eq. 'O' .or. H2 .eq. 'N') THEN
                      mark = '*  '
                      IF (D .lt. Dmark) mark = '***'
                   END IF
                END IF
                ANGL = 0.0
                MsourceResid = IresidNums(KtargetAtmNum)
C
C  if 1107
                IF (ModeAtomFlag .eq. 'HO') THEN
C
C  if 1108
                 IF (ResidNames(KtargetAtmNum) .ne. 'WAT ' .and. 
     +               ResidNames(KtargetAtmNum) .ne. 'HOH ') THEN
C
C----  If target is an oxygen, calculate the source...target
C      O---bonded carbon angle. If it is a nitrogen and the
C      hydrogen position is fixed, calculate
C      the source...hydrogen---nitrogen angle.
C
      IF (AtomNames(KtargetAtmNum) (1:1) .eq. 'O') THEN
        CALL GETO(AtomNames(KtargetAtmNum),ResidNames(KtargetAtmNum),
     +            MsourceResid,KtargetAtmNum,KsymmTarget,
     +            XYZ(1,MsourceAtom),XYZ(1,KtargetAtmNum),ANGL)
        Anglim = AngleOx
      ELSE IF (AtomNames(KtargetAtmNum) (1:1) .eq. 'N') THEN
        CALL GETH(AtomNames(KtargetAtmNum),ResidNames(KtargetAtmNum),
     +            MsourceResid,KtargetAtmNum,KsymmTarget,
     +            XYZ(1,MsourceAtom),XYZ(1,KtargetAtmNum),ANGL)
        Anglim = AngleHY
      END IF
C
C end if 1108
        END IF
C
C  end if 1107
        END IF
        MsourceResid = IresidNums(MsourceAtom)
        MtargetResid = IresidNums(KtargetAtmNum)
C
C----  compile statistics on water H bonds
C
        IF (ResidNames(MsourceAtom) .eq. 'WAT ' .and. 
     +      ModeAtomFlag .eq. 'HO') CALL WATER(MsourceResid,
     +      MtargetResid,AtomNames(KtargetAtmNum),
     +      ResidNames(KtargetAtmNum),D,ANGL,Iend,
     +      MsourceAtom,KtargetAtmNum,Bwater)
        IF (ResidNames(MsourceAtom) .eq. 'HOH ' .and. 
     +      ModeAtomFlag .eq. 'HO') CALL WATER(MsourceResid,
     +      MtargetResid,AtomNames(KtargetAtmNum),
     +      ResidNames(KtargetAtmNum),D,ANGL,Iend,
     +      MsourceAtom,KtargetAtmNum,Bwater)
C
C----  test angle
C  if 1109
       IF (ANGL .eq. 0.0 .or. ANGL .ge. Anglim) THEN
           Contct = .true.
           write (ResTarget,fmt='(I4,a)') IresidNums(KtargetAtmNum),
     +     ChainNames(KtargetAtmNum)
           Ework = ResidNames(KtargetAtmNum)(2:4)
           CALL CCPLWC(Ework)
           Twork = ResidNames(KtargetAtmNum)(1:1)//Ework
           IF (Twork .eq. 'Hoh') Twork = 'Wat '
C
C  if 1110
           IF (IPR .eq. 1) THEN
C
C  if 1111
             IF (ModeAtomFlag .eq. 'HO') THEN
             write (6,fmt=6016) Swork,ResSource,AtomNames(MsourceAtom),
     +       Twork,ResTarget,AtomNames(KtargetAtmNum),D,mark,ANGL
             ELSE
             write (6,fmt=6018) Swork,ResSource,AtomNames(MsourceAtom),
     +       Twork,ResTarget,AtomNames(KtargetAtmNum),D,mark
 6016        FORMAT (1X,A,1X,A,1X,A,'  ...   ',A,1X,A,1X,A,
     +               '   ...  ',f5.2,1X,A,F6.1)
 6018        FORMAT (1X,A,1X,A,1X,A,'  ...   ',A,1X,A,1X,A,
     +               '   ...  ',f5.2,1X,A)
C
C if 1112
             IF (mark .eq. '***' .and. LookForMetal) THEN
                 NumMetalAngles(Nmet) = NumMetalAngles(Nmet) + 1
                 IF (NumMetalAngles(Nmet) .gt. MaxAngles)
     +              CALL CCPERR(1,' too many angles')
                 SaveNames(NumMetalAngles(Nmet),Nmet) = Twork//
     +             ' '//ResTarget//' '//AtomNames(KtargetAtmNum)
                 DO 90 Jdo = 1,3
                    SaveAngles(Jdo,NumMetalAngles(Nmet),Nmet) = 
     +                  XYZS(Jdo,KtargetAtmNum)
   90            CONTINUE
C
C end if 1112
            END IF
C
C end if 1111
         END IF
         IPR = 0
         ELSE IF (ModeAtomFlag .eq. 'HO') THEN
          write (6,fmt=6020) Twork,ResTarget,
     +                       AtomNames(KtargetAtmNum),
     +                       D,mark,ANGL
         ELSE
          write (6,fmt=6022) Twork,ResTarget,AtomNames(KtargetAtmNum),
     +                        D,mark
 6020     FORMAT (18X,'...   ',A,1X,A,1X,A,'   ...  ',f5.2,1X,A,F6.1)
 6022     FORMAT (18X,'...   ',A,1X,A,1X,A,'   ...  ',f5.2,1X,A)
          IF (mark .eq. '***' .and. LookForMetal) THEN
              NumMetalAngles(Nmet) = NumMetalAngles(Nmet) + 1
              IF (NumMetalAngles(Nmet) .gt. MaxAngles) 
     +            CALL CCPERR(1,' too many angles')
              SaveNames(NumMetalAngles(Nmet),Nmet) = Twork//' '//
     +                   ResTarget//' '// AtomNames(KtargetAtmNum)
              DO 100 Jdo = 1,3
                SaveAngles(Jdo,NumMetalAngles(Nmet),Nmet) = 
     +              XYZ(Jdo,KtargetAtmNum)
  100         CONTINUE
          END IF
C
C end if 1110
            END IF
C
C end if 1109
             END IF
C
C end if 1106
              END IF
C
C end if 1105
               END IF
 111           continue
C
C end if 1104
ccx                END IF
C
C end if 1103
                 END IF
C
C end if 1102
                  END IF
C
C end if 1101
                   END IF
C
C
  110    mark = '   '
  120    CONTINUE
C
C end if 1003
         END IF
C
C end if 1002
         END IF
  130    CONTINUE
  140    CONTINUE
  150    CONTINUE
C
C----  if ModeFlag is HOH, see if a contact has been found
C
       IF (ModeAtomFlag .eq. 'HO' .and.  .not. Contct) THEN
              Nnocon = Nnocon + 1
              Inocon(Nnocon) = IresidNums(MsourceAtom)
       END IF
C
C end if 1001
       END IF
C
C end if 1000
       END IF
C
C
  160 CONTINUE
C
C----   loop over numangles  for zinc source number 1 value
C
C PJX: This prints out a table of angles
C These are the angles at the metal position for each pair of
C contacting atoms listed
C
      IF (Nmet .ge. 1) THEN
        write (6,fmt='(//,2(x,a,/),/,2(x,a,/))')
     +  'LIST OF ANGLES for metal atoms',
     +  '==============================',
     +  'These are the angles at each metal position formed by',
     +  'the metal and each pair of close-contacting atoms listed'
        DO 200 Ndo = 1,Nmet
          IF (NumMetalAngles(Ndo) .gt. 1) THEN
            write (6,fmt=6024) Metals(Ndo) (1:lenstr(Metals(Ndo)))
 6024       FORMAT (/,' for Metal atom ',a,/,' ',45 ('='),/)
            K = 1
            write (6,fmt=6026) K,SaveNames(1,Ndo) (1:20)
 6026       FORMAT (' ',i2,' ',a20,' ',20f8.1)
            DO 190 Jdo = 2,NumMetalAngles(Ndo)
              DO 180 Kdo = 1,Jdo - 1
                DO 170 Ldo = 1,3
                  temp1(Ldo) = SaveAngles(Ldo,Jdo,Ndo)
                  temp2(Ldo) = SaveAngles(Ldo,Kdo,Ndo)
                  temp3(Ldo) = Zcoord(Ldo,Ndo)
  170           CONTINUE
                CALL GETANGL(temp1,temp3,temp2,ZANG)
                Zangles(Kdo) = ZANG
  180         CONTINUE
              write (6,fmt=6026) Jdo,SaveNames(Jdo,Ndo) (1:20),
     +          (Zangles(Mdo),Mdo=1,Jdo-1)
  190       CONTINUE
            write (6,fmt=6028) (K,K=1,NumMetalAngles(Ndo)-1)
 6028       FORMAT (25X,20i8)
          ELSE
            write (6,fmt='(/,x,a,a,/)')
     +      'No angles found for metal atom ',
     +      Metals(Ndo) (1:lenstr(Metals(Ndo)))
          END IF
  200   CONTINUE
      END IF
C
C----  final call to WATER to get statistics
C
      IF (ModeAtomFlag .eq. 'HO') THEN
        IF (Nnocon .ne. 0) write (6,fmt=6030) Nnocon,
     +      (Inocon(I),I=1,Nnocon)
 6030   FORMAT (/,
     +' ',I4,' Source ATOMS (WATERS) have NO Contacts',/,
     +' The Residue Numbers are',/ (1X,16I5))
        Iend = 1
        I = 1
        Ktarget = 1
        MsourceResid = 1
        MtargetResid = 1
C
C
        CALL WATER(MsourceResid,MtargetResid,AtomNames(Ktarget),
     +             ResidNames(Ktarget),D,ANGL,Iend,I,Ktarget,Bwater)
C
C
      END IF
C
C
      IF ( .not. LookForMetal ) THEN
      IF ( AutoFlag ) THEN
      write (6,fmt=6032)
 6032 FORMAT (
     +'  Sorted summation for Number of Contacts for',/,
     +'  atoms between symmetry related molecules',/,
     +'  Excluding water molecules',/)
      Nhit = 0
      DO 220 Jdo = 1,MaxSymmetry
C>>        DO 210 Kdo = 1,27
         DO 210 Kdo = 1,125
C<<   modified HP990407 for 5 translations
          IF (NsymmCounter(Jdo,Kdo) .gt. 0) THEN
            Isymmetry = Jdo
            Itranslate = Kdo
            Ncons = NsymmCounter(Jdo,Kdo)
            Nhit = Nhit + 1
            NAhits(Nhit) = Ncons
            NBhits(Nhit) = Isymmetry
            NChits(Nhit) = Itranslate
          END IF
  210   CONTINUE
  220 CONTINUE
C
C--- next sort to max number
C    find max
C
      CALL sorthits(Nhit,NAhits,NBhits,NChits)
      IF (Nhit .gt. 1) THEN
        write (6,fmt=6034)
 6034   FORMAT (
     +' Num contacts TransSymm       Symmetry',/,
     +'=================================================')
        DO 230 Jdo = Nhit,1,-1
          Isymmetry = NBhits(Jdo)
          Itranslate = NChits(Jdo)
          IF (Isymmetry .le. 9) THEN
            write (6,fmt=6036) NAhits(Jdo),ortsym(Itranslate),
     +        Isymmetry,SymmTitles(NBhits(Jdo)) (1:
     +        lenstr(SymmTitles(NBhits(Jdo))))
          ELSE IF (Isymmetry .gt. 9 .and. Isymmetry .le. 99) THEN
            write (6,fmt=6038) NAhits(Jdo),ortsym(Itranslate),
     +        Isymmetry,SymmTitles(NBhits(Jdo)) (1:
     +        lenstr(SymmTitles(NBhits(Jdo))))
          ELSE IF (Isymmetry .gt. 99) THEN
            write (6,fmt=6040) NAhits(Jdo),ortsym(Itranslate),
     +        Isymmetry,SymmTitles(NBhits(Jdo)) (1:
     +        lenstr(SymmTitles(NBhits(Jdo))))
          END IF
 6036     FORMAT (i8,8X,a3,' 00',i1,3X,a)
 6038     FORMAT (i8,8X,a3,' 0',i2,3X,a)
 6040     FORMAT (i8,8X,a3,' ',i3,3X,a)
  230   CONTINUE
      END IF
      END IF
C
C
      END IF
C
C
      END
C
C
C
C      ==============================
      SUBROUTINE GETAMH(X1,X2,X3,XH)
C      ==============================
C
C----  calculate coords of amide hydrogen. Nitrogen coords in X2.
C
C     .. Array Arguments ..
      REAL X1(3),X2(3),X3(3),XH(3)
C     ..
C     .. Local Arrays ..
      REAL A(3),B(3),C(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL UNIT,VSUM,VDIF
C     ..
C
C          **************
      CALL VDIF(A,X1,X2)
      CALL UNIT(A)
      CALL VDIF(B,X3,X2)
      CALL UNIT(B)
      CALL VSUM(C,A,B)
      CALL UNIT(C)
      CALL VDIF(XH,X2,C)
C          **************
C
      END
C
C
C
C      ================================
      SUBROUTINE GETANGL(X1,X2,X3,ANG)
C      ================================
C
C     .. Scalar Arguments ..
      REAL ANG
C     ..
C     .. Array Arguments ..
      REAL X1(3),X2(3),X3(3)
C     ..
C     .. Local Scalars ..
      REAL CANG,X
C     ..
C     .. Local Arrays ..
      REAL A(3),B(3)
C     ..
C     .. External Functions ..
      REAL DOT
      EXTERNAL DOT
C     ..
C     .. Intrinsic Functions ..
C      INTRINSIC ACOSD
      INTRINSIC ACOS
C     ..
C     .. External Subroutines ..
      EXTERNAL UNIT,VDIF
C     ..
C     .. Statement Functions ..
      REAL ACOSD
C     ..
C     .. Statement Function definitions ..
      ACOSD(X) = ACOS(X)*57.29578
C     ..
C
C
      CALL VDIF(A,X1,X2)
      CALL VDIF(B,X3,X2)
      CALL UNIT(A)
      CALL UNIT(B)
      CANG = DOT(A,B)
      ANG = ACOSD(CANG)
C
      END
C
C
C
C      =====================================================
      SUBROUTINE GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZA,IERR)
C      =====================================================
C
C----  find next atom with name ATOM in file,
C      search forwards and backwards
C      for IR=1 but only
C      backwards for IR=-1, Residue number must be MsourceResid,
C      and if ModeFlag is IMOL
C      symmetry operation must be ISYM
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalar Arguments ..
      INTEGER IA,IERR,IR,ISYM,MsourceResid,NCH
      CHARACTER ATOM*4
C     ..
C     .. Array Arguments ..
      REAL XYZA(3)
C     ..
C     .. Scalars in Common ..
      REAL Dclose,Dlimit,FRdist
      INTEGER Icount,Icounts,Idents,NumSymm,
     +        IsitAtom,IsitResidue,IsitChains
      CHARACTER ModeAtomFlag*2,ModeFlag*2,Title*80
      CHARACTER ChainSelected*1
      LOGICAL DoChains,Hexclude
C     ..
C     .. Arrays in Common ..
      INTEGER KtoAtmBegin,KtoAtmEnd,KFromAtmBegin,KFromAtmEnd,
     +        KtoResidSelect,KFromResidSelect,KSelectAtmType
      REAL Batom,FRclose,RFin,Rsymm,S,XYZ,XYZS
      INTEGER*2 InumAtmSymm,IresidNums,Nsymop,NUM,NUS
      INTEGER*2 FromAtmSelected,ToAtomSelected
      LOGICAL SourceResFlag,TargetResFlag
      CHARACTER ChainNames*1,AtomNames*4,ResidNames*4,SymmTitles*80
C     ..
C     .. Local Scalars ..
      INTEGER I,IP,J
C     ..
C     .. Common blocks ..
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /InputNames/AtomNames(MaxAtoms),ResidNames(MaxAtoms),
     +       ChainNames(MaxAtoms)
      COMMON /Params/NUM(MaxAtoms),XYZ(3,MaxAtoms),Icount,
     +       Batom(MaxAtoms)
      COMMON /ResNumCom/IresidNums(MaxAtoms),
     +                  FromAtmSelected(MaxAtoms),
     +                  ToAtomSelected(MaxAtoms)
      COMMON /ATMCHR/ChainSelected(MaxSelectedChains,2)
      COMMON /SelectControl/IsitAtom,IsitResidue,IsitChains,DoChains,
     +       KtoAtmBegin(MaxSelectedAtoms),
     +       KtoAtmEnd(MaxSelectedAtoms),
     +       KFromAtmBegin(MaxSelectedAtoms),
     +       KFromAtmEnd(MaxSelectedAtoms),
     +       KtoResidSelect(MaxSelectedResidues),
     +       KFromResidSelect(MaxSelectedResidues),
     +       KSelectAtmType(MaxSelectedResidues)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
      IERR = 0
   10 CONTINUE
C
      J = IA
C
C----  atoms are stored with symmetry number changing fastest, hence
C      limits on DO loop. 17 is for 14 atoms in trp + 3 m/C
C ccx ????? what about residues > trp in num atoms
C
      DO 20 I = 1,17*NumSymm
        J = J + IR
C
C----  test for start of coordinates
C
        IF (J .lt. 1) THEN
          GO TO 30
        ELSE
C
          IF (ModeFlag .eq. 'IM') THEN
            IP = InumAtmSymm(J)
          ELSE
            IP = J
          END IF
C
C----  test for end of coordinates
C
          IF (IP .eq. 0) THEN
            GO TO 30
          ELSE IF (AtomNames(IP) (1:NCH) .eq. ATOM(1:NCH)) THEN
            IF (IresidNums(IP) .eq. MsourceResid) THEN
              IF (ModeFlag .ne. 'IM' .or. 
     +             Nsymop(J) .eq. ISYM) GO TO 40
            END IF
          END IF
        END IF
   20 CONTINUE
C
   30 IF (IR .eq. 1) THEN
        IR = -1
        GO TO 10
      ELSE
        GO TO 60
      END IF
C
   40 CONTINUE

C
      DO 50 I = 1,3
        IF (ModeFlag .eq. 'IM') THEN
          XYZA(I) = XYZS(I,J)
        ELSE
          XYZA(I) = XYZ(I,J)
        END IF
   50 CONTINUE
C
      RETURN
   60 IERR = -1
      write (6,fmt=6000) ATOM,IA,MsourceResid
 6000 FORMAT (' Cant find ATOM ',A,' for ATOM',I5,' Residue',I5)
C
      END
C
C
C     =================================================
      SUBROUTINE GETH(SATOM,RES,MsourceResid,IA,ISYM,XS,XN,ANG)
C     =================================================
C
C     .. Scalar Arguments ..
      REAL ANG
      INTEGER IA,ISYM,MsourceResid
      CHARACTER RES*4,SATOM*4
C     ..
C     .. Array Arguments ..
      REAL XN(3),XS(3)
C     ..
C     .. Local Scalars ..
      INTEGER I,IERR,IR,J,JRES,NCH
      CHARACTER ATOM*4,ATOM1*4,ATOM2*4
C     ..
C     .. Local Arrays ..
      REAL DSQ(2),XH(3),XH2(3,2),XYZ1(3),XYZ2(3),XYZ3(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL GETAMH,GETANGL,GETATOM,GETNH2
C     ..
C
      ANG = 0.0
C
C----   test for atom type. Possibilities are ...
C       single hydrogen... N (m/c), NE1 (trp), NE (arg)
C       NH2............... ND2 (asn), NE2 (gln), NH1,nh2 (arg)
C
C----  main chain NH
C
      NCH = 4
C
      IF (SATOM .eq. 'N   ') THEN
C
C----   find coords of carbonyl carbon of preceeding Residue and
C       alpha carbon of this Residue
C
        IR = -1
        ATOM = 'C   '
        JRES = MsourceResid - 1
        IF (JRES .ge. 1) THEN
          CALL GETATOM(ATOM,NCH,JRES,IA,IR,ISYM,XYZ1,IERR)
          IF (IERR .ne. -1) THEN
            IR = 1
            ATOM = 'CA  '
            CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ3,IERR)
            IF (IERR .ne. -1) THEN
              CALL GETAMH(XYZ1,XN,XYZ3,XH)
              CALL GETANGL(XS,XH,XN,ANG)
            END IF
          END IF
        END IF
C
C----  trp, NE on arg
C
      ELSE IF (RES .eq. 'TRP ' .or. (RES .eq. 'ARG '.and.
     +         SATOM .eq. 'NE  ')) THEN
        ATOM = 'CE2 '
        IF (RES .eq. 'ARG ') ATOM = 'CZ  '
        IR = 1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .ne. -1) THEN
          ATOM = 'CD1 '
          IF (RES .eq. 'ARG ') ATOM = 'CD  '
          IR = 1
          CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ3,IERR)
          IF (IERR .ne. -1) THEN
            CALL GETAMH(XYZ1,XN,XYZ3,XH)
            CALL GETANGL(XS,XH,XN,ANG)
          END IF
        END IF
C
C----  ND1 or NE2 on His (assume H lies on whichever is in the contact)
C
      ELSE IF (RES .eq. 'HIS ') THEN
        IF (SATOM .eq. 'NE2 ') THEN
          ATOM1 = 'CD2 '
          ATOM2 = 'CE1 '
        ELSE
          ATOM1 = 'CE1 '
          ATOM2 = 'CG  '
        END IF
        IR = 1
        CALL GETATOM(ATOM1,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .ne. -1) THEN
          IR = 1
          CALL GETATOM(ATOM2,NCH,MsourceResid,IA,IR,ISYM,XYZ3,IERR)
          IF (IERR .ne. -1) THEN
            CALL GETAMH(XYZ1,XN,XYZ3,XH)
            CALL GETANGL(XS,XH,XN,ANG)
          END IF
        END IF
C
C----  arg NH2
C
      ELSE IF (RES .eq. 'ARG ') THEN
        IF (SATOM .eq. 'NH1 ') THEN
          ATOM = 'NH2 '
          IR = 1
        ELSE
          ATOM = 'NH1 '
          IR = 1
        END IF
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ3,IERR)
        IF (IERR .ne. -1) THEN
          ATOM = 'CZ  '
          IR = 1
          CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ2,IERR)
          IF (IERR .ne. -1) THEN
            CALL GETNH2(XN,XYZ2,XYZ3,XH2)
            DO 20 I = 1,2
              DSQ(I) = 0.0
              DO 10 J = 1,3
                DSQ(I) = (XS(J)-XH2(J,I))**2 + DSQ(I)
   10         CONTINUE
   20       CONTINUE
            I = 1
            IF (DSQ(2) .lt. DSQ(1)) I = 2
            CALL GETANGL(XS,XH2(1,I),XN,ANG)
          END IF
        END IF
C
C----  asn or gln
C
      ELSE IF (RES .eq. 'ASN ' .or. RES .eq. 'GLN ') THEN
        IF (RES .eq. 'ASN ') THEN
          ATOM1 = 'CG  '
          ATOM2 = 'OD1 '
        ELSE
          ATOM1 = 'CD  '
          ATOM2 = 'OE1 '
        END IF
        IR = 1
        CALL GETATOM(ATOM1,NCH,MsourceResid,IA,IR,ISYM,XYZ2,IERR)
        IF (IERR .ne. -1) THEN
          CALL GETATOM(ATOM2,NCH,MsourceResid,IA,IR,ISYM,XYZ3,IERR)
          IF (IERR .ne. -1) THEN
            CALL GETNH2(XN,XYZ2,XYZ3,XH2)
C
C----  find H nearest source atom
C
            DO 40 I = 1,2
              DSQ(I) = 0.0
              DO 30 J = 1,3
                DSQ(I) = (XS(J)-XH2(J,I))**2 + DSQ(I)
   30         CONTINUE
   40       CONTINUE
            I = 1
            IF (DSQ(2) .lt. DSQ(1)) I = 2
            CALL GETANGL(XS,XH2(1,I),XN,ANG)
          END IF
        END IF
      ELSE
        write (6,fmt=6000) SATOM,RES
      END IF
 6000 FORMAT (' Cannot calculate H position for ATOM ',A,' Residue ',A)
C
      END
C
C
C     ===============================
      SUBROUTINE GETNH2(X1,X2,X3,XH2)
C     ===============================
C
C----  calculate coords for NH2 hydrogens, nitrogen in X1
C
C     .. Array Arguments ..
      REAL X1(3),X2(3),X3(3),XH2(3,2)
C     ..
C     .. Local Arrays ..
      REAL A(3),B(3),C(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL GETAMH,UNIT,VSUM,VDIF
C     ..
C
      CALL VDIF(A,X1,X2)
      CALL UNIT(A)
      CALL VDIF(B,X3,X2)
      CALL UNIT(B)
      CALL VSUM(C,A,B)
      CALL UNIT(C)
      CALL VSUM(XH2(1,1),X1,C)
      CALL GETAMH(X2,X1,XH2(1,1),XH2(1,2))
      END
C
C
C     =================================================
      SUBROUTINE GETO(SATOM,RES,MsourceResid,IA,ISYM,XS,XO,ANG)
C     =================================================
C
C     .. Scalar Arguments ..
      REAL ANG
      INTEGER IA,ISYM,MsourceResid
      CHARACTER RES*4,SATOM*4
C     ..
C     .. Array Arguments ..
      REAL XO(3),XS(3)
C     ..
C     .. Local Scalars ..
      INTEGER IERR,IR,NCH
      CHARACTER ATOM*4
C     ..
C     .. Local Arrays ..
      REAL XYZ1(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL GETANGL,GETATOM
C     ..
C
      ANG = 0.0
      NCH = 4
C
C----  See if this is a main-chain carbonyl oxygen or side-chain
C
      IF (SATOM .eq. 'O   ') THEN
        ATOM = 'C   '
        IR = -1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .eq. -1) RETURN
C
C----  test for sIdechain type
C      asp,asn get CG
C
      ELSE IF (RES .eq. 'ASP ' .or. RES .eq. 'ASN ') THEN
        ATOM = 'CG  '
        IR = 1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .eq. -1) RETURN
C
C----  glu,gln get CD
C
      ELSE IF (RES .eq. 'GLU ' .or. RES .eq. 'GLN ') THEN
        ATOM = 'CD  '
        IR = 1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .eq. -1) RETURN
C
C----  thr,ser get CB
C
      ELSE IF (RES .eq. 'THR ' .or. RES .eq. 'SER ') THEN
        ATOM = 'CB  '
        IR = 1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .eq. -1) RETURN
C
C----  tyr get CZ
C
      ELSE IF (RES .eq. 'TYR ') THEN
        ATOM = 'CZ  '
        IR = 1
        CALL GETATOM(ATOM,NCH,MsourceResid,IA,IR,ISYM,XYZ1,IERR)
        IF (IERR .eq. -1) RETURN
      ELSE
        write (6,fmt=6000) RES
        RETURN
      END IF
C
C----  calculate angle
C
      CALL GETANGL(XS,XO,XYZ1,ANG)
 6000 FORMAT (1X,'UNrecognised Residue Type ',A)
C
      END
C
C
C     ====================
      SUBROUTINE WaterInit
C     ====================
C
C     .. Arrays in Common ..
      REAL ANG,BTYP,D,DistMin,DMAX,DSQ
      INTEGER NANG,NB,NBD
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
C     .. Common blocks ..
      COMMON /SUM/D(5),DistMin(5),DMAX(5),DSQ(5),NB(5),NBD(11,5),
     +       NANG(11,5),ANG(11,5),BTYP(5)
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
      DO 20 I = 1,5
        D(I) = 0.0
        DistMin(I) = 100.0
        DMAX(I) = 0.0
        DSQ(I) = 0.0
        NB(I) = 0
C
        DO 10 J = 1,11
          NBD(J,I) = 0
          NANG(J,I) = 0
          ANG(J,I) = 0
   10   CONTINUE
   20 CONTINUE
C
      END
C
C
C     =======================================================
      SUBROUTINE SETNET(Xmin,Ymin,Zmin,NX,NY,NZ,NU,XY,IC,NET)
ccx         CALL SETNET(Xmin,Ymin,Zmin,NX,NY,NZ,NUS,XYZS,Icounts,NET)
C     =======================================================
C
      INTEGER MaxAtoms
      PARAMETER (MaxAtoms=48000)
      REAL Xmin,Ymin,Zmin
      INTEGER IC,NX,NY,NZ
      REAL XY(3,*)
      INTEGER*2 NET(51,NX,NY,NZ),NU(*)
      REAL X,Y,Z
      INTEGER I,IN,INMAX,IX,IY,IZ,J,K,Error_counter
      INTRINSIC INT
      EXTERNAL CCPERR
      data Error_counter/0/
      INMAX = -1
      DO 30 I = 1,NX
        DO 20 J = 1,NY
          DO 10 K = 1,NZ
            NET(1,I,J,K) = 0
   10     CONTINUE
   20   CONTINUE
   30 CONTINUE
C
C
      DO 40 I = 1,IC
        X = XY(1,I) - Xmin
        Y = XY(2,I) - Ymin
        Z = XY(3,I) - Zmin
C
        IF (Xmin .lt. 0.0) THEN
          IX = NINT( (X+6.1)/6.0) + 1
        ELSE
           IX = INT( X/6.0) + 2
        END IF
        IF (Ymin .lt. 0.0) THEN
           IY = NINT( (Y+6.1)/6.0) + 1
        ELSE
           IY = INT( Y/6.0) + 2
        END IF
        IF (Zmin .lt. 0.0) THEN
          IZ = NINT( (Z+6.1)/6.0) + 1
        ELSE
           IZ = INT( Z/6.0) + 2
        END IF
C
        NET(1,IX,IY,IZ) = NET(1,IX,IY,IZ) + 1
        IN = NET(1,IX,IY,IZ)
        IF (INMAX .lt. IN) INMAX = IN
c
c
        IF (IN+1 .gt. 51) THEN
ccx           write (6,6666) IN,IX,IY,IZ
 6666    format(4i10)
             error_counter=error_counter+1
ccx             write (6,8000) 'Num atoms per brick exceeded'
 8000       format(a)
        NET(1,IX,IY,IZ) = NET(1,IX,IY,IZ) - 1
        inmax = inmax-1
ccx         CALL CCPERR(1,' Num atoms per brick exceeded ')
         else
        NET(IN+1,IX,IY,IZ) = NU(I)
          end if
c
c
   40 CONTINUE
C
cpjx INMAX can never be greater than 50 because each time the brick limit
cpjx is exceeded, INMAX is decremented by 1
cpjx      IF (INMAX .gt. 51) THEN
      IF (error_counter .GT. 0) THEN
        write (6,fmt=6000) (INMAX+error_counter)
 6000   FORMAT (/,
     +' *** there are more than 50 atoms (',i3,') in one 6x',
     +         '6x6 A box',/,
     +' - this is unrealistic - check your matrices',/)
        CALL CCPERR(1,'Bad matrices')
      END IF
cpjx This line is redundant now - it is unreachable if error_couter is
cpjx greater than zero.
cpjx        if (error_counter.gt.50) call ' BRICK EXCEEDED over 50 times'
C
      END
C
C
C     =============================
      SUBROUTINE SIDENT(XMAT,Iflag)
C     =============================
C
C  test for unit matrix and zero translation
C
C     .. Scalar Arguments ..
      INTEGER Iflag
C     ..
C     .. Array Arguments ..
      REAL XMAT(12)
C     ..
C     .. Local Scalars ..
      REAL XSUM
      INTEGER I
C     ..
C
      Iflag = 0
      IF (XMAT(1) .eq. 1.0 .and. XMAT(6) .eq. 1.0 .and.
     +    XMAT(11) .eq. 1.0) THEN
        XSUM = 0.0
C
        DO 10 I = 1,12
          XSUM = XMAT(I) + XSUM
   10   CONTINUE
C
        IF (XSUM .eq. 3.0) Iflag = 1
      END IF
C
      END
C
C
C       ====================
      SUBROUTINE WaterStats
C       ====================
C
C     .. Arrays in Common ..
      REAL ANG,BTYP,D,DistMin,DMAX,DSQ
      INTEGER NANG,NB,NBD
C     ..
C     .. Local Scalars ..
      REAL ANGMEAN,BMEAN,DIS,RMS,SUMANG
      INTEGER I,K,N,NTOT
C     ..
C     .. Local Arrays ..
      CHARACTER CTYPE(5)*6
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC sqrt
C     ..
C     .. Common blocks ..
      COMMON /SUM/D(5),DistMin(5),DMAX(5),DSQ(5),NB(5),NBD(11,5),
     +       NANG(11,5),ANG(11,5),BTYP(5)
C     ..
C     .. Save statement ..
      SAVE
C     ..
C     .. Data statements ..
      DATA CTYPE/'M/C O','M/C N','S/C O','S/C N','WAT  '/
C     ..
C
      write (6,fmt=6000)
 6000 FORMAT (/,
     +' Analysis of Water Hydrogen Bonds',/,
     +' Atom  Number    <D>    RMS    MIN    MAX  <Angle>   <H2O B>')
C
      DO 20 I = 1,5
        DIS = 0.0
        RMS = 0.0
        SUMANG = 0.0
        ANGMEAN = 0.0
        NTOT = 0
        N = NB(I)
        IF (N .ne. 0) THEN
C
          DO 10 K = 1,11
            IF (NANG(K,I) .ne. 0) THEN
              SUMANG = ANG(K,I) + SUMANG
              NTOT = NANG(K,I) + NTOT
              ANG(K,I) = ANG(K,I)/NANG(K,I)
            END IF
   10     CONTINUE
C
          IF (NTOT .ne. 0) ANGMEAN = SUMANG/NTOT
          DIS = D(I)/N
          RMS = sqrt(DSQ(I)/N-DIS*DIS)
          BMEAN = BTYP(I)/N
        END IF
C
        write (6,fmt=6002) CTYPE(I),N,DIS,RMS,DistMin(I),DMAX(I),
     +    ANGMEAN,BMEAN
 6002   FORMAT (1X,A,I8,4F7.2,2F8.1)
   20 CONTINUE
C
      DO 30 I = 1,5
        write (6,fmt=6004) CTYPE(I), (NBD(K,I),K=1,11),
     +    (ANG(K,I),K=1,11)
 6004   FORMAT (/,
     +' Distribution of Bonds for ',A,' Atoms',/,
     +' Distance 2.4   2.5   2.6   2.7   2.8   ',
     +'2.9   3.0   3.1   3.2   ',
     +    '3.3   3.4   3.5',/,
     +' Number  ',11I6,/1X,'Angle    ',11F6.0)
   30 CONTINUE
C
C
      END
C
C
C     =====================================
      SUBROUTINE WaterSums(N,DIS,ANGL,BOLD)
C     =====================================
C
C     .. Scalar Arguments ..
      REAL ANGL,BOLD,DIS
      INTEGER N
C     ..
C     .. Arrays in Common ..
      REAL ANG,BTYP,D,DistMin,DMAX,DSQ
      INTEGER NANG,NB,NBD
C     ..
C     .. Local Scalars ..
      INTEGER IHIS
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MAX,MIN
C     ..
C     .. Common blocks ..
      COMMON /SUM/D(5),DistMin(5),DMAX(5),DSQ(5),NB(5),NBD(11,5),
     +       NANG(11,5),ANG(11,5),BTYP(5)
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
      IHIS = (DIS-2.4)*10 + 1
      IF (IHIS .lt. 1) IHIS = 1
      IF (IHIS .gt. 11) IHIS = 11
      NBD(IHIS,N) = NBD(IHIS,N) + 1
C
      IF (ANGL .ne. 0.0) THEN
        NANG(IHIS,N) = NANG(IHIS,N) + 1
        ANG(IHIS,N) = ANG(IHIS,N) + ANGL
      END IF
C
      NB(N) = NB(N) + 1
      D(N) = D(N) + DIS
      DistMin(N) = MIN(DIS,DistMin(N))
      DMAX(N) = MAX(DIS,DMAX(N))
      DSQ(N) = DSQ(N) + DIS*DIS
      BTYP(N) = BTYP(N) + BOLD
C
      END
C
C
C     =================
      SUBROUTINE SYMGEN
C     =================
C     Generate symmetry related atoms
C     NUS([atom index < MaxAtoms=48000]) internal number of symmetry generated atoms
C     InumAtmSymm([atom index]) -> internal atom number NUM([atom index]) of the
C       original atom from which it was generated
C     XYZS([dimension <=3],[atom index]) coordinates of generated atom
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalars in Common ..
      REAL Dclose,DistMetalLigands,Dlimit,FRdist
      INTEGER Icount,Icounts,Idents,MetalCounter,NumSymm
      LOGICAL AutoFlag,LookForMetal,Hexclude,BigSearch
      CHARACTER ModeAtomFlag*2,ModeFlag*2,Title*80
C     ..
C     .. Arrays in Common ..

      REAL Batom,FRclose,RFin,Rsymm,S,TR,XYZ,XYZbounds,XYZS
      INTEGER NumMetalAngles
      INTEGER*2 InumAtmSymm,Nsymop,NUM,NUS,TranslationCode
      LOGICAL SourceResFlag,TargetResFlag
      CHARACTER TranslationTitles*9,SymmTitles*80
C     ..
C     .. Local Scalars ..
      INTEGER I,J,K,L,LT,NTRANS,NIDENT
C     ..
C     .. Local Arrays ..
      REAL LocalBounds(2,3),X(3)
C     ..
C     .. Common blocks ..
      CHARACTER ChainNames*1,AtomNames*4,ResidNames*4
      COMMON /InputNames/AtomNames(MaxAtoms),ResidNames(MaxAtoms),
     +       ChainNames(MaxAtoms)
      COMMON /AutoCom/AutoFlag,TranslationCode(2*MaxAtoms),
     +       TR(4,4,MaxSymmetry),BigSearch
ccx      COMMON /BoundsCom/ Xmin,Xmax,Ymin,Ymax,Zmin,Zmax
      COMMON /BoundsCom/XYZbounds(2,3)
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /CharTT/TranslationTitles(125)
      COMMON /DatMetals/MetalCounter,DistMetalLigands,LookForMetal,
     +       NumMetalAngles(MaxMetals)
      COMMON /Params/NUM(MaxAtoms),XYZ(3,MaxAtoms),Icount,
     +       Batom(MaxAtoms)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
      Icounts = 0
      DO 10 K = 1,3
        LocalBounds(1,K) = XYZbounds(1,K) - Dlimit
        LocalBounds(2,K) = XYZbounds(2,K) + Dlimit
   10 CONTINUE
C
C---- AutoFlag ModeFlag...
C
      IF (AutoFlag) THEN
        DO 110 I = 1,Icount
C     Icount is from COMMON Params -> number of atoms selected
C
C---- Set number of translations required in total
C
           IF (BigSearch) THEN
             NTRANS = 125
             NIDENT = 63
           ELSE
             NTRANS = 27
             NIDENT = 14
           END IF
C
C---- simple translations first
C
C<<          DO 40 LT = NumSymm + 1,NumSymm + 27
           DO 40 LT = NumSymm + 1,NumSymm + NTRANS
C>> modified HP990407
            IF (LT-NumSymm .eq. NIDENT .and. 
     +          ( .not. LookForMetal) ) GO TO 40
            DO 20 K = 1,3
              X(K) = XYZ(K,I) + S(4,K,LT)
              IF (X(K) .lt. LocalBounds(1,K) .or.
     +            X(K) .gt. LocalBounds(2,K)) GO TO 40
   20       CONTINUE
            Icounts = Icounts + 1
            InumAtmSymm(Icounts) = NUM(I)
            NUS(Icounts) = Icounts
            DO 30 L = 1,3
              XYZS(L,Icounts) = X(L)
 30          CONTINUE
C
            Nsymop(Icounts) = NumSymm + 1
C
C---- lattice translation code
C
            TranslationCode(Icounts) = LT - NumSymm
   40     CONTINUE
C
C---- Symmetry + translations...
C
          DO 100 J = 1,NumSymm
C
C---- IDENTITY (translations already done)
C
            IF (J .eq. Idents) THEN
C
              DO 50 K = 1,3
                X(K) = XYZ(1,I)*S(1,K,J) + XYZ(2,I)*S(2,K,J) +
     +                 XYZ(3,I)*S(3,K,J) + S(4,K,J)
ccx           IF (X(K)  .ge.  LocalBounds(1,K) .and.
ccx     +           X(K)  .le.  LocalBounds(2,K)) GO TO 50
   50         CONTINUE
              Icounts = Icounts + 1
              InumAtmSymm(Icounts) = NUM(I)
              NUS(Icounts) = Icounts
              DO 60 L = 1,3
                XYZS(L,Icounts) = X(L)
   60         CONTINUE
              Nsymop(Icounts) = J
C
C---- lattice translation code
C
              TranslationCode(Icounts) = NIDENT
            ELSE
C
C----  other symmetries
C
C<<              DO 90 LT = NumSymm + 1,NumSymm + 27
               DO 90 LT = NumSymm + 1,NumSymm + NTRANS
C>> modified HP990407
                DO 70 K = 1,3
                  X(K) = XYZ(1,I)*S(1,K,J) + XYZ(2,I)*S(2,K,J) +
     +                   XYZ(3,I)*S(3,K,J) + S(4,K,J) + S(4,K,LT)
                  IF (X(K) .lt. LocalBounds(1,K) .or.
     +                X(K) .gt. LocalBounds(2,K)) GO TO 90
   70           CONTINUE
                Icounts = Icounts + 1
                InumAtmSymm(Icounts) = NUM(I)
                NUS(Icounts) = Icounts
                DO 80 L = 1,3
                  XYZS(L,Icounts) = X(L)
   80           CONTINUE
C
                Nsymop(Icounts) = J
C
C---- lattice translation code
C
                TranslationCode(Icounts) = LT - NumSymm
   90         CONTINUE
            END IF
  100     CONTINUE
  110   CONTINUE
      ELSE
C
C---- normal IMOL ModeFlag...
C
        DO 150 I = 1,Icount
          DO 140 J = 1,NumSymm
            DO 120 K = 1,3
              X(K) = XYZ(1,I)*S(1,K,J) + XYZ(2,I)*S(2,K,J) +
     +               XYZ(3,I)*S(3,K,J) + S(4,K,J)
              IF (X(K) .lt. LocalBounds(1,K) .or.
     +            X(K) .gt. LocalBounds(2,K)) GO TO 140
  120       CONTINUE
            Icounts = Icounts + 1
            InumAtmSymm(Icounts) = NUM(I)
            NUS(Icounts) = Icounts
            DO 130 L = 1,3
              XYZS(L,Icounts) = X(L)
  130       CONTINUE
            Nsymop(Icounts) = J
            TranslationCode(Icounts) = NIDENT
  140     CONTINUE
  150   CONTINUE
      END IF
      write (6,fmt=6002) Icounts
 6002 FORMAT (//' No. of symmetry related atoms in all boxes :',I6)
C
      END
C
C
C     =================
      SUBROUTINE SYMMAT
C     =================
C
C----  This subroutine prepares the [s] matrices (maximum MaxSymmetry) 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     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalars in Common ..
      REAL Dclose,Dlimit,FRdist
      INTEGER Icounts,Idents,NumSymm
      LOGICAL AutoFlag,Hexclude,BigSearch
      CHARACTER ModeAtomFlag*2,ModeFlag*2,Title*80
C     ..
C     .. Arrays in Common ..
      REAL FRclose,RFin,Rsymm,S,TR,XYZS
      INTEGER*2 InumAtmSymm,Nsymop,NUS,TranslationCode
      LOGICAL SourceResFlag,TargetResFlag
      CHARACTER TranslationTitles*9,SymmTitles*80
C     ..
C     .. Local Scalars ..
      REAL VOL
      INTEGER I,IADD,J,K
C     ..
C     .. Local Arrays ..
      REAL BRKMAT(4,4),cell(6),IBRKMAT(4,4),MAT1(4,4),MAT2(4,4),
     +     MAT3(4,4),RF(4,4),RO(4,4)
C     ..
C     .. External Subroutines ..
      EXTERNAL MATMUL4,RBcell,RBRORF
C     ..
C     .. Common blocks ..
      COMMON /AutoCom/AutoFlag,TranslationCode(2*MaxAtoms),
     +       TR(4,4,MaxSymmetry),BigSearch
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /CharTT/TranslationTitles(125)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
C
C---- the CRYST and SCALE records from BRK file :
C
      DO 10 I = 1,MaxSymmetry
        S(4,4,I) = 1.0
        TR(4,4,I) = 1.0
   10 CONTINUE
C
      BRKMAT(4,4) = 1.0
      RO(1,1) = 0.0
      RF(1,1) = 0.0
C
      CALL RBcell(cell,VOL)
      CALL RBRORF(RO,RF)
C
      DO 30 I = 1,4
        DO 20 J = 1,4
          BRKMAT(J,I) = RF(I,J)
          IBRKMAT(J,I) = RO(I,J)
   20   CONTINUE
   30 CONTINUE
C
C----   matrix multiplications:  [S] = [BRKMAT] * [TR] * [IBRKMAT]
C
      IADD = 0
      IF (AutoFlag) THEN
        IF (BigSearch) THEN
          IADD = 125
        ELSE
          IADD = 27
        END IF
      END IF
C
C
      DO 80 I = 1,NumSymm + IADD
        DO 50 J = 1,4
          DO 40 K = 1,4
            MAT2(K,J) = TR(K,J,I)
   40     CONTINUE
   50   CONTINUE
        CALL MATMUL4(MAT1,BRKMAT,MAT2)
        CALL MATMUL4(MAT3,MAT1,IBRKMAT)
        DO 70 J = 1,4
          DO 60 K = 1,4
            S(K,J,I) = MAT3(K,J)
   60     CONTINUE
   70   CONTINUE
   80 CONTINUE
C
      RETURN
C
C---- Format statements
C
C
      END
C
C
C     =======================================
      SUBROUTINE TRsymm(NumSymm,TR,BigSearch)
C     =======================================
C
C     .. Parameters ..
      INTEGER MaxSymmetry
      PARAMETER (MaxSymmetry=230)
C     ..
C     .. Scalar Arguments ..
      INTEGER NumSymm
      LOGICAL BigSearch
C     ..
C     .. Array Arguments ..
      REAL TR(4,4,MaxSymmetry)
C     ..
C     .. Local Scalars ..
      INTEGER I,J,K,N,NVEC,OFFS
C     ..
C     .. Local Arrays ..
      REAL TEMPL(4,4)
C     ..
C     .. Data statements ..
C
      DATA TEMPL/1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,
     +     0.0,0.0,1.0/
C     ..
C
      N = NumSymm
C
      DO 30 I = N + 1,N + 125
        DO 20 J = 1,4
          DO 10 K = 1,4
            TR(K,J,I) = TEMPL(K,J)
   10     CONTINUE
   20   CONTINUE
   30 CONTINUE
C
CPJX - For normal search I,J,K go 1 to 3, and TR(4,1,N)=float(I-2) etc
C      For BigSearch I,J,k go 1 to 5 and TR(4,1,N)=float(I-3) etc
C
C - this sets the parameters for generating the required lattice
C   translations
C
      IF (BigSearch) THEN
        NVEC = 5
        OFFS = 3
      ELSE
        NVEC = 3
        OFFS = 2
      END IF
C
      DO 60 I = 1,NVEC
        DO 50 J = 1,NVEC
          DO 40 K = 1,NVEC
            N = N + 1
C
ccx - check N to MaxSymmetry
            IF (N .gt. MaxSymmetry) CALL CCPERR(1,
     +      ' Maxsymmetry exceeded in routine TRSYMM')
C
C - this generates the required lattice translations
            TR(4,1,N) = float(I - OFFS)
            TR(4,2,N) = float(J - OFFS)
            TR(4,3,N) = float(K - OFFS)
   40     CONTINUE
   50   CONTINUE
   60 CONTINUE
C
      END
C
C
C     =============================================
      SUBROUTINE WATANG(IS,IR,IAPNT,ITRES,NC,ANGLW)
C     =============================================
C
C----  analyse bond angles at water by bonded atoms
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalar Arguments ..
      REAL ANGLW
      INTEGER IR,IS,NC
C     ..
C     .. Array Arguments ..
      INTEGER IAPNT(10),ITRES(10)
C     ..
C     .. Scalars in Common ..
      REAL Dclose,Dlimit,FRdist
      INTEGER Icount,Icounts,Idents,numBAD,NumSymm,TotalWaterAngles
      CHARACTER ModeAtomFlag*2,ModeFlag*2,Title*80
      LOGICAL Hexclude
C     ..
C     .. Arrays in Common ..
      REAL Batom,FRclose,RFin,Rsymm,S,XYZ,XYZS
      INTEGER BadWaters,HistogramAngles
      INTEGER*2 InumAtmSymm,Nsymop,NUM,NUS
      LOGICAL SourceResFlag,TargetResFlag
      CHARACTER SymmTitles*80
C     ..
C     .. Local Scalars ..
      REAL ANG,SANG
      INTEGER I,IHIS,J,K,NANG
C     ..
C     .. Local Arrays ..
      REAL X1(3),X2(3)
C     ..
C     .. External Subroutines ..
      EXTERNAL GETANGL
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC NINT
C     ..
C     .. Common blocks ..
      COMMON /BADCOMM/numBAD,BadWaters(4,MaxResidues),
     +       HistogramAngles(13),TotalWaterAngles
      COMMON /CardsCom/SourceResFlag(MaxResidues),
     +       TargetResFlag(MaxResidues),Dlimit,Dclose,Hexclude
      COMMON /CharSymm/SymmTitles(MaxSymmetry)
      COMMON /Params/NUM(MaxAtoms),XYZ(3,MaxAtoms),Icount,
     +       Batom(MaxAtoms)
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
      COMMON /Xcards/ModeAtomFlag,ModeFlag,Title
C     ..
C     .. Save statement ..
      SAVE
C     ..
C
      ANGLW = 0.0
      SANG = 0.0
      NANG = 0
      DO 40 I = 1,NC - 1
        DO 10 K = 1,3
          IF (ModeFlag .eq. 'IM') THEN
C     IM -> intermolecular contacts
            X1(K) = XYZS(K,IAPNT(I))
          ELSE
            X1(K) = XYZ(K,IAPNT(I))
          END IF
   10   CONTINUE
        DO 30 J = I + 1,NC
          DO 20 K = 1,3
            IF (ModeFlag .eq. 'IM') THEN
              X2(K) = XYZS(K,IAPNT(J))
            ELSE
              X2(K) = XYZ(K,IAPNT(J))
            END IF
   20     CONTINUE
          CALL GETANGL(X1,XYZ(1,IS),X2,ANG)
C
C----  bin 0-40,40-50,50-60......140-150,150-180
C
          IHIS = (ANG-30.0)*0.1 + 1
          IF (IHIS .lt. 1) IHIS = 1
          IF (IHIS .gt. 13) IHIS = 13
          HistogramAngles(IHIS) = HistogramAngles(IHIS) + 1
          TotalWaterAngles = TotalWaterAngles + 1
          IF ((ANG .lt. 70.0) .or. (ANG .gt. 150)) THEN
            numBAD = numBAD + 1
            IF (numBAD .le. MaxResidues) THEN
              BadWaters(1,numBAD) = ITRES(I)
              BadWaters(2,numBAD) = IR
              BadWaters(3,numBAD) = ITRES(J)
              BadWaters(4,numBAD) = NINT(ANG)
            END IF
          ELSE
            NANG = NANG + 1
            SANG = SANG + ANG
          END IF
   30   CONTINUE
   40 CONTINUE
C
      IF (NANG .ne. 0) ANGLW = SANG/NANG
C
      END
C
C
C     ===============================================================
      SUBROUTINE WATER(IRESS,MtargetResid,ATYP,RES,D,ANGL,Iend,ISPNT,
     +                 ITPNT,Bwater)
C     ===============================================================
C
C     .. Parameters ..
      INTEGER MaxResidues
      PARAMETER (MaxResidues=9000)
C     ..
C     .. Scalar Arguments ..
      REAL ANGL,Bwater,D
      INTEGER Iend,IRESS,ISPNT,ITPNT,MtargetResid
      CHARACTER ATYP*4,RES*4
C     ..
C     .. Scalars in Common ..
      REAL AngleHY,AngleOx,Bmax,DistMin
      INTEGER maxNB,numBAD,TotalWaterAngles
C     ..
C     .. Arrays in Common ..
      INTEGER BadWaters,HistogramAngles
C     ..
C     .. Local Scalars ..
      REAL Anglim,ANGLW,BOLD,DMAX
      INTEGER I,IR,IROLD,ISOLD,J,K,N,NC,NMIN,NMORE,NPR,NPW,NSEC,NWAT
      LOGICAL FIRST,MORE,SECOND
      CHARACTER line*80
C     ..
C     .. Local Arrays ..
      REAL ANGLA(10),ANGLE(10),BCO(10),BPR(11,11),BPRTOT(11),
     +     DDistMin(MaxResidues),DIS(10)
      INTEGER IAPNT(10),IIRESS(MaxResidues),IIREST(MaxResidues),
     +        ITRES(10),NANGLE(10),NCO(10),NPIRES(MaxResidues,11,11),
     +        NPROT(11,11),NPTOT(11)
      CHARACTER ATYPA(10)*4,CTYPE(4)*4,RESA(10)*4
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR,WaterInit,WaterStats,WaterSums,WATANG
C     ..
C     .. Common blocks ..
      COMMON /BADCOMM/numBAD,BadWaters(4,MaxResidues),
     +       HistogramAngles(13),TotalWaterAngles
      COMMON /WatCom/AngleHY,AngleOx,DistMin,maxNB,Bmax
C     ..
C     .. Save statement ..
      SAVE
C     ..
C     .. Data statements ..
C
      DATA CTYPE/'O   ','N   ','OD1 ','ND1 '/
      DATA FIRST/.true./
C     ..
C
      MORE = .false.
C
C----  set minimum distance, maximum number of neighbours
C      remember CONTACT may give distances of 0.0 for atoms in special
C      positions
C
      IF (FIRST) THEN
        NMIN = 0
        FIRST = .false.
        IROLD = IRESS
        ISOLD = ISPNT
        BOLD = Bwater
        NC = 0
C
C            **********
        CALL WaterInit
C            **********
C
        DO 20 I = 1,11
          NPTOT(I) = 0
          BPRTOT(I) = 0.0
          IF (I .ne. 11) THEN
            NCO(I) = 0
            NANGLE(I) = 0
            ANGLE(I) = 0.0
            BCO(I) = 0.0
          END IF
C
          DO 10 J = 1,11
            BPR(J,I) = 0.0
            NPROT(J,I) = 0
   10     CONTINUE
   20   CONTINUE
C
C---- Test for new Residue
C
      ELSE IF (IRESS .ne. IROLD .or. Iend .eq. 1) THEN
C
C----  new Residue, deal with previous one
C      first limit number of neighbours to maxNB
C
        MORE = .false.
   30   CONTINUE
        IF (NC .gt. maxNB) THEN
          MORE = .true.
          DMAX = 0.0
C
          DO 40 I = 1,NC
            IF (DIS(I) .ge. DMAX) THEN
              DMAX = DIS(I)
              IR = I
            END IF
   40     CONTINUE
C
          J = 0
C
          DO 50 I = 1,NC
            IF (I .ne. IR) THEN
              J = J + 1
              DIS(J) = DIS(I)
              RESA(J) = RESA(I)
              ATYPA(J) = ATYPA(I)
              ITRES(J) = ITRES(I)
              ANGLA(J) = ANGLA(I)
              IAPNT(J) = IAPNT(I)
            END IF
   50     CONTINUE
C
          NC = NC - 1
          GO TO 30
        END IF
C
C----  compile sums for this source water
C
        SECOND = .true.
        IF (MORE) NMORE = NMORE + 1
        NPR = 1
        NPW = 1
C
        DO 90 I = 1,NC
C----  assign type of bond
C
          DO 60 J = 1,2
            IF (ATYPA(I) .eq. CTYPE(J)) GO TO 80
   60     CONTINUE
C
          DO 70 J = 3,4
            IF (ATYPA(I) (1:1) .eq. CTYPE(J) (1:1)) GO TO 80
   70     CONTINUE
C
          write (6,fmt=6000) ATYPA(I)
          GO TO 90
   80     IF (J .eq. 1 .and. RESA(I) .eq. 'WAT ') J = 5
          IF (J .eq. 1 .and. RESA(I) .eq. 'HOH ') J = 5
C
C----  count number of protein H bonds
C
          IF (J .lt. 5) NPR = NPR + 1
          IF (J .eq. 5) NPW = NPW + 1
C
C----   do not call WATERSUMS for water-water contacts where target Residue
C     number is less than source, to avoid counting the same bond twice.
C
          IF (J .ne. 5 .or. ITRES(I) .ge. IROLD) THEN
C
C                ****************************
            CALL WaterSums(J,DIS(I),ANGLA(I),BOLD)
C                ****************************
C
            IF (J .le. 4) SECOND = .false.
          END IF
   90   CONTINUE
C
C----  number of waters forming bonds and also those
C      in second coordination shell
C
        NWAT = NWAT + 1
        IF (SECOND) NSEC = NSEC + 1
        IF (NC.GT.0) THEN
           NCO(NC) = NCO(NC) + 1
           BCO(NC) = BCO(NC) + BOLD
        END IF
        NPROT(NPR,NPW) = NPROT(NPR,NPW) + 1
        BPR(NPR,NPW) = BPR(NPR,NPW) + BOLD
        NPIRES(NPROT(NPR,NPW),NPR,NPW) = IROLD
C
C----  analyse bond angles at water
C
C                    *********************************************
        IF (NC .gt. 1) CALL WATANG(ISOLD,IROLD,IAPNT,ITRES,NC,ANGLW)
C                    *********************************************
C
        IF ((ANGLW .ne. 0).AND.(NC.GT.0)) THEN
          ANGLE(NC) = ANGLE(NC) + ANGLW
          NANGLE(NC) = NANGLE(NC) + 1
        END IF
C
        NC = 0
        IROLD = IRESS
        ISOLD = ISPNT
        BOLD = Bwater
        IF (Iend .eq. 1) THEN
C
C----  end of contact list, do statistics
C
          IF (NMORE .ne. 0) write (6,fmt=6006) NMORE,maxNB,maxNB
          IF (NMIN .ne. 0) write (6,fmt=6008) NMIN,DistMin,
     +        (IIRESS(I),IIREST(I),DDistMin(I),I=1,NMIN)
C
          DO 100 I = 1,maxNB
            IF (NCO(I) .ne. 0) BCO(I) = BCO(I)/NCO(I)
  100     CONTINUE
C
          write (6,fmt=6010) NWAT,NSEC, (NCO(I),I=1,maxNB)
          write (6,fmt=6012) (BCO(I),I=1,maxNB)
C
          DO 120 I = 1,maxNB + 1
            DO 110 J = 1,maxNB + 1
              BPRTOT(I) = BPRTOT(I) + BPR(I,J)
              NPTOT(I) = NPTOT(I) + NPROT(I,J)
  110       CONTINUE
C
            IF (NPTOT(I) .ne. 0) BPRTOT(I) = BPRTOT(I)/NPTOT(I)
  120     CONTINUE
C
          write (6,fmt=6014) (I,I=0,maxNB)
          write (6,fmt=6016) (NPTOT(I),I=1,maxNB+1)
          write (6,fmt=6018) (BPRTOT(I),I=1,maxNB+1)
C
          DO 140 K = maxNB + 1,1,-1
            DO 130 N = 1,maxNB - K + 2
C
              IF (NPROT(K,N) .ne. 0) THEN
                write (6,fmt=6020) K - 1,N - 1,
     +            (NPIRES(I,K,N),I=1,NPROT(K,N))
                write (6,fmt=6022) BPR(K,N)/NPROT(K,N)
              END IF
  130       CONTINUE
  140     CONTINUE
C
C              ***********
          CALL WaterStats
C              ***********
C
          DO 150 I = 2,maxNB
            IF (NANGLE(I) .ne. 0) ANGLE(I) = ANGLE(I)/NANGLE(I)
  150     CONTINUE
C
          write (6,fmt=6024) (I,I=2,maxNB)
          write (6,fmt=6026) (ANGLE(I),I=2,maxNB)
C
          write (6,fmt=6028) TotalWaterAngles,HistogramAngles
C
          IF (numBAD .ne. 0) write (6,fmt=6030) numBAD,
     +        ((BadWaters(I,J),I=1,4),J=1,numBAD)
          RETURN
        END IF
      END IF
C
C----   store details of this interaction
C       omit distances less than DistMin
C
      IF (D .ge. DistMin) THEN
C
        IF (ANGL .ne. 0.0) THEN
C
C----  reject unacceptable angles
C
          IF (ATYP(1:1) .eq. 'N') THEN
            Anglim = AngleHY
          ELSE
            Anglim = AngleOx
          END IF
C
          IF (ANGL .lt. Anglim) RETURN
        END IF
        NC = NC + 1
        IF (NC .gt. 10) THEN
          write (6,fmt=6004) IRESS
          CALL CCPERR(1,'More than 10 Bonds for Residue ')
        ELSE
C
          DIS(NC) = D
          ATYPA(NC) = ATYP
          RESA(NC) = RES
          ANGLA(NC) = ANGL
          ITRES(NC) = MtargetResid
          IAPNT(NC) = ITPNT
        END IF
C
C----  trivial case of zero distance
C
      ELSE IF (D .ge. 0.01) THEN
        NMIN = NMIN + 1
C
        IF (NMIN .gt. MaxResidues) THEN
          write (6,fmt=6002) DistMin
          CALL CCPERR(1,'TOO Many bond distances less than Dmin')
        ELSE
C
          IIRESS(NMIN) = IRESS
          IIREST(NMIN) = MtargetResid
          DDistMin(NMIN) = D
        END IF
      END IF
C
C---- Format statements
C
 6000 FORMAT (' ** Unrecognised Atom Type **  ',A)
 6002 FORMAT (//,' *** ERROR, TOO Many bond distances less than ',
     +           'supplied  value for DistMin of',F6.2)
 6004 FORMAT (' *** More than 10 Bonds for Residue ',I5,' ***',/,
     +        ' Decrease Dlimit (MAX Contact distance)')
 6006 FORMAT (/,' ',I5,' Waters had more than',I3,' Bonds',/,
     +          ' Only the',I3,' Shortest were considered')
 6008 FORMAT (/,' ',I4,' Bonds were less than',F5.2,
     +          ' The Residues involved are listed below',/,
     +          (1X,5 (2I5,F5.2)))
 6010 FORMAT (/,' Number of Waters forming H Bonds',I5,/,
     +          ' Number in Second coordination shell',I5,/,
     +          ' Number of waters making 1,2,3,4 etc H Bonds',10I5)
 6012 FORMAT (/' Mean B',38X,10F5.0)
 6014 FORMAT (//,' Number of Waters making a given number of Hydrogen',
     +       ' bonds to protein',/,
     +       ' Number of bonds to protein ',11I5)
 6016 FORMAT (' Number of Waters',11X,11I5)
 6018 FORMAT (' Mean B',22X,11F5.0)
 6020 FORMAT (//,' The following Waters make',I3,
     +     ' Bonds to protein and',I3,' to Water',/ (1X,16I5))
 6022 FORMAT (/,' Mean B for this group of Waters',F6.0)
 6024 FORMAT (//,' Analysis of Bond Angles at Water atom',/,
     +' The following values exclude the rejected angles listed below',
     +       /,' Number of Bonds ',10I5)
 6026 FORMAT (' Mean Bond Angle ',10F5.0)
 6028 FORMAT (//,' Histogram of the ',I5,' Angles Analysed',/,
     +' Angle    0   40   50   60   70   80   90  100  110  120  130',
     +       '  140  150  180',/,' Number  ',13I5)
 6030 FORMAT (//,' A Total of ',I4,' Angles were outside range 70 to ',
     +  '150 degrees',/,' These are listed below as RES1,Water,R',
     +       'ES2 and Angle',/ (1X,4 (4I4,2X)))
C
      END
C
C
C     ==========================================
      SUBROUTINE SYMTRN2(NSYM,RSM,SYMCHS)
C     ==========================================
C
C---- This translates the Symmetry matrices into INT TAB
C     character strings for each symmetry operation and prints the real
C     and reciprocal space operators on standard output.
C
C     It gives the real and reciprocal space operations.
C   eg     X,Y,Z    H,K,L
C   eg     -Y,X-Y, Z   -H-K, H, L  etc
C   That is more complicated than you might think!!
C
C---- Inverse symmetry needed to test systematic absences -
C     copy Rsmm Rsmtt this common block.
C
C      COMMON /SYSabs/ NsymT,RsmM(4,4,MaxSym),RsmTT(4,4,MaxSym)
C
C
C
C     .. Parameters ..
      INTEGER MaxSymmetry
      PARAMETER (MaxSymmetry=230)
C     ..
C     .. Scalar Arguments ..
      INTEGER NSYM
C     ..
C     .. Array Arguments ..
      REAL RSM(4,4,*)
      CHARACTER SYMCHS(*)*80
C     ..
C     .. Local Scalars ..
      INTEGER I,ICH,IST,J,NS,TEMPSYM
C     ..
C     .. Local Arrays ..
      REAL RSMT(4,4,MaxSymmetry)
      CHARACTER HKLCR(3)*1
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR,RBRINV,SYMTR4
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC LEN
C     ..
C     .. Data statements ..
C
      DATA HKLCR/'H','K','L'/
C     ..
C
C
      TEMPSYM = NSYM
      CALL SYMTR4(NSYM,RSM,SYMCHS)
      DO 30 NS = 1,NSYM
C
C---- H K L   - get inverse symmetry operation
C
        CALL RBRINV(RSM(1,1,NS),RSMT(1,1,NS))
        ICH = 40
        DO 20 J = 1,3
          IST = 0
          DO 10 I = 1,3
            IF (RSMT(I,J,NS) .ne. 0) THEN
              IF (RSMT(I,J,NS) .gt. 0 .and. IST .gt. 0) THEN
                IF (ICH .gt. LEN(SYMCHS(1))) CALL CCPERR(1,
     +              'SYMTRN: character array too short')
                SYMCHS(NS) (ICH:ICH) = '+'
                ICH = ICH + 1
              END IF
C
C
              IF (RSMT(I,J,NS) .lt. 0) THEN
                IF (ICH .gt. LEN(SYMCHS(1))) CALL CCPERR(1,
     +              'SYMTRN: character array too short')
                SYMCHS(NS) (ICH:ICH) = '-'
                IST = 1
                ICH = ICH + 1
              END IF
C
C
              IF (ICH .gt. LEN(SYMCHS(1))) CALL CCPERR(1,
     +            'SYMTRN: character array too short')
              SYMCHS(NS) (ICH:ICH) = HKLCR(I)
              ICH = ICH + 1
              IST = 1
            END IF
   10     CONTINUE
C
C---- ADD COMMA space
C
          IF (J .ne. 3) THEN
            IF (ICH+2 .gt. LEN(SYMCHS(1))) CALL CCPERR(1,
     +          'SYMTRN: character array too short')
            SYMCHS(NS) (ICH:ICH+2) = ',  '
            ICH = ICH + 3
          END IF
   20   CONTINUE
   30 CONTINUE
C
C
      END
C
C     =========================
      SUBROUTINE sorthits(NhitS,NAhits,NBhits,NChits)
C     =========================
c
c
C
C
C     .. Scalar Arguments ..
      INTEGER NhitS
C     ..
C     .. Array Arguments ..
      INTEGER NAhits(NhitS),NBhits(NhitS),NChits(NhitS)
C     ..
C     .. Local Scalars ..
      INTEGER I,IH,IP1,J,K,NA,NCC,NCX,ND,NO2,NP1
C     ..
      NO2 = NhitS/2
c
c
      DO 40 I = 1,NO2
        NCC = NAhits(I)
        ND = NCC
        IH = I
        K = I
        IP1 = I + 1
        NP1 = NhitS + 1 - I
C
C
        DO 20 J = IP1,NP1
          NA = NAhits(J)
          IF (NA .ge. NCC) GO TO 10
          IH = J
          NCC = NA
          GO TO 20
   10     IF (ND .ge. NA) GO TO 20
          ND = NA
          K = J
   20   CONTINUE
C
C
        IF (NCC .eq. ND) GO TO 50
        IF (IH .eq. I) GO TO 30
        NCX = NBhits(IH)
        NBhits(IH) = NBhits(I)
        NBhits(I) = NCX

        NCC = NAhits(IH)
        NAhits(IH) = NAhits(I)
        NAhits(I) = NCC

        NCC = NChits(IH)
        NChits(IH) = NChits(I)
        NChits(I) = NCC

   30   IF (K .eq. I) K = IH
        IF (K .eq. NP1) GO TO 40

        NCX = NBhits(K)
        NBhits(K) = NBhits(NP1)
        NBhits(NP1) = NCX

        ND = NAhits(K)
        NAhits(K) = NAhits(NP1)
        NAhits(NP1) = ND

        ND = NChits(K)
        NChits(K) = NChits(NP1)
        NChits(NP1) = ND
   40 CONTINUE
C
C
   50 RETURN
      END


C
C     ============================================
      SUBROUTINE getspecial(XF,Lspecial,Lspecdone,cutoff)
C     ============================================
C
C---- This subroutine finds what coordinates occupy special positions
C     ie have occupancies less than 1.0
C     from consideration of the Symmetry Operations.
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains,MaxSpecial
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30,
     +           MaxSpecial=2000)
C     ..
C     .. Array Arguments ..
      REAL XF(3)
C     ..
C     .. Scalar Arguments ..
      REAL cutoff
      LOGICAL Lspecdone,Lspecial
C     ..
C     .. Local Scalars ..
      REAL D,XXF,YYF,ZZF
      INTEGER J,Jdo,Kdo,Nspec,NumSpecial
C     ..
C     .. Scalars in Common ..
      REAL FRdist
      INTEGER Icounts,Idents,NumSymm
C     ..
C     .. Arrays in Common ..
      REAL FRclose,RFin,Rsymm,S,XYZS
      INTEGER*2 InumAtmSymm,Nsymop,NUS
C     ..
C     .. Local Arrays ..
      REAL delta(3),XYZspecial(3,MaxSpecial),Rdiff
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPERR
C     ..
C     .. Common blocks ..
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
C     ..
C     .. Save statement ..
      SAVE
C     ..
C     .. Data statements ..
      DATA NumSpecial/0/
C     ..
      Lspecdone = .false.
      Lspecial = .false.
      Nspec = 0
C
C
C---- Generate symm equivs
C
C---- test whether xxf yyf zzf equals xf yf zf  - possibly with a unit
C---- cell translation...
C
      IF (Numsymm.le.1) then
          LLw= 1
          LLq=1
        else
          LLq=2
          llw=Numsymm + 1
        end if

      DO 10 J = llq,llw
        XXF = Rsymm(1,1,J)*XF(1) + Rsymm(1,2,J)*XF(2) +
     +        Rsymm(1,3,J)*XF(3) + Rsymm(1,4,J)
C
C
        IF (XXF .lt. XF(1)+cutoff .and. 
     +      XXF .gt. XF(1)-cutoff) THEN
          YYF = Rsymm(2,1,J)*XF(1) + Rsymm(2,2,J)*XF(2) +
     +          Rsymm(2,3,J)*XF(3) + Rsymm(2,4,J)
C
C
          IF (YYF .lt. XF(2)+cutoff .and. 
     +        YYF .gt. XF(2)-cutoff) THEN
            ZZF = Rsymm(3,1,J)*XF(1) + Rsymm(3,2,J)*XF(2) +
     +            Rsymm(3,3,J)*XF(3) + Rsymm(3,4,J)
C
C
            IF (ZZF .lt. XF(3)+cutoff .and. 
     +          ZZF .gt. XF(3)-cutoff) THEN
               Rdiff = (XXF - XF(1))**2 + 
     +                 (YYF - XF(2))**2 + 
     +                  (ZZF - XF(3))**2
            IF (Rdiff .lt. FRdist) Nspec = Nspec + 1
            END IF
          END IF
        END IF
   10 CONTINUE
C
C
      IF (Nspec .ne. 0) Lspecial = .true.
C
C
      IF (Lspecial) THEN
        NumSpecial = NumSpecial + 1
        IF (NumSpecial .gt. MaxSpecial) THEN
c            WRITE (6,6000) 
c 6000       FORMAT(' WARNING:: MaxSpecialatom exceeded')
            NumSpecial = NumSpecial - 1
        ELSE
        DO 20 Jdo = 1,3
          XYZspecial(Jdo,NumSpecial) = XF(Jdo)
   20   CONTINUE
        END IF
        IF (NumSpecial .gt. 1) THEN
          DO 40 Jdo = 1,NumSpecial - 1
            DO 30 Kdo = 1,3
              delta(Kdo) = XYZspecial(Kdo,Jdo) - XF(Kdo)
   30       CONTINUE
            D = delta(1)*delta(1) + delta(2)*delta(2) +
     +          delta(3)*delta(3)
            IF (D .le. 0.01) THEN
              Lspecdone = .true.
              RETURN
            END IF
   40     CONTINUE
        END IF
      END IF
C
      END
C
C     ========================================================
      SUBROUTINE gettooclose(XF,Lclose)
C     ========================================================
C
C     .. Parameters ..
      INTEGER MaxAtoms,MaxResidues,MaxSymmetry,MaxSize,MaxMetals,
     +        MaxSelectedAtoms,MaxSelectedResidues,
     +        MaxSelectedChains
      PARAMETER (MaxSelectedAtoms=25,MaxSelectedResidues=1998,
     +           MaxAtoms=48000,MaxResidues=9000,MaxSymmetry=230,
     +           MaxSize=8000000,MaxMetals=20,MaxSelectedChains=30)
C     ..
C     .. Scalar Arguments ..
      LOGICAL Lclose
C     ..
C     .. Array Arguments ..
      REAL XF(3)
C     ..
C     .. Scalars in Common ..
      REAL FRdist
      INTEGER Icounts,Idents,NumSymm
C     ..
C     .. Arrays in Common ..
      REAL FRclose,RFin,Rsymm,S,XYZS
      INTEGER*2 InumAtmSymm,Nsymop,NUS
C     ..
C     .. Local Scalars ..
      REAL Rdiff,XXF,YYF,ZZF
      INTEGER J,Nclose
C     ..
C     .. Common blocks ..
      COMMON /SymmCom/NumSymm,Idents,S(4,4,MaxSymmetry),NUS(2*MaxAtoms),
     +       InumAtmSymm(2*MaxAtoms),XYZS(3,2*MaxAtoms),
     +       Nsymop(2*MaxAtoms),Icounts,Rsymm(4,4,MaxSymmetry),
     +       RFin(4,4),FRclose(3),FRdist
C     ..
      SAVE
C
C
      Lclose = .false.
      Nclose = 0
      IF (Numsymm.le.1) then
          LLw= 1
           llq=1
        else
          llq=2
          llw=Numsymm + 1
        end if

      DO 10 J = llq,llw
        XXF = Rsymm(1,1,J)*XF(1) + Rsymm(1,2,J)*XF(2) +
     +        Rsymm(1,3,J)*XF(3) + Rsymm(1,4,J)
C
C
        IF (XXF .lt. XF(1)+FRclose(1) .and. 
     +       XXF .gt. XF(1)-FRclose(1)) THEN
          YYF = Rsymm(2,1,J)*XF(1) + Rsymm(2,2,J)*XF(2) +
     +          Rsymm(2,3,J)*XF(3) + Rsymm(2,4,J)
          IF (YYF .lt. XF(2)+FRclose(2) .and.
     +        YYF .gt. XF(2)-FRclose(2)) THEN
            ZZF = Rsymm(3,1,J)*XF(1) + Rsymm(3,2,J)*XF(2) +
     +            Rsymm(3,3,J)*XF(3) + Rsymm(3,4,J)
            IF (ZZF .lt. XF(3)+FRclose(3) .and.
     +          ZZF .gt. XF(3)-FRclose(3)) THEN
              Rdiff = (XXF-XF(1))**2 + (YYF-XF(2))**2 + (ZZF-XF(3))**2
              IF (Rdiff .lt. FRdist) Nclose = Nclose + 1
            END IF
          END IF
        END IF
   10 CONTINUE
C
C
      IF (Nclose .ne. 0) Lclose = .true.
      RETURN
      END
C
C     =======================================
      SUBROUTINE GETRecipcell(cell,Recipcell)
C     =======================================
C
C  Calculation of reciprocal cell parameters
C
C
C
C     .. Array Arguments ..
      REAL cell(6),Recipcell(6)
C     ..
C     .. Local Scalars ..
      REAL ALPHA,AR,BETA,BR,CONST,COSA,COSAS,COSB,COSBS,COSG,COSGS,CR,
     +     GAMMA,S,SINA,SINB,SINBS,SING,VOL
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC COS,SIN,sqrt
C     ..
C
C
      CONST = 3.1415926/180.0
      ALPHA = cell(4)*CONST
      BETA = cell(5)*CONST
      GAMMA = cell(6)*CONST
      COSA = COS(ALPHA)
      SINA = SIN(ALPHA)
      COSB = COS(BETA)
      SINB = SIN(BETA)
      COSG = COS(GAMMA)
      SING = SIN(GAMMA)
C
      IF (SINA.EQ.0 .OR. SINB.EQ.0 .OR. SING.EQ.0) THEN
        WRITE (6,FMT=1010)
 1010   FORMAT (/,' * Warning: One or more cell angles have been set',
     +  ' to zero *',//,
     +  ' This may be because XYZIN does not contain CRYST1 or',
     +  ' SCALE cards,',/,
     +  ' which are required for MODEs AUTO or IMOL.',/,
     +  ' Check that the input file contains the correct cell',
     +  ' information.',/)
        CALL CCPERR (1, 'GETRecipcell: bad or missing cell angle(s)')
      END IF
C
      COSAS = (COSB*COSG-COSA)/ (SINB*SING)
      COSBS = (COSA*COSG-COSB)/ (SINA*SING)
      COSGS = (COSA*COSB-COSG)/ (SINA*SINB)
      SINBS = sqrt(1.0-COSBS*COSBS)
      AR = 1.0/ (cell(1)*SINB*sqrt(1.0-COSGS*COSGS))
      BR = 1.0/ (cell(2)*SING*sqrt(1.0-COSAS*COSAS))
      CR = 1.0/ (cell(3)*SINA*SINBS)
      Recipcell(1) = AR
      Recipcell(2) = BR
      Recipcell(3) = CR
      Recipcell(4) = COSAS
      Recipcell(5) = COSBS
      Recipcell(6) = COSGS
C
C----  cell volume
C
      S = (ALPHA+BETA+GAMMA)*0.5
      VOL = cell(1)*2.0*cell(2)*cell(3)*
     +      sqrt(SIN(S-ALPHA)*SIN(S)*SIN(S-BETA)*SIN(S-GAMMA))
C
C
      END
C
C     ======================================
      SUBROUTINE FRCLIM(Recipcell,DELD,FDEL,FRdist)
C     ======================================
C
C Set fractional limits FDEL equivalent to distance DELD
C   = DELD . a*, .b*, .c*
C
C
C
C     .. Scalar Arguments ..
      REAL DELD,FRdist
C     ..
C     .. Array Arguments ..
      REAL FDEL(3),Recipcell(6)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
C
C
      DO 10 I = 1,3
        FDEL(I) = Recipcell(I)*DELD
   10 CONTINUE
C
      FRdist = FDEL(1)**2 + FDEL(2)**2 + FDEL(3)**3
C
      END
C
C
C     ========================
      SUBROUTINE RELOAD(ARRAY)
C     ========================
C
C     Reloads the 125-element arrays so that the first 27 elements
C     contain the information about +/-1 lattice translations.
C
C     It is needed to get the right titles etc for the different
C     translations.
C
      IMPLICIT NONE
C
C     ..Array Arguments
      CHARACTER*(*) ARRAY(125)
C
C     ..Local scalars
      INTEGER I,J,K,INDX1,INDX2
C
      DO 10 I=1,3
        DO 20 J=1,3
          DO 30 K=1,3
C
            INDX1 = (I-1)*9 + (J-1)*3 + K
            INDX2 = ((I*5)+J)*5 + K + 1
C
            ARRAY(INDX1) = ARRAY(INDX2)
C
 30       CONTINUE
 20     CONTINUE
 10   CONTINUE
C
      RETURN
      END
C
C