C
C $Id: bpair.F,v 1.5 1999/04/06 20:36:41 arjan Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE BPAIR(IERROR)
C
C     BONDED ATOM PAIRLIST GENERATION ROUTINE.  TWO ATOMS ARE STORED
C     IN THE BONDED ATOM PAIRLIST IF THEY EVER APPEAR TOGETHER IN THE
C     SAME SUBSYSTEM AND THEY ARE SEPARATED BY LESS THAN THE USER
C     PARAMETER CUTBOND.
C
C     FOR I=2,3,...,NATOMS, THE ATOMS THAT "BOND" WITH I AND WHICH HAVE
C     LOWER ATOM NUMBERS THAN I ARE:
C
C         IPAIR(K), FOR K=IP1(I),IP1(I)+1,...,IP1(I+1)-1
C
C     IF THE MAXIMUM ALLOWED NUMBER OF PAIRS IS EXCEEDED, IERROR IS
C     RETURNED WITH THE ATOM NUMBER I AT WHICH THE MAXIMUM WAS REACHED.
C     THIS ROUTINE SHOULD BE CALLED EACH TIME A NEW SET OF SUBSYSTEMS
C     IS DEFINED.
C
C     NOTE THAT THE COMMON BLOCK /PAIRIJ/ IS USED TEMPORARILY HERE TO
C     STORE THE OLD PAIRLIST AND POINTERS TO DIATOMIC BLOCKS IN THE
C     OLD DENSITY MATRIX.  THIS INFORMATION WILL BE USED IN INITP TO
C     CONSTRUCT A NEW DENSITY MATRIX FROM THE OLD ONE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
C
C
C     LOCAL:
C
      LOGICAL FIRST
      DATA FIRST /.TRUE./
      SAVE FIRST
C
      IF(FIRST)THEN
        FIRST = .FALSE.
C
C
      ELSE
c
c temporary -- may be used to test whether or not the old density
c              matrix is getting transferred properly into the new
c              one.  Changing bcut will cause a fairly drastic change
c              in the bonded atom pairlist, so that it will be obvious
c              whether the algorithm is working.
c
C-RDC c       write(0,*) '********* changing bcut ********'
C-RDC c       write(8,*) '********* changing bcut ********'
c       bcut = bcut - 1.0d0
c       bcut2 = bcut*bcut
c
c  end temorary
C
C       TEMPORARILY STORE OLD PAIRLIST AND DIATOMIC BLOCK POINTERS.
C
        DO 50 IATM=2,NATOMS
          J1 = IP1(IATM)
          J2 = IP1(IATM+1)-1
          IP1OLD(IATM) = J1
          IF(J2.LT.J1) GO TO 50
          DO 40 J=J1,J2
            IPOLD(J) = IPAIR(J)
            IJOLD(J) = IJMAT(J)
 40       CONTINUE
 50     CONTINUE
        IP1OLD(NATOMS+1) = IP1(NATOMS+1)
        IJMAX = IJMAT(IP1(NATOMS+1))-1
        IJOLD(IJMAX) = IJMAT(IJMAX)
      ENDIF
C
      IERROR = 0
      NPAIRS = 0
      DO 100 I=2,NATOMS
        IP1(I) = NPAIRS + 1
        IF(IATNUM(I).EQ.0) GO TO 100
C
C       POINTERS FOR SUBSYSTEMS BELONGING TO ATOM I:
C
        I1 = IS1ATM(I)
        I2 = IS1ATM(I+1)-1
        IF(I2.LT.I1) GO TO 100
C
C       LOCATE ATOMS < I THAT SHARE A SUBSYSTEM WITH ATOM I.
C
        DO 90 J=1,I-1
          IF(IATNUM(J).EQ.0) GO TO 90
C
C         SKIP I,J PAIR IF SEPARATED BY MORE THAN THE USER-SPECIFIED
C         CUTBOND PARAMETER.
C
          IF(CUTBND)THEN
            IF(PBC)THEN
              CALL PBCXYZ(I,J,XJ,YJ,ZJ)
            ELSE
              XJ = XYZ(1,J)
              YJ = XYZ(2,J)
              ZJ = XYZ(3,J)
            ENDIF
            RIJSQR = (XYZ(1,I)-XJ)**2
     .             + (XYZ(2,I)-YJ)**2
     .             + (XYZ(3,I)-ZJ)**2
            IF(RIJSQR.GT.BCUT2) GO TO 90
          ENDIF
C
C         POINTERS FOR SUBSYSTEMS BELONGING TO ATOM J:
C
          J1 = IS1ATM(J)
          J2 = IS1ATM(J+1)-1
          IF(J2.LT.J1) GO TO 90
C
C         SKIP I,J PAIR IF SUBSYSTEM LISTS HAVE NO CHANCE OF OVERLAP.
C
          IF(IASUBS(J2).LT.IASUBS(I1).OR.
     .       IASUBS(I2).LT.IASUBS(J1)) GO TO 90
C
C         LOOP OVER I-J SUBSYSTEMS AND LOOK FOR ONE IN COMMON.
C
          DO 80 ILIST=I1,I2
            ISUB = IASUBS(ILIST)
            DO 70 JLIST=J1,J2
              JSUB = IASUBS(JLIST)
              IF(ISUB.EQ.JSUB)THEN
C
C               WE HAVE A MATCH -- STORE THE I,J PAIR.
C
                NPAIRS = NPAIRS + 1
                IF(NPAIRS.GT.MBPAIR)THEN
                  IERROR = I
                  RETURN
                ENDIF
                IPAIR(NPAIRS) = J
                GO TO 90
              ENDIF
 70         CONTINUE
 80       CONTINUE
 90     CONTINUE
 100  CONTINUE
C
C     STORE EXTRA POINTER FOR CONVENIENCE IN LOOPING.
C
      IP1(NATOMS+1) = NPAIRS + 1
C-RDC       WRITE(IOUT,'(" NPAIRS = ",I8)') NPAIRS
C-RDC       IF(SCREEN) WRITE(ISCR,'(" NPAIRS = ",I8)') NPAIRS
      RETURN
      END