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
      PROGRAM VECREF
C
C====   VECREF Vector space refinement of heavy atom parameters using SIR
C       difference Patterson.
C
C====   Ian Tickle
C
C       Phil Evans                MRC LMB                   13-Nov-1989
C       Changed input to read lattice type as character.
C       Occupancies are now true occupancies for centred lattices,
C       form-factors are multiplied by lattice multiplicity.
C
C       Ian Tickle   24-Jan-1990
C       Improvements to group form-factor option; print curve for
C       use in PHASE; print vector peak profile.
C       New peak radius algorithm, old one had problems with group
C       form-factors and B's > 45.  Radius calculations use same
C       function RADIUS.
C       PKLIST routine moved to more logical place; fixed bug in
C       asymmetric unit check.
C       Residuals & correlation coefficients computed for each atom.
C
C       Ian Tickle   26-Apr-1990
C       Correct bug in atom statistics.
C       Make atom labels CHARACTER*4.
C       Change MGULP to MGULPR.
C
C       Ian Tickle   8-Jan-1991
C       Fixed bug in Harker vector weights in trigonal and hexagonal.
C
C       Ian Tickle   19-May-1992
C       Added keywording.
C
C       Ian Tickle   1-Nov-1992
C       Added low resolution cutoff correction.
C
C       Ian Tickle   27-Jun-1994
C       Added RCUT & DAMP options.  Improved radius calculation.
C       Used ATOMSF.LIB instead of SCATCOEF.LIB
C
C       Ian Tickle   24-Aug-1994
C       Reject atoms with B <= 0.
C
C       Ian Tickle   29-Sep-1994
C       Allow negative occupancies.  Added THRE & BLIM options.
C
C       Ian Tickle   11-Oct-1994
C       Added dynamic memory allocation.
C
C       Ian Tickle   11-Dec-1995
C       Added shift factor for coordinates to fix problem of instability in
C       high symmetry space groups.
C
C
        IMPLICIT NONE
        INTEGER MPTD,MAT,MCF,MST
        PARAMETER (MPTD=200000,MAT=50,MCF=10,MST=8)
C
        CHARACTER EVAL*80,TARR1(3),TARR2(1)
        INTEGER I,MAXPTS,LARR1(3),LARR2(1)
C
        INTEGER MPT,NPT,NCYC(3),IVF(4)
        COMMON /REFCOM/ MPT,NPT,NCYC,IVF
C
        INTEGER IUVW(3),NUVW(3),NXYZ(3)
        INTEGER LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        REAL VOL,OBF,RCG(3),TM(3,3)
        COMMON /CELCOM/ VOL,OBF,RCG,TM
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
C       Block data: to make sure it's linked in:
        EXTERNAL VR
        EXTERNAL GENREF,PKLIST
C
        DATA TARR1,TARR2 /'I','I','R','R'/
        DATA NCYC /3*0/
C
        CALL CCPRCS(6,'VECREF','$Date: 1997/10/14 14:39:35 $')
        WRITE (6,'(/3(/20X,A)//)')
     &  ' ###############################################',
     &  ' # VECTOR SPACE SINGLE ISOMORPHOUS REFINEMENT. #',
     &  ' ###############################################'
        CALL CCPFYP
C
C====   Read all input data.
C
        CALL GETINP
C
C====   Set up lengths of dynamic arrays.
C
        CALL UGTENV('VECREF_MAXPTS',EVAL)
        MAXPTS=MPTD
        IF (EVAL.NE.' ') READ (EVAL,*,ERR=80) MAXPTS
        IF (MAXPTS.LE.0) GOTO 80
        MPT=MAXPTS
        LARR1(1)=(9*MPT+1)/2
        LARR1(2)=2*MPT
        LARR1(3)=NUVW(1)*NUVW(2)
        LARR2(1)=NUVW(1)*NUVW(2)*NUVW(3)
C
C====   Generate list of points and do refinement cycles on occupancies first.
C
C        OPEN (2,FILE='POINTS',STATUS='UNKNOWN')
        IVF(1)=0
        IVF(2)=0
        IVF(3)=1
        IVF(4)=0
C
10      CALL CCPALC(GENREF,3,TARR1,LARR1)
        IF (NPT.GE.MPT) THEN
          MPT=NPT+1
          LARR1(1)=(9*MPT+1)/2
          LARR1(2)=2*MPT
          GOTO 10
        ENDIF
C
        IF (NCYC(1).LT.ABS(NCYC(3))) THEN
          CALL CCPALC(PKLIST,1,TARR2,LARR2)
C
C====   Now set coordinates to be refined.
C
          IVF(1)=1
          DO 30 I=NCYC(1)+1,NCYC(2)
            WRITE (6,'(///A,I3)') ' Cycle',I
20          CALL CCPALC(GENREF,3,TARR1,LARR1)
            IF (NPT.GE.MPT) THEN
              MPT=NPT+1
              LARR1(1)=(9*MPT+1)/2
              LARR1(2)=2*MPT
              GOTO 20
            ENDIF
30        CONTINUE
C
C====   Now set all parameters to be refined.
C
          IF (NCYC(3).LT.0) THEN
            IVF(2)=1
          ELSE
            IVF(4)=1
          ENDIF
C
          DO 50 I=NCYC(2)+1,ABS(NCYC(3))
            WRITE (6,'(///,A,I3)') ' Cycle',I
40          CALL CCPALC(GENREF,3,TARR1,LARR1)
            IF (NPT.GE.MPT) THEN
              MPT=NPT+1
              LARR1(1)=(9*MPT+1)/2
              LARR1(2)=2*MPT
              GOTO 40
            ENDIF
50        CONTINUE
        ENDIF
C        CLOSE (2)
C
        IF (NCYC(3).NE.0) CALL WRITAT
C
C====   Reset flags and do a final cycle for statistics.
C
        DO 60 I=1,4
60        IVF(I)=0
C
70      CALL CCPALC(GENREF,3,TARR1,LARR1)
        IF (NPT.GE.MPT) THEN
          MPT=NPT+1
          LARR1(1)=(9*MPT+1)/2
          LARR1(2)=2*MPT
          GOTO 70
        ENDIF
C
        CALL CCPALC(PKLIST,1,TARR2,LARR2)
C
        IF (MPT.GT.MAXPTS) WRITE (6,'(/2(/1X,A),I9/)')
     &  '                    ********** NOTE **********',
     &  'To save time in future runs, set VECREF_MAXPTS to at least',MPT
        CALL CCPERR(0,'**** NORMAL TERMINATION ****')
80      CALL CCPERR(1,'*** ERROR in VECREF: Bad VECREF_MAXPTS')
        END
C
C
        BLOCK DATA VR
        PARAMETER (MAT=50,MCF=10,MST=8)
        INTEGER IUVW(3),NUVW(3),NXYZ(3),LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        DATA SXYZ /3*1./
        DATA IOR /.TRUE./
        DATA ITC /
     &  12*0,
     &  3*0, 3*6, 6*0,
     &  3*0, -4,4,4, 4,-4,-4, 3*0,
     &  3*0, 0,6,6, 6,0,6, 6,6,0,
     &  3*0, 0,6,6, 6*0,
     &  3*0, 6,0,6, 6*0,
     &  3*0, 6,6,0, 6*0/
        DATA NAT /0/
        DATA DMIN,DMAX,RCUT /3*0./
        DATA DAMP /.25/
        DATA TRMS /2./
        DATA BLIM /0., 200./
        DATA F0 /MCF*0./
        END
C
C
        SUBROUTINE GENREF(L1,IPTS,L2,ISORT,L3,SECT)
        IMPLICIT NONE
        INTEGER I,L1,L2,L3
C
        INTEGER MPT,NPT,NCYC(3),IVF(4)
        COMMON /REFCOM/ MPT,NPT,NCYC,IVF
        INTEGER*2 IPTS(3,3,MPT)
        INTEGER ISORT(2,MPT)
C
        INTEGER IUVW(3),NUVW(3),NXYZ(3),LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
        REAL SECT(LU:MU,LV:MV)
C
        CALL GENPTS(IPTS,ISORT)
        IF (NPT.GE.MPT) RETURN
C
        IF (IVF(1).EQ.0 .AND. IVF(3).EQ.1) THEN
          DO 10 I=1,NCYC(1)
            WRITE (6,'(///A,I3)') ' Cycle',I
10          CALL REFCYC(IPTS,ISORT,SECT)
        ELSE
          CALL REFCYC(IPTS,ISORT,SECT)
        ENDIF
        END
C
C
        SUBROUTINE GETINP
C
        IMPLICIT NONE
        INTEGER MAT,MCF,MST,MTOK,NKEY
        PARAMETER (MAT=50,MCF=10,MST=8,MTOK=20,NKEY=15)
        REAL PI,PI4,DTR
        PARAMETER (PI=3.141592654,PI4=4.*PI,DTR=PI/180.)
C
        CHARACTER LINE*400,KEY*4,CVAL(MTOK)*4,KEYA(NKEY)*4,TITLE*100,
     &  SGN*8,LATTIC*7,GS*4,AS1*4,AS2*4,AST(MCF)*4
        LOGICAL EOF
        INTEGER IFAIL,KE,MODE,I,J,IB,NSG,IT,N,IE,IS,NCF,IC1,IC2,IA,IAG,
     &  IC,NTOK
        REAL RL,RM,RS,RR,S,F,W1,R1,W2,R2,PI4R1,PI4R2,R,G1,A,G2,F1,F2
        INTEGER LENSTR,MULTIP
        INTEGER NTOKS(NKEY),IBEG(MTOK),IEND(MTOK),ITYP(MTOK),IDEC(MTOK),
     &  NXYZH(3),IXH(3),IVH(3)
        REAL FVAL(MTOK),CELL(6),CG(3),CA(3),SA(3),RTS(4,4,96),
     &  ACF(4,MCF),BCF(4,MCF),CCF(MCF),FT(0:29)
C
        INTEGER MPT,NPT,NCYC(3),IVF(4)
        COMMON /REFCOM/ MPT,NPT,NCYC,IVF
C
        INTEGER IUVW(3),NUVW(3),NXYZ(3)
        INTEGER LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
C
        REAL VOL,OBF,RCG(3),TM(3,3)
        COMMON /CELCOM/ VOL,OBF,RCG,TM
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        DATA KEYA/
     &  'TITL','SPAC','RESO','RCUT','GROU','CYCL','ORIG','SCAL','ATOM',
     &  'ANIS','BREF','DAMP','THRE','BLIM','END '/
        DATA NTOKS/
     &      1,    -2,     2,    -2,     2,     2,    -1,    -4,    -8,
     &     -9,     1,    -2,    -2,    -3,     1/
        DATA TITLE /' '/, SGN /' '/, GS /' '/
        DATA NSG /0/, NXYZH /3*20/, IXH /1,3,9/, NCF /0/
        DATA LATTIC /'PIRFABC'/
        DATA IB /0/
        DATA AS1 /' '/, W1 /1./, R1 /0./, AS2 /' '/, W2 /0./, R2 /0./
        DATA IAG /0/
C
        KE=0
C
C====   Open Patterson density map for reading.
C
C IJT   I dummied out space group number in next call because I realised
C       that for Pattersons FFT now writes the Laue group NOT the true
C       space group! (which I suppose makes some sense).  This meant that
C       if the user didn't specify SPACEGROUP he got wrong results!
C
        CALL MRDHDR(1,'MAPIN',LINE,NUVW(3),IUVW,NXYZ,LW,
     &  LU,MU,LV,MV,CELL,I,MODE,RL,RM,RS,RR)
        NUVW(1)=MU-LU+1
        NUVW(2)=MV-LV+1
        MW=LW+NUVW(3)-1
C
C====   Get cos & sines of cell angles.
C
        S=1.
        F=2.
        DO 10 I=1,3
          CG(I)=CELL(I)/NXYZ(I)
          IF (CELL(I+3).EQ.90.) THEN
            CA(I)=0.
            SA(I)=1.
          ELSE
            CA(I)=COS(DTR*CELL(I+3))
            SA(I)=SIN(DTR*CELL(I+3))
          ENDIF
          S=S-CA(I)**2
10        F=F*CA(I)
        F=SQRT(S+F)
C
C====   Get cell volume, reciprocal cell and metric tensor.
C
        VOL=F
        DO 20 I=1,3
          RCG(I)=SA(I)/(F*CG(I))
          VOL=VOL*CELL(I)
          TM(I,I)=CG(I)**2
          DO 20 J=1,3
20          IF (I.NE.J) TM(I,J)=CG(I)*CG(J)*CA(6-I-J)
CD      WRITE (6,'(/A,3F8.4,F10.1)') ' RCG, VOL',RCG,VOL
CD      WRITE (6,'(A,3(/1X,3F10.4))') ' TM',((TM(I,J),J=1,3),I=1,3)
C
C====   Read control data.
C
220     LINE=' '
        NTOK=MTOK
        CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,FVAL,CVAL,IDEC,NTOK,EOF,
     &              .FALSE.)
        IF (EOF) GOTO 1100
        IF (NTOK.EQ.0) GOTO 220
        WRITE (6,'(/2A)') ' >>>> Input line >>>>  ',LINE(:LENSTR(LINE))
        DO 230 I=1,NKEY
230       IF (KEY.EQ.KEYA(I)) GOTO 240
        WRITE (6,*) '*** ERROR: Last line has invalid keyword.'
        KE=KE+1
        GOTO 220
240     IF ((NTOKS(I).LT.0 .AND. NTOK.NE.-NTOKS(I)).OR.
     &  NTOK.LT.NTOKS(I)) THEN
          WRITE (6,'(/A)') ' **** ERROR: bad argument count.'
          KE = KE + 1
        ENDIF
        GOTO (1001,1002,1003,1004,1005,1006,1007,1008,1009,1010,1011,
     &  1012,1013,1014,1100),I
C
C====   Title.
C
1001    IF (IBEG(2).GT.0) TITLE=LINE(IBEG(2):)
        GOTO 220
C
C====   Space group.
C
1002    IF (ITYP(2).EQ.1) THEN
          SGN=LINE(IBEG(2):IEND(2))
          CALL CCPUPC(SGN)
          NSG=0
        ELSE
          CALL GTPINT(2,NSG,NTOK,ITYP,FVAL)
          SGN=' '
        ENDIF
        GOTO 220
C
C====   Resolution.
C
1003    CALL GTPREA(2,DMIN,NTOK,ITYP,FVAL)
        IF (NTOK.GT.2) THEN
          CALL GTPREA(3,DMAX,NTOK,ITYP,FVAL)
          IF (DMIN.GT.DMAX) THEN
            S=DMIN
            DMIN=DMAX
            DMAX=S
          ENDIF
        ENDIF
        GOTO 220
C
C====   Radius cutoff.
C
1004    CALL GTPREA(2,RCUT,NTOK,ITYP,FVAL)
        GOTO 220
C
C====   Group definition.
C
1005    GS=CVAL(2)
        CALL CCPUPC(GS)
        AS1=CVAL(3)
        CALL CCPUPC(AS1)
        CALL GTPREA(4,W1,NTOK,ITYP,FVAL)
        CALL GTPREA(5,R1,NTOK,ITYP,FVAL)
        IF (W1.EQ.0.) R1=0.
        AS2=CVAL(6)
        CALL CCPUPC(AS2)
        CALL GTPREA(7,W2,NTOK,ITYP,FVAL)
        CALL GTPREA(8,R2,NTOK,ITYP,FVAL)
        IF (W2.EQ.0.) R2=0.
        GOTO 220
C
C====   No of cycles.
C
1006    CALL GTPINT(2,NCYC(1),NTOK,ITYP,FVAL)
        CALL GTPINT(3,NCYC(2),NTOK,ITYP,FVAL)
        CALL GTPINT(4,NCYC(3),NTOK,ITYP,FVAL)
        GOTO 220
C
C====   Do not exclude origin.
C
1007    IOR=.FALSE.
        GOTO 220
C
C====   Coordinate scales.
C
1008    CALL GTPREA(2,SXYZ(1),NTOK,ITYP,FVAL)
        CALL GTPREA(3,SXYZ(2),NTOK,ITYP,FVAL)
        CALL GTPREA(4,SXYZ(3),NTOK,ITYP,FVAL)
        GOTO 220
C
C====   Atom record.
C
1009    DO 260 I=2,8
          IF ((I.EQ.2 .AND. ITYP(I).EQ.0) .OR.
     &        (I.GT.2 .AND. ITYP(I).NE.2)) THEN
            WRITE (6,*) '*** ERROR: Atom record format.'
            KE=KE+1
            GOTO 220
          ENDIF
260     CONTINUE
        NAT=NAT+1
        IF (NAT.LE.MAT) THEN
          AS(NAT)=LINE(IBEG(2):IEND(2))
          CALL CCPUPC(AS(NAT))
          CALL GTPINT(3,IAT(NAT),NTOK,ITYP,FVAL)
          CALL GTPREA(4,OCC(NAT),NTOK,ITYP,FVAL)
          CALL GTPREA(5,XAT(1,NAT),NTOK,ITYP,FVAL)
          CALL GTPREA(6,XAT(2,NAT),NTOK,ITYP,FVAL)
          CALL GTPREA(7,XAT(3,NAT),NTOK,ITYP,FVAL)
          CALL GTPREA(8,BAT(NAT),NTOK,ITYP,FVAL)
          IF (OCC(NAT).EQ.0. .OR. (TRMS.GT.0. .AND. OCC(NAT).LT.0.) .OR.
     &        BAT(NAT).LE.BLIM(1) .OR. BAT(NAT).GT.BLIM(2)) THEN
            WRITE (6,'(/2A,I4,A)') ' Atom  ',AS(NAT),IAT(NAT),
     &      '  has occupancy or B out of range.  *** REJECTED ***'
            NAT=NAT-1
          ENDIF
        ENDIF
        GOTO 220
C
C====   Anisotropic B factors.
C
1010    CONTINUE
        GOTO 220
C
C====   B factor refine flags.
C
1011    IB=1
        GOTO 220
C
C====   DAMPing factor..
C
1012    CALL GTPREA(2,DAMP,NTOK,ITYP,FVAL)
        GOTO 220
C
C====   THREshold RMS for occupancy rejection.
C
 1013   CALL GTPREA(2,TRMS,NTOK,ITYP,FVAL)
        GOTO 220
C
C====   B rejection limits.
 1014   CALL GTPREA(2,BLIM(1),NTOK,ITYP,FVAL)
        CALL GTPREA(3,BLIM(2),NTOK,ITYP,FVAL)
        GOTO 220
C
C====   End/EOF
C
 1100   WRITE (6,'(//1X,A)') TITLE
C
C====   Decode first character of space group as lattice type
C
        IF (SGN.EQ.' ') THEN
          IF (NSG.LE.0) THEN
            WRITE (6,*)
     &      ' *** ERROR: Invalid or missing value for space group.'
            KE=KE+1
            NSG=1
          ENDIF
        ENDIF
        CALL MSYMLB(1,NSG,SGN,KEY,NGEP,I,RTS)
        LTY=INDEX(LATTIC,SGN(:1))
C
        NTC=LTY
        IF (NTC.GT.4) NTC=2
        WRITE (6,'(/3A,I2,A,I4/)')
     &  ' Space group = ',SGN,'    Multiplicity =',NTC,
     &  '    Number of g.e.p.s =',NGEP
        IF (NGEP.GT.24) THEN
          WRITE (6,*) '*** ERROR: Invalid number of geps.'
          KE=KE+1
        ENDIF
        DO 30 IT=1,NTC
          DO 30 I=1,3
            N=ITC(I,IT,LTY)*NXYZ(I)
            IF (MOD(N,12).NE.0) THEN
              WRITE (6,*) '*** ERROR: Centring translation not integer.'
              KE=KE+1
            ENDIF
30          ITC(I,IT,LTY)=N/12
CD      WRITE (6,'(A,3I6)') (' TC',(ITC(I,IT,LTY),I=1,3),IT=1,NTC)
C
        DO 40 IE=1,NGEP
          DO 40 I=1,3
            RTE(I,4,IE)=RTS(I,4,IE)*NXYZ(I)
            DO 40 J=1,3
40            RTE(I,J,IE)=RTS(I,J,IE)
C
C====   Find identity position and compute multiplicities of Harker vectors.
C
        IDEN=0
        DO 70 IE=1,NGEP
          DO 60 I=1,3
            IVH(I)=IXH(I)
            DO 60 J=1,3
60            IVH(I)=IVH(I)-NINT(RTE(I,J,IE))*IXH(J)
          IF (IVH(1).EQ.0 .AND. IVH(2).EQ.0 .AND. IVH(3).EQ.0) THEN
            IF (IDEN.NE.0) KE=KE+1
            IDEN=IE
            WH(IE)=.5
          ELSE
            WH(IE)=1./MULTIP(NXYZH,NGEP,RTE,IVH)
          ENDIF
70      CONTINUE
C
        IF (IDEN.EQ.0) THEN
          WRITE (6,*) '*** ERROR: Identity position not found.'
          KE=KE+1
        ENDIF
        WRITE (6,'(/A,3(/1X,8F10.4))')
     &  ' Harker vector weights :',(WH(IE),IE=1,NGEP)
C
C====   Integration intervals.
C
        IF (DMIN.LE.0.) THEN
          WRITE (6,*) '*** ERROR: Invalid or missing value for Dmin.'
          KE=KE+1
          DMIN=3.
        ENDIF
        IF (DMAX.EQ.DMIN) THEN
          WRITE (6,*) '*** ERROR: Invalid value for Dmax.'
          KE=KE+1
          DMAX=0.
        ENDIF
C
        DO 270 I=1,3
          IF (CG(I).GT.DMIN/3.) THEN
            WRITE (6,'(/A)') ' *** ERROR: Map grid too coarse.'
            KE=KE+1
            GOTO 280
          ENDIF
270     CONTINUE
C
280     DSTMIN=0.
        IF (DMAX.GT.0.) DSTMIN=.5/DMAX
        DSTMAX=.5/DMIN
        STH=(DSTMAX-DSTMIN)/MST
        S=1./SQRT(3.)
        DO 80 IS=1,MST-1,2
          ST(IS)=DSTMIN+(IS-S)*STH
          SS(IS)=ST(IS)**2
          ST(IS+1)=DSTMIN+(IS+S)*STH
80        SS(IS+1)=ST(IS+1)**2
CD        WRITE (6,'(/A,8F9.4)') ' ST',ST
CD        WRITE (6,'(/A,8F9.5)') ' SS',SS
        STH=8.*STH
        DSTMIN=PI4*DSTMIN
        DSTMAX=PI4*DSTMAX
        OBF=1./DMIN**3
        IF (DMAX.GT.0.) OBF=OBF-1./DMAX**3
        OBF=PI4*VOL*OBF/(3.*(NXYZ(1)*NXYZ(2)*NXYZ(3)))
C
        NCYC(1)=MIN(NCYC(1),10)
        NCYC(2)=MIN(NCYC(2),10)
        WRITE (6,'(3(/1X,A,F8.2),5(/1X,A,I4))')
     &  '                                  Dmax =',DMAX,
     &  '                                  Dmin =',DMIN,
     &  '                                  Rcut =',RCUT,
     &  'Number of cycles to refine occupancies =',NCYC(1),
     &  'Number of cycles to refine coordinates =',NCYC(2),
     &  'Number of cycles to refine all params  =',NCYC(3)
        NCYC(2)=NCYC(1)+NCYC(2)
        NCYC(3)=NCYC(2)+NCYC(3)
C
        IF (IOR) THEN
          WRITE (6,'(/A)')
     &    ' Patterson origin to be excluded in refinement.'
        ELSE
          WRITE (6,'(/A)')
     &    ' Patterson origin to be included in refinement.'
        ENDIF
C
        IF (NCYC(3).GT.NCYC(2)) THEN
          IF (IB.EQ.0) THEN
            WRITE (6,'(/A)')
     &      ' Overall temperature factor to be refined.'
          ELSE
            WRITE (6,'(/A)')
     &      ' Individual atomic temperature factors to be refined.'
            NCYC(3)=-NCYC(3)
          ENDIF
          WRITE (6,'(/1X,A,F8.2)')
     &    '                                  DAMP =',DAMP
        ENDIF
        WRITE (6,'(2(/1X,A,F8.2),F8.2)')
     &  '                   Occupancy threshold =',TRMS,
     &  '                    B rejection limits =',BLIM
C
        IF (GS.NE.' ') THEN
          WRITE (6,'(/2A)') ' Group ',GS
          CALL ADDAST(AS1,AST,MCF,NCF,IC1,KE)
          WRITE (6,'(2A,4X,A,F8.1,4X,A,F8.2)')
     &    ' Atomic symbol ',AS1,'Multiplicity =',W1,'Radius =',R1
          PI4R1=PI4*R1
          IF (AS2.NE.' ') THEN
            CALL ADDAST(AS2,AST,MCF,NCF,IC2,KE)
            WRITE (6,'(2A,4X,A,F8.1,4X,A,F8.2)')
     &      ' Atomic symbol ',AS2,'Multiplicity =',W2,'Radius =',R2
            PI4R2=PI4*R2
          ENDIF
        ENDIF
C
        WRITE (6,'(/A,3F10.2)')
     &  ' Grid intervals for atomic coordinates =',SXYZ
        FRACT=.TRUE.
        DO 90 I=1,3
          IF (ABS(SXYZ(I)).GE.10.) FRACT=.FALSE.
90        SXYZ(I)=SXYZ(I)/NXYZ(I)
C
        IF (NAT.LE.0 .OR. NAT.GT.MAT) THEN
          WRITE (6,*) '*** ERROR: Invalid number of atoms.'
          NAT=MIN(NAT,MAT)
          KE=KE+1
        ENDIF
C
        DO 110 IA=1,NAT
          CALL ADDAST(AS(IA),AST,MCF,NCF,ICF(IA),KE)
          IF (FRACT) THEN
            WRITE (6,'(/A,I4,8X,A,I4,4F8.4,F8.2)') ' Atom',ICF(IA),
     &      AS(IA),IAT(IA),OCC(IA),(XAT(I,IA),I=1,3),BAT(IA)
          ELSE
            WRITE (6,'(/A,I4,8X,A,I4,F8.4,4F8.2)') ' Atom',ICF(IA),
     &      AS(IA),IAT(IA),OCC(IA),(XAT(I,IA),I=1,3),BAT(IA)
          ENDIF
          DO 100 I=1,3
            XAT(I,IA)=XAT(I,IA)/SXYZ(I)
100         XATS(I,IA)=XAT(I,IA)
          BATS(IA)=BAT(IA)
          OCCS(IA)=OCC(IA)
110       IF (IAG.EQ.0 .AND. AS(IA).EQ.GS) IAG=IA
        WRITE (6,'(/A,I5)') ' Number of atoms in asymmetric unit =',NAT
C
C====   Form-factors.
C
        IFAIL=0
        CALL CCPDPN(1,'ATOMSF','READONLY','F',0,IFAIL)
C
C BUG: Form factors should be fudged by sqrt(multiplicity) because there
C      are N times as many vectors, not vectors with N times the peak height.
C
        F=SQRT(REAL(NTC))
120     READ (1,'(A)',END=170) KEY
        CALL CCPUPC(KEY)
        DO 130 IC=1,NCF
130       IF (KEY.EQ.AST(IC)) GOTO 140
        GOTO 120
C
140     READ (1,'(22X,F14.6)') CCF(IC)
        READ (1,'(4(2X,F14.6))') (ACF(I,IC),I=1,4)
        READ (1,'(4(2X,F14.6))') (BCF(I,IC),I=1,4)
        READ (1,'(2X,F14.6)') F1
        DO 145 I=1,4
145       ACF(I,IC)=F*ACF(I,IC)
        CCF(IC)=F*(CCF(IC)+F1)
C       WRITE (6,'(1X,A,5F14.6/5X,4F14.6)')
C     & AST(IC),(ACF(I,IC),I=1,4),CCF(IC),(BCF(I,IC),I=1,4)
C
        F0(IC)=CCF(IC)
        DO 150 I=1,4
150       F0(IC)=F0(IC)+ACF(I,IC)
        F0(IC)=NGEP*F0(IC)
        DO 160 IS=1,MST
          FF(IS,IC)=CCF(IC)
          DO 160 I=1,4
160         FF(IS,IC)=FF(IS,IC)+ACF(I,IC)*EXP(-SS(IS)*BCF(I,IC))
C        WRITE (6,'(2A,F10.5/1X,8F10.5)')
C     &  ' FF ',AST(IC),F0(IC),(FF(IS,IC),IS=1,MST)
        GOTO 120
170     CLOSE(1)
C
        R=.5*SQRT(.5)*DMIN
        DO 210 IC=1,NCF
          R0(IC)=R
          IF (AST(IC).EQ.GS) THEN
            R0(IC)=MAX(R,SQRT(.5)*MAX(ABS(R1),ABS(R2)))
            F0(IC)=W1*F0(IC1)+W2*F0(IC2)
            DO 180 IS=1,MST
              G1=W1
              A=PI4R1*ST(IS)
              IF (A.NE.0.) G1=G1*SIN(A)/A
              G2=W2
              A=PI4R2*ST(IS)
              IF (A.NE.0.) G2=G2*SIN(A)/A
180           FF(IS,IC)=G1*FF(IS,IC1)+G2*FF(IS,IC2)
CD        WRITE (6,'(2A,9F8.2)')
CD    &   ' FF ',AST(IC),F0(IC),(FF(IS,IC),IS=1,MST)
C
            N=MIN(NINT(25./DMIN)+2,29)
            DO 200 J=0,N
              S=.02*J
              A=S*S
              F1=CCF(IC1)
              F2=CCF(IC2)
              DO 190 I=1,4
                F1=F1+ACF(I,IC1)*EXP(-A*BCF(I,IC1))
190             F2=F2+ACF(I,IC2)*EXP(-A*BCF(I,IC2))
              G1=W1
              A=PI4R1*S
              IF (A.NE.0.) G1=G1*SIN(A)/A
              G2=W2
              A=PI4R2*S
              IF (A.NE.0.) G2=G2*SIN(A)/A
200           FT(J)=G1*F1+G2*F2
            WRITE (6,'(////A/4(/1X,5F10.2))')
     &      ' GROUP FORM-FACTOR CURVE AT '//
     &      'INTERVALS OF 0.02 SIN(THETA)/LAMBDA',(FT(J),J=0,N)
          ENDIF
C
          IF (F0(IC).EQ.0.) THEN
            WRITE (6,*)
     &      '*** ERROR: Atom symbol not in scat. fact. list = ',AST(IC)
            KE=KE+1
          ENDIF
210     CONTINUE
C
        IF (KE.GT.0) CALL CCPERR(1,
     &  '*** ERROR(S) IN GETINP: See preceding error message(s).')
        IF (IAG.GT.0) CALL CCP4PROFIL(IAG)
        END
C
C
        INTEGER FUNCTION MULTIP(NXYZ,NGEP,RTE,IX)
C
C====   Compute multiplicity of point.
C
        IMPLICIT NONE
        INTEGER NGEP
        INTEGER NXYZ(3),IX(3)
        REAL RTE(3,4,24)
        INTEGER I,IE,J,K,M
        INTEGER JX(3),KX(3)
C
        DO 10 I=1,3
          JX(I)=MOD(IX(I),NXYZ(I))
10        IF (JX(I).LT.0) JX(I)=JX(I)+NXYZ(I)
C
        M=0
        DO 60 IE=1,NGEP
          DO 30 I=1,3
            KX(I)=0
            DO 20 J=1,3
20            KX(I)=KX(I)+NINT(RTE(I,J,IE))*JX(J)
30          KX(I)=MOD(KX(I),NXYZ(I))
          DO 60 K=1,2
            DO 40 I=1,3
              IF (KX(I).LT.0) KX(I)=KX(I)+NXYZ(I)
40            IF (JX(I).NE.KX(I)) GOTO 50
            M=M+1
50          DO 60 I=1,3
60            KX(I)=-KX(I)
        MULTIP=M
        END
C
C
        SUBROUTINE ADDAST(AS,AST,MCF,NCF,IC,KE)
C
C====   Add atom (or group) symbol to table.
C
        IMPLICIT NONE
        CHARACTER AS*(*),AST(*)*(*)
        INTEGER MCF,NCF,IC,KE
C
        DO 10 IC=1,NCF
10        IF (AS.EQ.AST(IC)) RETURN
        IF (NCF.EQ.MCF) THEN
          WRITE (6,*) '*** ERROR IN ADDAST: array bound check (MCF)'
          KE=KE+1
          RETURN
        ENDIF
        NCF=NCF+1
        AST(NCF)=AS
        IC=NCF
        END
C
C
        SUBROUTINE CCP4PROFIL(IA)
C
        IMPLICIT NONE
        INTEGER MPC,MAT,MCF,MST
        PARAMETER (MPC=120,MAT=50,MCF=10,MST=8)
        REAL PI,PI4
        PARAMETER (PI=3.141592654,PI4=4.*PI)
        INTEGER IA
        CHARACTER C*(MPC)
        INTEGER IC,IS,I,M,N,J,L
        REAL D,DR,PMIN,PMAX,R,PI4R,S,DP
        REAL RADIUS
        REAL PC(MPC)
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
C====   Print profile of Patterson peak for group form-factor option.
C
        IC=ICF(IA)
        DO 10 IS=1,MST
10        TF(IS,IA)=FF(IS,IC)*EXP(-MIN(BAT(IA)*SS(IS),88.))
        D=RADIUS(IA,IA)
        WRITE (6,'(////2A,I4,2(A,F8.2)/A,F8.2/)')
     &  ' Patterson peak profile of group ',AS(IA),IAT(IA),
     &  ' for dmin =',DMIN,'    B-factor =',BAT(IA),
     &  ' Horizontal line is zero density;'//
     &  ' vertical line is effective radius =',D
C
        DR=4./3.*D/MPC
        PMIN=0.
        PMAX=0.
        DO 30 I=1,MPC
          R=(I-.5)*DR
          PI4R=PI4*R
          S=0.
          DO 20 IS=1,MST
20          S=S+TF(IS,IA)**2*SIN(PI4R*ST(IS))*ST(IS)
          PC(I)=STH*S/R
C         WRITE (6,'(A,2F10.2)') ' R, PC =',R,PC(I)
          PMIN=MIN(PMIN,PC(I))
30        PMAX=MAX(PMAX,PC(I))
        DP=(PMAX-PMIN)/40
        M=NINT(.75*MPC)
        N=NINT(PMIN/DP)
C
        DO 50 J=39+N,N,-1
          DO 40 I=1,MPC
            IF (NINT(PC(I)/DP).GT.J) THEN
              C(I:I)='*'
              L=I
            ELSE
              IF (I.EQ.M) THEN
                C(I:I)='|'
                L=MAX(L,M)
              ELSE
                C(I:I)=' '
              ENDIF
            ENDIF
40        CONTINUE
          WRITE (6,'(1X,A)') C(:L)
50        IF (J.EQ.0) WRITE (6,'(133A)') '+',('_',I=1,MPC)
        END
C
C
        REAL FUNCTION RADIUS(IA,JA)
C
        IMPLICIT NONE
        INTEGER MAT,MCF,MST
        PARAMETER (MAT=50,MCF=10,MST=8)
        REAL PI,PI4
        PARAMETER (PI=3.141592654,PI4=4.*PI)
        INTEGER IA,JA
        INTEGER IS,N
        REAL P0,PI4R,R,S,SA,SP,SPP,T,A,TS,P
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
C====   If RCUT < 0, just take |RCUT| as radius; useful if iteration below
C       converges to wrong zero!
C
        IF (RCUT.LT.0.) THEN
          RADIUS=-RCUT
          RETURN
        ENDIF
C
C====   Compute radius of peak; make initial estimate then do Newton-Raphson
C       iteration to get minimum P^2 (ie zero P or zero dP/dR).
C
CD        WRITE (6,'(A,I4,A/1X,8F10.5)') ' RADIUS:',IA,'  TF =',
CD     &  (TF(IS,IA),IS=1,MST)
CD        WRITE (6,'(A,I4,A/1X,8F10.5)') ' RADIUS:',JA,'  TF =',
CD     &  (TF(IS,JA),IS=1,MST)
C
        P0=0.
        DO 10 IS=1,MST
10        P0=P0+TF(IS,IA)*TF(IS,JA)*SS(IS)
        P0=1.E-5*(PI4*P0)**2
C
        R=R0(ICF(IA))+R0(ICF(JA))
CD        WRITE (6,'(A,F8.2)') ' R0 =',R
        DO 30 N=1,100
          PI4R=PI4*R
          S=0.
          SP=0.
          SPP=0.
          DO 20 IS=1,MST
            A=PI4R*ST(IS)
            T=TF(IS,IA)*TF(IS,JA)
            TS=T*SIN(A)*ST(IS)
            S=S+TS
            SP=SP+T*COS(A)*SS(IS)
20          SPP=SPP+TS*SS(IS)
          S=S/R
          SP=S-PI4*SP
          SPP=MAX(SP*SP+S*(2.*SP-PI4*PI4R*SPP),500.)
          R=MIN(MAX(R*ABS(1.+S*SP/SPP),R-1.,1.),R+.5)
CD          T=S*SP/P0
CD          WRITE (6,'(A,I3,3F10.2,F10.0,1PE12.3)')
CD     &    ' N, S, SP, SPP, R, T =',N,S,SP,SPP,R,T
30        IF (ABS(S*SP).LT.P0) GOTO 40
C
C====   If RCUT > 0, take it as factor*RMS radius to cut off.
C
40      IF (RCUT.GT.0.) THEN
          PI4R=PI4*R
          S=0.
          SP=0.
          DO 50 IS=1,MST
            A=PI4R*ST(IS)
            SA=SIN(A)
            T=(TF(IS,IA)*TF(IS,JA))/A
            P=SA-A*COS(A)
            S=S+T*P
50          SP=SP+T*(P*(1.-6./A**2)+2.*SA)
          R=R*MIN(RCUT*SQRT(SP/S),1.)
CD          WRITE (6,'(A,2F10.2,2F8.2)') ' S, SP, RCUT, R =',S,SP,RCUT,R
        ENDIF
        RADIUS=R
        END
C
C
        SUBROUTINE GENPTS(IPTS,ISORT)
C
C====   Generate list of points within peaks in the Patterson.
C
        IMPLICIT NONE
        INTEGER MAT,MCF,MST,MAXI
        PARAMETER (MAT=50,MCF=10,MST=8)
        PARAMETER (MAXI=(2**30-1)+2**30)
        REAL PI,PI4
        PARAMETER (PI=3.141592654,PI4=4.*PI)
C
        INTEGER IA,IC,IS,JA,IE,I,J,I3,I2,I1,IM,IT,II,JE,IN
        REAL D,DD,S,RR
        INTEGER ICEIL,IFLOOR
        REAL RADIUS
        INTEGER LB(3),MB(3),IXP(3),JXP(3),IUP(3),JUP(3)
        REAL XV(3),XP(3)
C
        INTEGER MPT,NPT,NCYC(3),IVF(4)
        COMMON /REFCOM/ MPT,NPT,NCYC,IVF
        INTEGER*2 IPTS(3,3,MPT)
        INTEGER ISORT(2,MPT)
C
        INTEGER IUVW(3),LUVW(3),MUVW(3),NUVW(3),NXYZ(3)
        INTEGER LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
C
        REAL VOL,OBF,RCG(3),TM(3,3)
        COMMON /CELCOM/ VOL,OBF,RCG,TM
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        LUVW(1) = LU
        LUVW(2) = LV
        LUVW(3) = LW
        MUVW(1) = MU
        MUVW(2) = MV
        MUVW(3) = MW
C
C====   Compute the atomic temperature factors.
C
        DO 10 IA=1,NAT
          IC=ICF(IA)
          DO 10 IS=1,MST
10          TF(IS,IA)=FF(IS,IC)*EXP(-MIN(BAT(IA)*SS(IS),88.))
C
C====   Compute the interatomic vectors.
C
        NPT=0
        DO 150 IA=1,NAT
          DO 150 JA=IA,NAT
            D=RADIUS(IA,JA)
            DD=D*D
C
C====   Generate equivalent positions for the second atom.
C
            DO 150 IE=1,NGEP
C
C====   Compute coordinates of vector.
C
              DO 70 I=1,3
                XV(I)=XAT(I,IA)-RTE(I,4,IE)
                DO 60 J=1,3
60                XV(I)=XV(I)-RTE(I,J,IE)*XAT(J,JA)
C
C====   Compute size and limits of box bounding this peak.
C
                S=D*RCG(I)
                LB(I)=ICEIL(XV(I)-S)
70              MB(I)=IFLOOR(XV(I)+S)
C              WRITE (6,'(A,3I4,4F8.2,6I4)')
C     &        ' VECTOR',IA,JA,IE,D,XV,(LB(I),MB(I),I=1,3)
C
C====   Generate points inside box.
C
CD            RN=DD+.1
              DO 140 I3=LB(3),MB(3)
                IXP(3)=I3
                XP(3)=IXP(3)-XV(3)
                DO 140 I2=LB(2),MB(2)
                  IXP(2)=I2
                  XP(2)=IXP(2)-XV(2)
                  DO 140 I1=LB(1),MB(1)
                    IXP(1)=I1
                    XP(1)=IXP(1)-XV(1)
C
C====   Compute distance of point from peak centre.
C
                    RR=0.
                    DO 80 I=1,3
                      DO 80 J=1,3
80                      RR=RR+TM(I,J)*XP(I)*XP(J)
C
                    IF (RR.LE.DD) THEN
C
C====   Compute asymmetric unit map coordinates of point.
C
                      IM=MAXI
                      DO 120 IT=1,NTC
                        DO 120 II=-1,1,2
                          DO 120 JE=1,NGEP
                            DO 90 I=1,3
                              JXP(I)=0
                              DO 90 J=1,3
90                              JXP(I)=JXP(I)+NINT(RTE(I,J,JE))*IXP(J)
C
C====   Apply inversion centre and centring translations to get point
C       in asymmetric unit within map boundaries.
C
                            DO 100 I=1,3
                              J=IUVW(I)
                              JUP(I)=II*JXP(J)+ITC(J,IT,LTY)
                              IF (JUP(I).LT.LUVW(I)) THEN
                                JUP(I)=JUP(I)+NXYZ(J)*
     &                          ((LUVW(I)-JUP(I)-1)/NXYZ(J)+1)
                                IF (JUP(I).GT.MUVW(I)) GOTO 120
                              ELSEIF (JUP(I).GT.MUVW(I)) THEN
                                JUP(I)=JUP(I)-NXYZ(J)*
     &                          ((JUP(I)-MUVW(I)-1)/NXYZ(J)+1)
                                IF (JUP(I).LT.LUVW(I)) GOTO 120
                              ENDIF
100                         CONTINUE
C
C====   Pack map coordinates into one word and find minimum to make it unique.
C
                            IN=JUP(1)+NUVW(1)*(JUP(2)+NUVW(2)*JUP(3))
                            IF (IN.LT.IM) THEN
                              IM=IN
                              DO 110 I=1,3
110                             IUP(I)=JUP(I)
                            ENDIF
120                   CONTINUE
                      IF (IM.EQ.MAXI) THEN
                        WRITE (6,'(A,3I4)')
     +                       '*** ERROR: Point not found,'
     &                  //' check map symmetry',IXP
                        CALL CCPERR(1,'*** ERROR IN GENPTS.')
                      ENDIF
C
C====   Point is inside peak, store pointers.
C
                      NPT=NPT+1
                      IF (NPT.LT.MPT) THEN
                        IPTS(1,1,NPT)=IA
                        IPTS(2,1,NPT)=JA
                        IPTS(3,1,NPT)=IE
                        DO 130 I=1,3
                          IPTS(I,2,NPT)=IXP(I)
130                       IPTS(I,3,NPT)=IUP(I)
C
C====   Store sort key and index pointer.
C
                        ISORT(1,NPT)=IM
                        ISORT(2,NPT)=NPT
                      ENDIF
                    ENDIF
140           CONTINUE
150     CONTINUE
C
C====   Sort the points in the order they will come in the map.
C
        IF (NPT.EQ.0) CALL CCPERR(1,'*** ERROR IN GENPTS: no points.')
C
        IF (NPT.LT.MPT) THEN
          CALL BTSORT(NPT,ISORT)
C
C====   Put marker at end of sorted array.
C
          ISORT(1,NPT+1)=ISORT(1,1)-1
        ELSE
          WRITE (6,'(/1X,A,I9,8X,A,I9/1X,A)')
     &    'Dynamic array dimension =',MPT,'Number of points =',NPT,
     &    'Dimension will be increased & last step repeated.'
        ENDIF
        END
C
C
        SUBROUTINE BTSORT(N,IA)
C
C====   Binary tree sort.  IA(1,I) is sort key; IA(2,I) is usually table index.
C
        IMPLICIT NONE
        INTEGER N
        INTEGER IA(2,N)
        INTEGER L,M,I1,I2,I,J
C
        L=N/2+1
        M=N
10      IF (L.GT.1) THEN
          L=L-1
          I1=IA(1,L)
          I2=IA(2,L)
        ELSE
          I1=IA(1,M)
          IA(1,M)=IA(1,1)
          I2=IA(2,M)
          IA(2,M)=IA(2,1)
          M=M-1
          IF (M.LE.1) THEN
            IA(1,1)=I1
            IA(2,1)=I2
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF (J.LE.M) THEN
          IF (J.LT.M) THEN
            IF (IA(1,J).LT.IA(1,J+1)) J=J+1
          ENDIF
          IF (I1.LT.IA(1,J)) THEN
            IA(1,I)=IA(1,J)
            IA(2,I)=IA(2,J)
            I=J
            J=J+J
          ELSE
            J=M+1
          ENDIF
          GOTO 20
        ENDIF
        IA(1,I)=I1
        IA(2,I)=I2
        GOTO 10
        END
C
C
        SUBROUTINE REFCYC(IPTS,ISORT,SECT)
C
C====   Do a refinement cycle
C
        IMPLICIT NONE
        INTEGER MAT,MCF,MST,MVAR
        PARAMETER (MAT=50,MCF=10,MST=8,MVAR=5*MAT)
        REAL PI,PI4
        PARAMETER (PI=3.141592654,PI4=4.*PI)
C
        CHARACTER*3 CXYZ
        LOGICAL WF
        INTEGER IA,IC,NPAR,NVAR,NOB1,NOB2,NPU,NOB,IM,IV,JV,IR,IW,JT,KT,
     &  N,IT,JA,IE,I,J,IS,MP,IP,IER,NOB3
        REAL P000,SWOS1,SWES1,SWOS2,SWES2,SWES,PS,RR,R,PI4R,S,T,W,F,P,X,
     &  DPDR,DPDB,PC,PO,WO,WP,E,WE,WPS,ES,POC,POS,PCS,PS1,PS2,WS,XS,A,B,
     &  C,RMO1,RMC1,RMO2,RMC2,RMO3,RMC3,S1,S2,S3,C1,C2,C3,RMAX,SMAX,
     &  DRDX,SWESP,SFAC
        DOUBLEPRECISION SO1,SC1,SOC1,SOS1,SCS1,SO2,SC2,SOC2,SOS2,SCS2,
     &  SO3,SC3,SOC3,SOS3,SCS3
        INTEGER MULTIP
        REAL WEIGHT
        INTEGER IXP(3),IUP(3),NOBA(MAT)
        REAL XP(3),DRDV(3),DX(3),ER(MAT),SWESA(MAT),SWOSA(MAT),
     &  PA(MAT),DP000(MAT),AM(MVAR*(MVAR+1)/2),BV(MVAR),DV(MVAR),
     &  ESD(MVAR)
        DOUBLE PRECISION SOA(MAT),SCA(MAT),SOCA(MAT),SOSA(MAT),SCSA(MAT)
C
        INTEGER MPT,NPT,NCYC(3),IVF(4)
        COMMON /REFCOM/ MPT,NPT,NCYC,IVF
        INTEGER*2 IPTS(3,3,MPT)
        INTEGER ISORT(2,MPT)
C
        INTEGER IUVW(3),NUVW(3),NXYZ(3),LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
        REAL SECT(LU:MU,LV:MV)
C
        REAL VOL,OBF,RCG(3),TM(3,3)
        COMMON /CELCOM/ VOL,OBF,RCG,TM
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        SAVE SWESP,SFAC
C
        DATA ER/MAT*0./
        DATA CXYZ/'XYZ'/
        DATA SWESP,SFAC/1.7E38,1./
C
C====   Compute F000, F000^2/V and derivatives.
C
        P000=0.
        DO 10 IA=1,NAT
          IC=ICF(IA)
10        P000=P000+OCC(IA)*F0(IC)
C
        IF (IVF(3).EQ.1) THEN
          DO 20 IA=1,NAT
            IC=ICF(IA)
20          DP000(IA)=2.*P000*F0(IC)/VOL
        ENDIF
        P000=P000*P000/VOL
        WRITE (6,'(/A,F10.2)') ' F000 term =',P000
CD      WRITE (6,'(//1X,A/)')
CD    & 'Atoms                 g.e.p.     ix   iy   iz'//
CD    & '           Pobs       Pcalc'
C
C====   Number of variables.
C
        NPAR=3*IVF(1)+IVF(2)+IVF(3)
        NVAR=NPAR*NAT+IVF(4)
C
C====   Zero counters and sums for statistics.
C
        NOB1=0
C       SWO1=0.
C       SWE1=0.
        SWOS1=0.
        SWES1=0.
        SO1=0.
        SC1=0.
        SOC1=0.
        SOS1=0.
        SCS1=0.
C
        NOB2=0
C       SWO2=0.
C       SWE2=0.
        SWOS2=0.
        SWES2=0.
        SO2=0.
        SC2=0.
        SOC2=0.
        SOS2=0.
        SCS2=0.
C
        NPU=0
        NOB=0
        SWES=0.
C
        DO 25 IA=1,MAT
          NOBA(IA)=0
          SWESA(IA)=0.
          SWOSA(IA)=0.
          SOA(IA)=0.
          SCA(IA)=0.
          SOCA(IA)=0.
          SOSA(IA)=0.
25        SCSA(IA)=0.
C
C====   Zero normal matrix and vector.
C
        IM=0
        DO 30 IV=1,NVAR
          BV(IV)=0.
          DO 30 JV=1,IV
            IM=IM+1
30          AM(IM)=0.
C
C====   Read back points array.
C
        IR=ISORT(2,1)
        IW=IPTS(3,3,IR)-1
        JT=1
        DO 220 KT=1,NPT
C
C====   Check for overlaps.
C
          IF (ISORT(1,JT).EQ.ISORT(1,KT+1)) GOTO 220
C
C====   Initialise calculated Patterson density and derivative vector.
C
          PS=0.
          DO 35 IA=1,NAT
35          PA(IA)=0.
          DO 40 IV=1,NVAR
40          DV(IV)=0.
C
C====   Loop over overlapping points and get atom pointers.
C
          WF=.FALSE.
          N=0
          DO 140 IT=JT,KT
            IR=ISORT(2,IT)
            IA=IPTS(1,1,IR)
            JA=IPTS(2,1,IR)
            IE=IPTS(3,1,IR)
            IF (IA.EQ.JA .AND. IE.EQ.IDEN) WF=.TRUE.
C
C====   Get coordinates of point and compute vector.
C
            DO 50 I=1,3
              IXP(I)=IPTS(I,2,IR)
              XP(I)=XAT(I,IA)-RTE(I,4,IE)-IXP(I)
              DO 50 J=1,3
50              XP(I)=XP(I)-RTE(I,J,IE)*XAT(J,JA)
C
C====   Compute distance of point from peak centre.
C
            RR=0.
            DO 70 I=1,3
              DO 70 J=1,3
70              RR=RR+TM(I,J)*XP(I)*XP(J)
            R=SQRT(RR)
            PI4R=PI4*R
C
C====   Compute Patterson density at point.
C
            S=0.
            DO 80 IS=1,MST
              T=TF(IS,IA)*TF(IS,JA)
              IF (R.EQ.0.) THEN
                S=S+T*SS(IS)
              ELSE
                S=S+T*SIN(PI4R*ST(IS))*ST(IS)
              ENDIF
80          CONTINUE
C
            IF (S.LE.0.) GOTO 140
            N=N+1
C
            IF (R.EQ.0.) THEN
              S=PI4*STH*S
            ELSE
              S=STH*S/R
            ENDIF
C
            IF (IA.EQ.JA) THEN
C
C====   Get general Harker vector multiplicity for this symmetry operator.
C
              W=.5
            ELSE
              W=1.
            ENDIF
C
            F=W*OCC(IA)*OCC(JA)
            P=F*S
C            WRITE (2,'(10I5,3F10.2)') 1,IT,N,IR,IA,JA,IE,IXP,W,R,P
            PS=PS+P
            PA(IA)=PA(IA)+P*OCC(IA)
            PA(JA)=PA(JA)+P*OCC(JA)
C
C====   Test if excluding origin.
C
            IF (IOR .AND. WF) GOTO 140
C
C====   Derivatives with respect to atomic coordinates.
C
            IV=NPAR*(IA-1)
            JV=NPAR*(JA-1)
            IF (F.EQ.0.) THEN
              IF (IVF(1).EQ.1) THEN
                IV=IV+3
                JV=JV+3
              ENDIF
              IF (IVF(2).EQ.1) THEN
                IV=IV+1
                JV=JV+1
              ENDIF
            ELSE
              IF (IVF(1).EQ.1) THEN
                IF (RR.EQ.0.) THEN
                  IV=IV+3
                  JV=JV+3
                ELSE
                  DPDR=0.
                  DO 90 IS=1,MST
90                  DPDR=DPDR+TF(IS,IA)*TF(IS,JA)*COS(PI4R*ST(IS))*
     +                   SS(IS)
                  DPDR=(PI4*F*STH*DPDR-P)/RR
C
                  DO 100 I=1,3
                    DRDV(I)=0.
                    DO 100 J=1,3
100                   DRDV(I)=DRDV(I)+TM(I,J)*XP(J)
C
                  DO 120 I=1,3
                    IV=IV+1
                    DV(IV)=DV(IV)+DPDR*DRDV(I)
                    DRDX=0.
                    DO 110 J=1,3
110                   DRDX=DRDX+RTE(J,I,IE)*DRDV(J)
                    JV=JV+1
120                 DV(JV)=DV(JV)-DPDR*DRDX
                ENDIF
              ENDIF
C
C====   Derivatives with respect to thermal parameters.
C
              IF (IVF(2).EQ.1 .OR. IVF(4).EQ.1) THEN
                DPDB=0.
                DO 130 IS=1,MST
                  T=TF(IS,IA)*TF(IS,JA)
                  IF (R.EQ.0.) THEN
                    DPDB=DPDB+T*SS(IS)**2
                  ELSE
                    DPDB=DPDB+T*SIN(PI4R*ST(IS))*ST(IS)*SS(IS)
                  ENDIF
130             CONTINUE
C
                IF (R.EQ.0.) THEN
                  DPDB=-PI4*F*STH*DPDB
                ELSE
                  DPDB=-F*STH*DPDB/R
                ENDIF
C
                IF (IVF(2).EQ.1) THEN
                  IV=IV+1
                  DV(IV)=DV(IV)+DPDB
                  JV=JV+1
                  DV(JV)=DV(JV)+DPDB
                ELSE
                  DV(NVAR)=DV(NVAR)+DPDB
                ENDIF
              ENDIF
            ENDIF
C
C====   Derivatives with respect to occupancies.
C
            IF (IVF(3).EQ.1) THEN
              IV=IV+1
              DV(IV)=DV(IV)+W*OCC(JA)*S
              JV=JV+1
              DV(JV)=DV(JV)+W*OCC(IA)*S
            ENDIF
140       CONTINUE
C
C====   Compute multiplicity.
C
          IF (N.EQ.0) GOTO 210
          MP=MULTIP(NXYZ,NGEP,RTE,IXP)
          PC=MP*PS-P000
C
C====   Correct derivatives with multiplicity
C       and subtract origin contribution from occupancy derivatives.
C
          IF (NPAR.GT.0) THEN
            IV=0
            DO 160 IA=1,NAT
              DO 150 IP=1,NPAR
                IV=IV+1
150             DV(IV)=MP*DV(IV)
160           IF (IVF(3).EQ.1) DV(IV)=DV(IV)-DP000(IA)
          ENDIF
          IF (IVF(4).EQ.1) DV(NVAR)=MP*DV(NVAR)
C
C====   Get observed Patterson density.
C
          DO 170 I=1,3
170         IUP(I)=IPTS(I,3,IR)
          IF (IUP(3).NE.IW) THEN
            IW=IUP(3)
            CALL MPOSN(1,IW)
            CALL MGULPR(1,SECT,IER)
            IF (IER.NE.0) CALL CCPERR(1,
     &      '*** ERROR IN REFCYC: MGULPR error.')
C           WRITE (6,*) 'SECTION',IW,' READ'
          ENDIF
          PO=SECT(IUP(1),IUP(2))
C
          IF (ISORT(1,JT).EQ.0) THEN
            WRITE (6,'(/29X,A,44X,A)')
     &      'g.e.p.','Pobs       Pcalc'
            WRITE (6,'(1X,A,I6,39X,2F12.2)')
     &      'Origin peak               ',IDEN,PO,PC
          ENDIF
C
C===    Compute weight.
C
          WO=WEIGHT(NXYZ,NGEP,RTE,TM,DSTMIN,DSTMAX,IXP)
C          WRITE (2,'(10I5,F10.5,2F10.2)') 2,JT,KT,MP,IUP,IXP,WO,PO,PC
C
C====   Accumulate sums for least squares.
C
          WP=WO**2*ABS(PO+P000)
          E=WO*(PO-PC)
          WE=WO*ABS(E)
          WPS=(WO*PO)**2
          ES=E*E
          POC=PO*PC
          POS=PO*PO
          PCS=PC*PC
          IF (WF) THEN
            IF (PS.GT.0.) THEN
              NOB1=NOB1+1
C             SWO1=SWO1+WP
C             SWE1=SWE1+WE
              SWOS1=SWOS1+WPS
              SWES1=SWES1+ES
C
              SO1=SO1+PO
              SC1=SC1+PC
              SOC1=SOC1+POC
              SOS1=SOS1+POS
              SCS1=SCS1+PCS
            ENDIF
            IF (IOR) GOTO 210
          ELSE
            IF (PS.GT.0.) THEN
              NOB2=NOB2+1
C             SWO2=SWO2+WP
C             SWE2=SWE2+WE
              SWOS2=SWOS2+WPS
              SWES2=SWES2+ES
C
              SO2=SO2+PO
              SC2=SC2+PC
              SOC2=SOC2+POC
              SOS2=SOS2+POS
              SCS2=SCS2+PCS
C
              PS1=0.
              PS2=0.
              DO 194 IA=1,NAT
                PS1=PS1+PA(IA)
194             PS2=PS2+PA(IA)**2
C
              DO 195 IA=1,NAT
                W=PA(IA)/PS1
                IF (W.GT.0.) THEN
                  WS=PA(IA)**2/PS2
                  NOBA(IA)=NOBA(IA)+1
                  SWOSA(IA)=SWOSA(IA)+WS*WPS
                  SWESA(IA)=SWESA(IA)+WS*ES
C
                  SOA(IA)=SOA(IA)+W*PO
                  SCA(IA)=SCA(IA)+W*PC
                  SOCA(IA)=SOCA(IA)+WS*POC
                  SOSA(IA)=SOSA(IA)+WS*POS
                  SCSA(IA)=SCSA(IA)+WS*PCS
                ENDIF
195           CONTINUE
C
            ENDIF
          ENDIF
C
          NPU=NPU+N
          NOB=NOB+1
          SWES=SWES+ES
C
          IM=0
          DO 200 IV=1,NVAR
            DV(IV)=WO*DV(IV)
            BV(IV)=BV(IV)+DV(IV)*E
            DO 200 JV=1,IV
              IM=IM+1
200           AM(IM)=AM(IM)+DV(IV)*DV(JV)
C
210       JT=KT+1
220     CONTINUE
C
C====   Statistics.
C
        IF (NOB.EQ.0) CALL CCPERR(1,
     &  '*** ERROR IN REFCYC: no observations.')
        WRITE (6,'(//A,2I8)') ' Number of points sorted, used =',NPT,NPU
        NOB3=NOB1+NOB2
C
        IF (NOB3.GT.0) THEN
          SO3=(SO1+SO2)/NOB3
          SC3=(SC1+SC2)/NOB3
          SOC3=(SOC1+SOC2)/NOB3
          SOS3=(SOS1+SOS2)/NOB3
          SCS3=(SCS1+SCS2)/NOB3
C
          IF (NOB1.GT.0) THEN
            SO1=SO1/NOB1
            SC1=SC1/NOB1
            SOC1=SOC1/NOB1
            SOS1=SOS1/NOB1
            SCS1=SCS1/NOB1
          ENDIF
C
          IF (NOB2.GT.0) THEN
            SO2=SO2/NOB2
            SC2=SC2/NOB2
            SOC2=SOC2/NOB2
            SOS2=SOS2/NOB2
            SCS2=SCS2/NOB2
          ENDIF
C
          RMO1=SQRT(SOS1)
          RMC1=SQRT(SCS1)
          RMO2=SQRT(SOS2)
          RMC2=SQRT(SCS2)
          RMO3=SQRT(SOS3)
          RMC3=SQRT(SCS3)
C
C         R1=1.
C         IF (SWO1.GT.0.) R1=SWE1/SWO1
C         R2=SWE2/SWO2
C         R3=(SWE1+SWE2)/(SWO1+SWO2)
C
          S1=1.
          IF (SWOS1.GT.0.) S1=SQRT(SWES1/SWOS1)
          S2=SQRT(SWES2/SWOS2)
          S3=SQRT((SWES1+SWES2)/(SWOS1+SWOS2))
C
          C1=(SOS1-SO1*SO1)*(SCS1-SC1*SC1)
          IF (C1.GT.0.) C1=(SOC1-SO1*SC1)/SQRT(C1)
          C2=(SOS2-SO2*SO2)*(SCS2-SC2*SC2)
          IF (C2.GT.0.) C2=(SOC2-SO2*SC2)/SQRT(C2)
          C3=(SOS3-SO3*SO3)*(SCS3-SC3*SC3)
          IF (C3.GT.0.) C3=(SOC3-SO3*SC3)/SQRT(C3)
C
          WRITE (6,'(///2(/1X,A)/)')
     &  'Statistic                origin only   Patterson    Patterson',
     &  '                                      exc. origin  inc. origin'
          WRITE (6,'(1X,A,3(5X,I8))')
     &    'Number of observations ',NOB1,NOB2,NOB3
          WRITE (6,'(1X,A,3(3X,F10.2))')
     &    'Mean observed density  ',SO1,SO2,SO3
          WRITE (6,'(1X,A,3(3X,F10.2))')
     &    'Mean calculated density',SC1,SC2,SC3
          WRITE (6,'(1X,A,3(3X,F10.2))')
     &    'Rms observed density   ',RMO1,RMO2,RMO3
          WRITE (6,'(1X,A,3(3X,F10.2))')
     &    'Rms calculated density ',RMC1,RMC2,RMC3
C         WRITE (6,'(1X,A,3(5X,F8.3))')
C     &   'R-factor               ',R1,R2,R3
          WRITE (6,'(1X,A,3(5X,F8.3))')
     &    'RMS Residual           ',S1,S2,S3
          WRITE (6,'(1X,A,3(5X,F8.3))')
     &    'Correlation coefficient',C1,C2,C3
        ENDIF
C
C====   Solve equations and output new parameters.
C
        IF (NVAR.EQ.0) RETURN
        NOB=NINT(OBF*NOB)
        CALL EIGSOL(.0004,SWES,NOB,.FALSE.,NVAR,BV,IM,AM,DV,ESD)
C
C====   Test for increasing sum weighted square error and apply golden
C       section shift factor.
C
        IF (SWES.GE.SWESP) THEN
          SFAC=.382*SFAC
        ELSE
          SFAC=MIN(1.618*SFAC,1.)
        ENDIF
        SWESP=SWES
        WRITE (6,'(/A,F8.3)') ' Shift factor =',SFAC
C
C====   Write out results.
C
        WRITE (6,'(///1X,A/)')
     &  'Atom    Parameter      Init       Old     Shift'//
     &  '    Change       New       Esd     Residual  Correlation'
        RR=0.
        RMAX=2.*DMIN**2
        SMAX=0.
        SO1=0.
        SO2=0.
        IV=0
        IM=0
        DO 290 IA=1,NAT
C
C         R1=1.
          S1=1.
          C1=0.
          IF (NOBA(IA).GT.0) THEN
            IF (SWOSA(IA).GT.0.) S1=SQRT(SWESA(IA)/SWOSA(IA))
            SOA(IA)=SOA(IA)/NOBA(IA)
            SCA(IA)=SCA(IA)/NOBA(IA)
            C1=MAX((SOSA(IA)/NOBA(IA)-SOA(IA)**2)*
     &      (SCSA(IA)/NOBA(IA)-SCA(IA)**2),0.D0)
            IF (C1.GT.0.)
     &      C1=MIN(MAX((SOCA(IA)/NOBA(IA)-SOA(IA)*SCA(IA))/SQRT(C1),
     &      -1.D0),1.D0)
          ENDIF
C
          IF (IVF(1).EQ.0) THEN
            WRITE (6,'(1X,A,I4,71X,2F10.3)') AS(IA),IAT(IA),S1,C1
          ELSE
            IF (OCC(IA).EQ.0.) THEN
              WRITE (6,'(/1X,A,I4,71X,2F10.3)') AS(IA),IAT(IA),S1,C1
              IM=IM+3*IV+6
              IV=IV+3
            ELSE
C
C====   Limit total change to Dmin*sqrt(2).
C
              DO 230 I=1,3
                DV(IV+I)=SFAC*DV(IV+I)
230             XP(I)=XAT(I,IA)-XATS(I,IA)
C
              A=0.
              B=0.
              C=0.
              DO 250 I=1,3
                DX(I)=0.
                DO 240 J=1,3
                  DX(I)=DX(I)+TM(I,J)*DV(IV+J)
240               C=C+TM(I,J)*XP(I)*XP(J)
                A=A+DV(IV+I)*DX(I)
250             B=B+XP(I)*DX(I)
C
              RR=A+2.*B+C
CD            WRITE (6,'(1X,6(1PE14.6))')
CD    &       XP,(DV(IV+I),I=1,3),A,B,C,RMAX,RR,OCC(IA)
              IF (A.EQ.0.) THEN
                S=0.
                R=SQRT(RR)
                RR=MIN(RR,RMAX)
              ELSEIF (RR.LE.RMAX) THEN
                S=1.
                R=SQRT(RR)
              ELSE
                S=(SQRT(MAX(B*B+A*(RMAX-C),0.))-B)/A
                R=SQRT(RMAX)
              ENDIF
C
              E=0.
              DO 270 I=1,3
                IM=IM+IV
                DO 260 J=1,I-1
                  IM=IM+1
CD                WRITE (6,*) IV,IM,I,J
260               E=E+2.*AM(IM)*DX(I)*DX(J)
                IM=IM+1
CD              WRITE (6,*) IV,IM,I
270             E=E+AM(IM)*DX(I)**2
              IF (A.GT.0.) THEN
                ER(IA)=S*SQRT(E/A)
                A=S*SQRT(A)
              ENDIF
              WRITE (6,'(/1X,A,I4,29X,2F10.2,10X,F10.2,2X,2F10.3)')
     &        AS(IA),IAT(IA),A,R,ER(IA),S1,C1
C
              DO 280 I=1,3
                IV=IV+1
                DV(IV)=S*DV(IV)
                IF (ESD(IV).GT.0.) SMAX=MAX(SMAX,ABS(DV(IV))/ESD(IV))
                XS=SXYZ(I)*XATS(I,IA)
                X=SXYZ(I)*XAT(I,IA)
                XAT(I,IA)=XAT(I,IA)+DV(IV)
                DV(IV)=SXYZ(I)*DV(IV)
                B=X+DV(IV)
                C=B-XS
                ESD(IV)=SXYZ(I)*ESD(IV)
                IF (FRACT) THEN
                  WRITE (6,'(15X,A,2X,6F10.4)')
     &            CXYZ(I:I),XS,X,DV(IV),C,B,ESD(IV)
                ELSE
                  WRITE (6,'(15X,A,2X,6F10.2)')
     &            CXYZ(I:I),XS,X,DV(IV),C,B,ESD(IV)
                ENDIF
280           CONTINUE
            ENDIF
          ENDIF
C
          IF (IVF(2).EQ.1 .OR. IVF(4).EQ.1) THEN
            IF (IVF(2).EQ.1) THEN
              IV=IV+1
              JV=IV
              IM=IM+IV
            ELSE
              JV=NVAR
            ENDIF
C
            IF (OCC(IA).NE.0. .OR. IVF(4).EQ.1) THEN
              A=SFAC*DAMP*DV(JV)
              IF (ESD(JV).GT.0.) SMAX=MAX(SMAX,ABS(A)/ESD(JV))
              B=MIN(MAX(BAT(IA)+A,-99.),999.)
              C=B-BATS(IA)
              WRITE (6,'(15X,A,2X,6F10.2)')
     &        'B',BATS(IA),BAT(IA),A,C,B,ESD(JV)
              BAT(IA)=B
            ENDIF
          ENDIF
C
          SO1=SO1+ABS(OCC(IA))
          IF (IVF(3).EQ.1) THEN
            IV=IV+1
            DV(IV)=SFAC*DV(IV)
            IF (IVF(2).EQ.1) DV(IV)=DAMP*DV(IV)
            IF (ESD(IV).GT.0.) SMAX=MAX(SMAX,ABS(DV(IV))/ESD(IV))
            IF (RR.GT.RMAX) THEN
              B=0.
            ELSE
              B=OCC(IA)+DV(IV)
            ENDIF
            C=B-OCCS(IA)
            WRITE (6,'(15X,A,6F10.2)')
     &      'Occ',OCCS(IA),OCC(IA),DV(IV),C,B,ESD(IV)
            OCC(IA)=B
            ESDOCC(IA)=ESD(IV)
            SO2=SO2+ABS(OCC(IA))
            IM=IM+IV
          ENDIF
290     CONTINUE
        IF (SO1.EQ.0.0 .AND. SO2.EQ.0.0) CALL CCPERR(1,
     &  '*** ERROR IN REFCYC: all occupancies zero.')
        WRITE (6,'(/A,F10.2)') ' Maximum shift/esd =',SMAX
        END
C
C
        REAL FUNCTION WEIGHT(NXYZ,NGEP,RTE,TM,DSTMIN,DSTMAX,IX)
C
C===    Compute reciprocal of sqrt(variance) of density.
C
        IMPLICIT NONE
        INTEGER NGEP
        REAL DSTMIN,DSTMAX
        INTEGER NXYZ(3),IX(3)
        REAL RTE(3,4,24),TM(3,3)
        INTEGER I,IE,J,K
        REAL G,RR,S,XMIN,XMAX
        INTEGER NH(3),JX(3),IV(3)
C
        DO 10 I=1,3
10        NH(I)=NXYZ(I)/2
C
        S=0.
        DO 50 IE=1,NGEP
          DO 20 I=1,3
            JX(I)=0
            DO 20 J=1,3
20            JX(I)=JX(I)+NINT(RTE(I,J,IE))*IX(J)
          DO 50 K=-1,1,2
            DO 30 I=1,3
              IV(I)=MOD(IX(I)+K*JX(I),NXYZ(I))
              IF (IV(I).LT.-NH(I)) THEN
                IV(I)=IV(I)+NXYZ(I)
              ELSEIF (IV(I).GT.NH(I)) THEN
                IV(I)=IV(I)-NXYZ(I)
              ENDIF
30          CONTINUE
            RR=0.
            DO 40 I=1,3
              DO 40 J=1,3
40              RR=RR+TM(I,J)*IV(I)*IV(J)
            RR=SQRT(RR)
            XMIN=DSTMIN*RR
            XMAX=DSTMAX*RR
            G=XMAX**3-XMIN**3
            IF (G.EQ.0.) THEN
              G=1.
            ELSE
              G=3.*((SIN(XMAX)-XMAX*COS(XMAX))-
     &              (SIN(XMIN)-XMIN*COS(XMIN)))/G
            ENDIF
50          S=S+G
        WEIGHT=1./SQRT(S)
        END
C
C
        SUBROUTINE EIGSOL
     &  (ELIM, SWES, NOB, DBUG, NVAR, BV, NA, AM, DV, ESD)
C
C====   Solution of linear equations by eigenvalue filtering.
C
C
C       Ian J. Tickle           Birkbeck
C
C
C       Imported arguments:
C       ==================
C
C       ELIM    Eigenvalue limit for filtering.
C
C       SWES    Sum of weighted squared differences.
C
C       NOB     Number of observations.
C
C       DBUG    Flag to turn on some debugging output (LOGICAL).
C
C       NVAR    Number of variables.
C
C       BV      Normal equation vector.
C
C       NA      Dimension of matrix AM.
C
C       Modified arguments.
C       ==================
C
C       AM      Array used to store normal matrix and to return the
C               variance/covariance matrix.  Only one triangle is stored,
C               so the dimension of A must be at least NVAR * (NVAR + 1) / 2 .
C               The storage order is ((AM(I,J), J = 1, I), I = 1, NVAR) .
C
C       Exported arguments.
C       ==================
C
C       DV      Parameter shift vector.
C
C       ESD     Standard deviation vector.
C
C
        IMPLICIT NONE
        INTEGER         MVAR
        PARAMETER      (MVAR=250)
        INTEGER         NOB, NVAR, NA
        LOGICAL         DBUG
        REAL            ELIM, SWES
        REAL            BV(NVAR), AM(NA), DV(NVAR), ESD(NVAR)
C
        REAL            RM(MVAR*MVAR), SV(MVAR), EV(MVAR)
        INTEGER         IB(10), JB(10)
        REAL            CB(10)
        INTEGER         IA, IV, JA, JV, J, I, NNZ, K, KV, KB
        REAL            DM, SUMP, AIJ, RMSWD, CC
C
C
C       Solve matrix equation.  First get the diagonal elements of the normal
C       matrix and use them to pre-condition the vector and matrix to minimise
C       the effects of round-off errors due to large differences in the
C       magnitudes of the eigenvalues.
C
        DM = 0.
        IA = 0
        DO 10 IV = 1, NVAR
          IA = IA + IV
10        DM = MAX(DM, AM(IA))
        DM = DM * 1.E-12
C
CD      WRITE(6,'(/A)') ' VECTOR AND MATRIX'
        IA = 0
        DO 20 IV = 1, NVAR
          JA = IA + IV
CD        WRITE(6,'(1X,I2,(1X,10(1PE11.3)))') IV,BV(IV),(AM(I),I=IA+1,JA)
          IF (AM(JA) .LE. DM) THEN
            SV(IV) = 0.
          ELSE
            SV(IV) = 1. / SQRT(AM(JA))
          ENDIF
          BV(IV) = SV(IV) * BV(IV)
          DO 20 JV = 1, IV
            IA = IA + 1
20          AM(IA) = SV(IV) * SV(JV) * AM(IA)
C
CD      WRITE(6,'(/A)') ' SCALED VECTOR AND MATRIX'
CD      IA = 0
CD      DO 30 IV = 1, NVAR
CD        JA = IA + IV
CD        WRITE(6,'(1X,I2,(1X,10(1PE11.3)))') IV,BV(IV),(AM(I),I=IA+1,JA)
CD30      IA = JA
C
C       Now calc eigenvalues & eigenvectors.
C
        CALL EIGEN(AM, RM, NVAR, 0)
C
C       Find largest eigenvalue.
C
        DM = 0.
        IA = 0
        DO 40 IV = 1, NVAR
          IA =IA + IV
          EV(IV) = AM(IA)
40        DM = MAX(DM, ABS(EV(IV)))
        DM=ELIM*DM
C
        IF (DBUG) THEN
          WRITE(6,'(/A)') ' EIGENVALUES'
          WRITE(6,'((1X,10(1PE11.3)))') (EV(IV), IV = 1, NVAR)
          WRITE(6,'(/A)') ' EIGENVECTORS'
          DO 50 J = 0, NVAR-1
50          WRITE(6,'(1X,10F11.4)') (RM(I+J), I=1,NVAR**2,NVAR)
        ENDIF
C
C       Invert the eigenvalues with filtering.
C
        NNZ = 0
        DO 60 IV = 1, NVAR
          IF (ABS(EV(IV)) .LE. DM) THEN
            EV(IV) = 0.
          ELSE
            EV(IV) = 1. / EV(IV)
            NNZ = NNZ + 1
        ENDIF
60      CONTINUE
C
        IF (DBUG) THEN
          WRITE(6,'(/A)') ' INVERSE EIGENVALUES'
          WRITE(6,'((1X,10(1PE11.3)))') (EV(IV), IV = 1, NVAR)
        ENDIF
C
C       Form the elements of the inverse normal matrix.
C
        SUMP = 0.
        IA = 0
        DO 90 IV = 1, NVAR
          DV(IV) = 0.
          DO 80 JV = 1, NVAR
            AIJ = 0.
            K = 0
            DO 70 KV = 1, NVAR
              AIJ = AIJ + RM(IV + K) * EV(KV) * RM(JV + K)
70            K= K + NVAR
C
C       Save inverse normal matrix for correlations and e.s.d.'s
C
            IF (IV .GE. JV) THEN
              IA = IA + 1
              AM(IA) = AIJ
            ENDIF
C
C       Shift vector = inverse matrix . difference vector
C
80          DV(IV) = DV(IV) + AIJ * BV(JV)
C
          SUMP = SUMP + DV(IV) * BV(IV)
90        DV(IV) = SV(IV) * DV(IV)
C
C       Calculate standard error of fit.
C
CD      WRITE(6, *) 'SWES, SUMP =', SWES, SUMP
        RMSWD = SQRT(MAX(SWES - SUMP, 0.) / MAX(NOB - NNZ,1))
C
CD      WRITE(6,'(/A)') ' INVERSE MATRIX'
CD      IA = 0
CD      DO 100 IV = 1, NVAR
CD        JA = IA + IV
CD        WRITE(6,'(1X,I2,(1X,10F10.4))') IV,(AM(I),I=IA+1,JA)
CD100     IA = JA
C
C       Calculate e.s.d.'s and the correlation matrix.
C
        IA = 0
        KB = -1
        DO 110 IV = 1, NVAR
          ESD(IV) = MAX( SQRT(AM(IA + IV)), 1.E-6)
          DO 110 JV = 1, IV
            IA = IA + 1
            CC = AM(IA) / (ESD(IV) * ESD(JV))
            IF (IV .NE. JV .AND. ABS(CC) .GT. 0.9) THEN
              IF (KB .LT. 0) THEN
                WRITE(6,'(/A)') ' CORRELATION MATRIX ELEMENTS > 0.9'
                KB = 0
              ENDIF
              KB = KB + 1
              IB(KB) = IV
              JB(KB) = JV
              CB(KB) = CC
              IF (KB .EQ. 6) THEN
                WRITE(6, '(1X,6(4X,2I4,F8.4))')
     &          (JB(KB), IB(KB), CB(KB), KB = 1, 6)
                KB = 0
              ENDIF
            ENDIF
110     CONTINUE
        IF (KB .GT. 0) WRITE(6, '(1X,6(4X,2I4,F8.4))')
     &  (JB(I), IB(I), CB(I), I = 1, KB)
C
C       Scale the variance / covariance matrix.
C
        IA = 0
        DO 120 IV = 1, NVAR
          SV(IV) = RMSWD * SV(IV)
          ESD(IV) = SV(IV) * ESD(IV)
          DO 120 JV = 1, IV
            IA = IA + 1
120         AM(IA) = SV(IV) * SV(JV) * AM(IA)
        END
C
C
        SUBROUTINE EIGEN(A, R, N, MV)
C
C---- SUBROUTINE TO COMPUTE EIGENVALUES & EIGENVECTORS OF A REAL SYMMETRIC
C---- MATRIX, FROM IBM SSP MANUAL (SEE P165).
C---- DESCRIPTION OF PARAMETERS -
C---- A - ORIGINAL MATRIX STORED COLUMNWISE AS UPPER TRIANGLE ONLY,
C---- I.E. "STORAGE MODE" = 1.  EIGENVALUES ARE WRITTEN INTO DIAGONAL
C---- ELEMENTS OF A  I.E.  A(1)  A(3)  A(6)  FOR A 3*3 MATRIX.
C---- R - RESULTANT MATRIX OF EIGENVECTORS STORED COLUMNWISE IN SAME
C---- ORDER AS EIGENVALUES.
C---- N - ORDER OF MATRICES A & R.
C---- MV = 0 TO COMPUTE EIGENVALUES & EIGENVECTORS.
C
        IMPLICIT NONE
        REAL            A(*), R(*)
        INTEGER         N, MV
C
        INTEGER         IQ, J, I, IJ, IA, IND, L, M, MQ, LQ, LM, LL, MM
        INTEGER         ILQ, IMQ, IM, IL, ILR, IMR
        REAL            ANORM, ANRMX, RANGE, THR, X, Y, SINX, SINX2,
     &                  COSX, COSX2, SINCS
C
C       DATA RANGE/1.D-12/
        DATA RANGE/1.E-6/
C
        IF (MV.EQ.0) THEN
          IQ=-N
          DO 20 J=1,N
            IQ=IQ+N
            DO 20 I=1,N
              IJ=IQ+I
              IF (I.EQ.J) THEN
                R(IJ)=1.
              ELSE
                R(IJ)=0.
              ENDIF
20        CONTINUE
        ENDIF
C
C---- INITIAL AND FINAL NORMS (ANORM & ANRMX)
        IA=0
        ANORM=0.
        DO 35 I=1,N
          DO 35 J=1,I
            IA=IA+1
35          IF (J.NE.I) ANORM=ANORM+A(IA)**2
C
        IF (ANORM.LE.0.) GOTO 165
        ANORM=SQRT(2.*ANORM)
        ANRMX=ANORM*RANGE/N
C
C---- INITIALIZE INDICATORS AND COMPUTE THRESHOLD
        IND=0
        THR=ANORM
45      THR=THR/N
50      L=1
55      LQ=L*(L-1)/2
        LL=L+LQ
        M=L+1
        ILQ=N*(L-1)
C
C---- COMPUTE SIN & COS
60      MQ=M*(M-1)/2
        LM=L+MQ
        IF (ABS(A(LM))-THR.LT.0.) GOTO 130
        IND=1
        MM=M+MQ
        X=.5*(A(LL)-A(MM))
        Y=-A(LM)/SQRT(A(LM)**2+X*X)
        IF (X.LT.0.) Y=-Y
        SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y))))
        SINX2=SINX**2
        COSX=SQRT(1.-SINX2)
        COSX2=COSX**2
        SINCS=SINX*COSX
C
C---- ROTATE L & M COLUMNS
        IMQ=N*(M-1)
        DO 125 I=1,N
          IQ=I*(I-1)/2
          IF (I.NE.L .AND. I.NE.M) THEN
            IF(I.LT.M) THEN
              IM=I+MQ
            ELSE
              IM=M+IQ
            ENDIF
            IF (I.LT.L) THEN
              IL=I+LQ
            ELSE
              IL=L+IQ
            ENDIF
            X=A(IL)*COSX-A(IM)*SINX
            A(IM)=A(IL)*SINX+A(IM)*COSX
            A(IL)=X
          ENDIF
          IF (MV.EQ.0) THEN
            ILR=ILQ+I
            IMR=IMQ+I
            X=R(ILR)*COSX-R(IMR)*SINX
            R(IMR)=R(ILR)*SINX+R(IMR)*COSX
            R(ILR)=X
          ENDIF
125     CONTINUE
C
        X=2.*A(LM)*SINCS
        Y=A(LL)*COSX2+A(MM)*SINX2-X
        X=A(LL)*SINX2+A(MM)*COSX2+X
        A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
        A(LL)=Y
        A(MM)=X
C
C---- TESTS FOR COMPLETION
C---- TEST FOR M = LAST COLUMN
130     IF (M.NE.N) THEN
          M=M+1
          GOTO 60
        ENDIF
C
C---- TEST FOR L =PENULTIMATE COLUMN
        IF (L.NE.N-1) THEN
          L=L+1
          GOTO55
        ENDIF
        IF (IND.EQ.1) THEN
          IND=0
          GOTO50
        ENDIF
C
C---- COMPARE THRESHOLD WITH FINAL NORM
        IF (THR.GT.ANRMX) GOTO 45
165     RETURN
        END
C
C
        SUBROUTINE PKLIST(L,RHO)
C
C====   Write list of peaks in the Patterson.
C
        IMPLICIT NONE
        INTEGER MAT,MCF,MST,MPL,MAXI
        PARAMETER (MAT=50,MCF=10,MST=8,MPL=64000)
        PARAMETER (MAXI=(2**30-1)+2**30)
        REAL PI,PI4
        PARAMETER (PI=3.141592654,PI4=4.*PI)
C
        INTEGER IER,IA,IC,JA,IE,NP,NPL,I,J,I3,I2,I1,
     &  IS,IM,IT,II,JE,IN,IW,MP,IL,K,L
        REAL P000,DZ,D,DD,F,POM,PCM,SPO,SPC,S,RR,R,PI4R,T,PC,PO,FM,FN
        INTEGER ICEIL,IFLOOR,MULTIP
        REAL RADIUS
        INTEGER IPL(MPL),IXV(3),LB(3),MB(3),IXP(3),JXP(3),IUP(3),JUP(3)
        REAL XV(3),XP(3),XM(3),UP(3)
C
        INTEGER IUVW(3),LUVW(3),MUVW(3),NUVW(3),NXYZ(3)
        INTEGER LU,LV,LW,MU,MV,MW
        COMMON /MAPCOM/ IUVW,LU,LV,LW,MU,MV,MW,NUVW,NXYZ
        REAL RHO(LU:MU,LV:MV,LW:MW)
C
        REAL VOL,OBF,RCG(3),TM(3,3)
        COMMON /CELCOM/ VOL,OBF,RCG,TM
C
        LOGICAL IOR
        INTEGER LTY,NGEP,IDEN,NTC,ITC(3,4,7)
        REAL RTE(3,4,24),WH(24)
        COMMON /GEPCOM/ IOR,LTY,NGEP,IDEN,NTC,ITC,RTE,WH
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        INTEGER ICF(MAT)
        REAL DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST(MST),SS(MST),R0(MCF),
     &  F0(MCF),FF(MST,MCF),TF(MST,MAT)
        COMMON /SCFCOM/ ICF,DMIN,DMAX,DSTMIN,DSTMAX,STH,RCUT,ST,SS,R0,
     &  F0,FF,TF
C
        LUVW(1) = LU
        LUVW(2) = LV
        LUVW(3) = LW
        MUVW(1) = MU
        MUVW(2) = MV
        MUVW(3) = MW
C
C====   Read in Patterson map.
C
        CALL MPOSN(1,LW)
        DO 10 IW=LW,MW
          CALL MGULPR(1,RHO(LU,LV,IW),IER)
10        IF (IER.NE.0) CALL CCPERR(1,
     &    '**** ERROR IN PKLIST: MGULPR error.')
C
C====   Compute F000, F000^2/V.
C
        P000=0.
        DO 20 IA=1,NAT
          IC=ICF(IA)
20        P000=P000+OCC(IA)*F0(IC)
        P000=P000*P000/VOL
        WRITE (6,'(//1X,A/)')
     &  'Atoms                       g.e.p.'//
     &  '       x      y      z         r            Pobs       Pcalc'//
     &  '           <Pobs>     <Pcalc>'
C
C====   Compute the interatomic vectors.
C
        DO 240 IA=1,NAT
          IF (OCC(IA).EQ.0.) GOTO 240
          DO 230 JA=IA,NAT
            IF (OCC(JA).EQ.0.) GOTO 230
            DZ=RADIUS(IA,JA)
            D=.5*DZ
            DD=D*D
C
C====   Generate equivalent positions for the second atom.
C
            F=OCC(IA)*OCC(JA)
            DO 220 IE=1,NGEP
              IF (IA.EQ.JA .AND. IE.EQ.IDEN) GOTO 220
C
              NP=0
              NPL=0
              POM=0.
              PCM=0.
              SPO=0.
              SPC=0.
C
C====   Compute coordinates of vector.
C
              DO 80 I=1,3
                XV(I)=XAT(I,IA)-RTE(I,4,IE)
                DO 70 J=1,3
70                XV(I)=XV(I)-RTE(I,J,IE)*XAT(J,JA)
                IXV(I)=NINT(XV(I))
C
C====   Compute size and limits of box bounding this peak.
C
                S=D*RCG(I)
                LB(I)=ICEIL(XV(I)-S)
80              MB(I)=IFLOOR(XV(I)+S)
CD
CD      IF (IA.EQ.1 .AND. JA.EQ.1)
CD     &        WRITE (6,'(A,F8.2,I4,3F8.2,3I4,4X,3I4)')
CD     &        ' D,IE,XV,LB,MB',D,IE,XV,LB,MB
CD
C
C====   Generate points inside box.
C
              DO 160 I3=LB(3),MB(3)
                IXP(3)=I3
                XP(3)=IXP(3)-XV(3)
                DO 160 I2=LB(2),MB(2)
                  IXP(2)=I2
                  XP(2)=IXP(2)-XV(2)
                  DO 160 I1=LB(1),MB(1)
                    IXP(1)=I1
                    XP(1)=IXP(1)-XV(1)
C
C====   Compute distance of point from peak centre.
C
                    RR=0.
                    DO 90 I=1,3
                      DO 90 J=1,3
90                      RR=RR+TM(I,J)*XP(I)*XP(J)
C
                    IF (RR.LE.DD) THEN
C
C====   Compute Patterson density at point.
C
                      R=SQRT(RR)
                      PI4R=PI4*R
                      S=0.
                      DO 100 IS=1,MST
                        T=TF(IS,IA)*TF(IS,JA)
                        IF (R.EQ.0.) THEN
                          S=S+T*SS(IS)
                        ELSE
                          S=S+T*SIN(PI4R*ST(IS))*ST(IS)
                        ENDIF
100                   CONTINUE
                      IF (S.LE.0.) GOTO 160
C
                      IF (R.EQ.0.) THEN
                        S=PI4*STH*S
                      ELSE
                        S=STH*S/R
                      ENDIF
C
                      PC=F*S
                      NP=NP+1
                      SPC=SPC+PC
C
C====   Compute asymmetric unit map coordinates of point.
C
                      IM=MAXI
                      DO 140 IT=1,NTC
                        DO 140 II=-1,1,2
                          DO 140 JE=1,NGEP
                            DO 110 I=1,3
                              JXP(I)=0
                              DO 110 J=1,3
110                             JXP(I)=JXP(I)+NINT(RTE(I,J,JE))*IXP(J)
C
C====   Apply inversion centre and centring translations to get point
C       in asymmetric unit within map boundaries.
C
                            DO 120 I=1,3
                              J=IUVW(I)
                              JUP(I)=II*JXP(J)+ITC(J,IT,LTY)
                              IF (JUP(I).LT.LUVW(I)) THEN
                                JUP(I)=JUP(I)+NXYZ(J)*
     &                          ((LUVW(I)-JUP(I)-1)/NXYZ(J)+1)
                                IF (JUP(I).GT.MUVW(I)) GOTO 140
                              ELSEIF (JUP(I).GT.MUVW(I)) THEN
                                JUP(I)=JUP(I)-NXYZ(J)*
     &                          ((JUP(I)-MUVW(I)-1)/NXYZ(J)+1)
                                IF (JUP(I).LT.LUVW(I)) GOTO 140
                              ENDIF
120                         CONTINUE
C
C====   Pack map coordinates into one word and find minimum to make it unique.
C
                            IN=JUP(1)+NUVW(1)*(JUP(2)+NUVW(2)*JUP(3))
                            IF (IN.LT.IM) THEN
                              IM=IN
                              DO 130 I=1,3
130                             IUP(I)=JUP(I)
                            ENDIF
140                   CONTINUE
                      IF (IM.EQ.MAXI) THEN
                        WRITE (6,'(A,3I4)')
     +                       '*** ERROR: Point not found,'
     &                  //' check map symmetry',IXP
                        CALL CCPERR(1,'*** ERROR IN PKLIST.')
                      ENDIF
C
                      MP=MULTIP(NXYZ,NGEP,RTE,IXP)
                      PO=RHO(IUP(1),IUP(2),IUP(3))
                      IF (IXP(1).EQ.IXV(1) .AND. IXP(2).EQ.IXV(2) .AND.
     &                    IXP(3).EQ.IXV(3)) THEN
                        IF (IA.EQ.JA) PC=WH(IE)*PC
                        PCM=MP*PC-P000
                        POM=PO
                      ENDIF
C
                      DO 150 IL=1,NPL
150                     IF (IPL(IL).EQ.IM) GOTO 160
                      IF (NPL.EQ.MPL) CALL CCPERR(1,
     &                '*** ERROR IN PKLIST: array bound check (MPL).')
                      NPL=NPL+1
                      IPL(NPL)=IM
                      SPO=SPO+(PO+P000)/MP
                    ENDIF
160           CONTINUE
C
C====   Compute asymmetric unit map coordinates of point.
C
              FM=MAXI
              DO 200 IT=1,NTC
                DO 200 II=-1,1,2
                  DO 200 JE=1,NGEP
                    DO 170 I=1,3
                      XP(I)=0
                      DO 170 J=1,3
170                     XP(I)=XP(I)+RTE(I,J,JE)*XV(J)
C
C====   Apply inversion centre and centring translations to get point
C       in asymmetric unit within map boundaries.
C
                    DO 180 I=1,3
                      J=IUVW(I)
                      UP(I)=II*XP(J)+ITC(J,IT,LTY)
                      JUP(I)=NINT(UP(I))
                      IF (JUP(I).LT.LUVW(I)) THEN
                        K=NXYZ(J)*((LUVW(I)-JUP(I)-1)/NXYZ(J)+1)
                        UP(I)=UP(I)+K
                        IF (JUP(I)+K.GT.MUVW(I)) GOTO 200
                      ELSEIF (JUP(I).GT.MUVW(I)) THEN
                        K=NXYZ(J)*((JUP(I)-MUVW(I)-1)/NXYZ(J)+1)
                        UP(I)=UP(I)-K
                        IF (JUP(I)-K.LT.LUVW(I)) GOTO 200
                      ENDIF
180                 CONTINUE
C
C====   Pack map coordinates into one word and find minimum to make it unique.
C
                    FN=UP(1)+NUVW(1)*(UP(2)+NUVW(2)*UP(3))
                    IF (FN.LT.FM) THEN
                      FM=FN
                      DO 190 I=1,3
190                     XM(IUVW(I))=UP(I)
                    ENDIF
200           CONTINUE
              IF (FM.EQ.MAXI) THEN
                WRITE (6,'(A,3F8.2)') '*** ERROR: Point not found,'//
     &          ' check map symmetry',XP
                CALL CCPERR(1,'*** ERROR IN PKLIST.')
              ENDIF
C
              IF (.NOT.FRACT) THEN
                DO 210 I=1,3
210               XM(I)=SXYZ(I)*XM(I)
              ENDIF
C
              IF (NP.GT.0) THEN
                IF (IA.EQ.JA) SPO=SPO/WH(IE)
                SPO=SPO/NP
                SPC=SPC/NP
              ENDIF
              WRITE (6,'(1X,2(A,I4,7X),I2,4X,3F7.1,F10.2,2(4X,2F12.2))')
     &        AS(IA),IAT(IA),AS(JA),IAT(JA),IE,XM,DZ,POM,PCM,SPO,SPC
220         CONTINUE
230       CONTINUE
240     CONTINUE
        END
C
C
        INTEGER FUNCTION ICEIL(X)
C
C====   Integer Ceiling function of real X.
C       Smallest integer algebraically >= X.
C
        IMPLICIT NONE
        REAL X
C
        ICEIL=X
        IF (ICEIL.NE.X .AND. X.GT.0.) ICEIL=ICEIL+1
        END
C
C
        INTEGER FUNCTION IFLOOR(X)
C
C====   Integer Floor function of real X.
C       Largest integer algebraically <= X.
C
        IMPLICIT NONE
        REAL X
C
        IFLOOR=X
        IF (IFLOOR.NE.X .AND. X.LT.0.) IFLOOR=IFLOOR-1
        END
C
C
        SUBROUTINE WRITAT
C
C====   Write out atomic coordinates to disk.
C
        IMPLICIT NONE
        INTEGER MAT
        PARAMETER (MAT=50)
        INTEGER IFAIL,IA,I,NR
        REAL X(3)
C
        CHARACTER*4 AS(MAT)
        COMMON /TXTCOM/ AS
C
        LOGICAL FRACT
        INTEGER NAT,IAT(MAT)
        REAL DAMP,TRMS,BLIM(2),SXYZ(3),XATS(3,MAT),BATS(MAT),OCCS(MAT),
     &  XAT(3,MAT),BAT(MAT),OCC(MAT),ESDOCC(MAT)
        COMMON /ATMCOM/ FRACT,NAT,IAT,DAMP,TRMS,BLIM,SXYZ,XATS,BATS,
     &  OCCS,XAT,BAT,OCC,ESDOCC
C
        IFAIL=0
        CALL CCPDPN(2,'ATOUT','NEW','F',0,IFAIL)
C
        NR=0
        DO 20 IA=1,NAT
          IF ((TRMS.LE.0. .AND. ABS(OCC(IA)).LE.-TRMS*ESDOCC(IA)) .OR.
     &        (TRMS.GT.0. .AND. OCC(IA).LE.TRMS*ESDOCC(IA)) .OR.
     &        BAT(IA).LE.BLIM(1) .OR. BAT(IA).GT.BLIM(2)) THEN
            WRITE (6,'(2A,I4,A)')
     &      ' >>>> ATOM ',AS(IA),IAT(IA),' REJECTED. <<<<'
            OCC(IA)=0.
            NR=NR+1
          ELSE
            DO 10 I=1,3
10            X(I)=SXYZ(I)*XAT(I,IA)
            IF (FRACT) THEN
              WRITE (2,'(2(A,1X),I4,4F10.4,F10.2)')
     &        'ATOM',AS(IA),IAT(IA),OCC(IA),X,BAT(IA)
            ELSE
              WRITE (2,'(2(A,1X),I4,F10.4,4F10.2)')
     &        'ATOM',AS(IA),IAT(IA),OCC(IA),X,BAT(IA)
            ENDIF
          ENDIF
20      CONTINUE
        CLOSE(2)
        IF (NR.EQ.NAT) CALL CCPERR(1,'*** ERROR: All atoms rejected!')
        END