C
C $Id: diat.F,v 1.3 1998/07/16 16:39:36 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE GDIAT(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,AI0,AI1,AI2,DI1,
     .                DI2,BSI,BPI,BDI,ACI,AGI,BGI,CGI,XI,YI,ZI,ZNUCI,
     .                IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,AJ0,AJ1,AJ2,DJ1,
     .                DJ2,BSJ,BPJ,BDJ,ACJ,AGJ,BGJ,CGJ,XJ,YJ,ZJ,ZNUCJ,
     .               RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,DXREPIJ,DYREPIJ,DZREPIJ,
     .                                      DXHIJ,DYHIJ,DZHIJ,
     .                                      DXENNIJ,DYENNIJ,DZENNIJ,     
     .                EENUC,resinter)
c    .                RIJ,ILEVEL,ENNIJ,HIJ,REPIJ,EENUC)
c resinter = 1 if interresidue, 0 if intraresidue interaction between I and J
C
C     ILEVEL CODES:  0 ==> DO ENNIJ, HIJ, REPIJ, EENUC
C                    1 ==> DO ENNIJ,      REPIJ, EENUC
C                    2 ==> DO ENNIJ, HIJ, REPIJ
C                    3 ==> DO ENNIJ,      REPIJ
C                    4 ==> DO             REPIJ
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
#include "constants.h"
      real*8 AGI(4),BGI(4),CGI(4),AGJ(4),BGJ(4),CGJ(4)
      real*8 HIJ(9,9),REPIJ(10,10),EENUC(2,10)
      real*8  DXREPIJ(10,10),DYREPIJ(10,10),DZREPIJ(10,10)
      real*8  DXHIJ(9,9),DYHIJ(9,9),DZHIJ(9,9)
      integer resinter
      logical lok
C
C     LOCAL ARRAYS:
C
      LOGICAL ORBSI,ORBSJ
      real*8 SIJ(9,9)
      integer IREPUL(10)
      real*8 DXSIJ(9,9),DYSIJ(9,9),DZSIJ(9,9)
      DATA IREPUL /1,2,5,3,6,8,4,7,9,10/
      SAVE IREPUL
C

         do i=1,10
          do j=1,10
           REPIJ(j,i) = 0.0d0
           DXREPIJ(j,i) = 0.0d0
           DYREPIJ(j,i) = 0.0d0
           DZREPIJ(j,i) = 0.0d0
          enddo
         enddo

         do i=1,9
          do j=1,9
           HIJ(j,i) = 0.0d0
           DXHIJ(j,i) = 0.0d0
           DYHIJ(j,i) = 0.0d0
           DZHIJ(j,i) = 0.0d0
           SIJ(j,i) = 0.0d0
           DXSIJ(j,i) = 0.0d0
           DYSIJ(j,i) = 0.0d0
           DZSIJ(j,i) = 0.0d0
          enddo
         enddo
          ENNIJ = 0.0d0
          DXENNIJ = 0.0d0
          DYENNIJ = 0.0d0
          DZENNIJ = 0.0d0

C
      ORBSI = NORBSI.GT.0
      ORBSJ = NORBSJ.GT.0
C
C     GET DIATOMIC OVERLAPS AND ASSIGN OFF-DIAGONAL BLOCK OF
C     1-ELECTRON MATRIX (NON-SPARKLE ATOMS ONLY).
C
      IF(ORBSI.AND.ORBSJ.AND.(ILEVEL.EQ.0.OR.ILEVEL.EQ.2))THEN
        CALL GOVERLP(IAI,NORBSI,NPQI,XISI,XIPI,XIDI,XI,YI,ZI,
     .              IAJ,NORBSJ,NPQJ,XISJ,XIPJ,XIDJ,XJ,YJ,ZJ,
     .              RIJ,SIJ,DXSIJ,DYSIJ,DZSIJ)
        BETAI = BSI
        DO 20 I=1,NORBSI
          IF(I.EQ.2)THEN
            BETAI = BPI
          ELSEIF(I.EQ.5)THEN
            BETAI = BDI
          ENDIF
          BETAJ = BSJ
          DO 10 J=1,NORBSJ
            IF(J.EQ.2)THEN
              BETAJ = BPJ
            ELSEIF(J.EQ.5)THEN
              BETAJ = BDJ
            ENDIF
            HIJ(I,J) = 0.5D0*(BETAI+BETAJ)*SIJ(I,J)
             DXHIJ(I,J) = 0.5D0*(BETAI+BETAJ)*DXSIJ(I,J)
             DYHIJ(I,J) = 0.5D0*(BETAI+BETAJ)*DYSIJ(I,J)
             DZHIJ(I,J) = 0.5D0*(BETAI+BETAJ)*DZSIJ(I,J)            
 10       CONTINUE
 20     CONTINUE
      ENDIF
C
C     GET ELECTRON-ELECTRON REPULSIONS.  TREAT SPARKLES AS MONOPOLES.
C
      NIMAX = MAX(NORBSI,1)
      NJMAX = MAX(NORBSJ,1)
      CALL GREPUL(IAI,NIMAX,AI0,AI1,AI2,DI1,DI2,XI,YI,ZI,
     .           IAJ,NJMAX,AJ0,AJ1,AJ2,DJ1,DJ2,XJ,YJ,ZJ,
     .           RIJ,REPIJ,DXREPIJ,DYREPIJ,DZREPIJ)
C
C     WE ARE DONE IF ONLY REPIJ HAS BEEN REQUESTED.
C
      IF(ILEVEL.EQ.4) RETURN
C
C     ASSEMBLE ELECTRON-NUCLEAR ATTRACTION INTEGRALS.
C
C     ELECTRONS ON I ATTRACTED TO CORE OF J:
C
      IF(ILEVEL.EQ.0.OR.ILEVEL.EQ.1)THEN
        NI = (NORBSI*(NORBSI+1))/2
        IF(NI.GT.0)THEN
          DO 130 K=1,NI
            EENUC(1,K) = -ZNUCJ*REPIJ(IREPUL(K),1)
 130      CONTINUE
        ENDIF
C
        NJ = (NORBSJ*(NORBSJ+1))/2
C
C       ELECTRONS ON J ATTRACTED TO CORE OF I:
C
        IF(NJ.GT.0)THEN
          DO 140 K=1,NJ
            EENUC(2,K) = -ZNUCI*REPIJ(1,IREPUL(K))
 140      CONTINUE
        ENDIF
      ENDIF
C
C     CORE-CORE REPULSION:
C

   

      IAMAX = MAX(IAI,IAJ)
      IAMIN = MIN(IAI,IAJ)
      ZIJ = ZNUCI*ZNUCJ
      REPSS = REPIJ(1,1)
       DXREPSS = DXREPIJ(1,1)*Bohr2Ang
       DYREPSS = DYREPIJ(1,1)*Bohr2Ang
       DZREPSS = DZREPIJ(1,1)*Bohr2Ang
      EIJ = ZIJ*REPSS
       DXEIJ = ZIJ*DXREPSS
       DYEIJ = ZIJ*DYREPSS
       DZEIJ = ZIJ*DZREPSS

      ARGI = ACI*RIJ
      IF(ARGI.GT.25.0D0)THEN
        EXPI = 0.0D0
        DEXPI = 0.0D0
      ELSE
        EXPI = EXP(-ARGI)
        DEXPI = -ACI*EXP(-ARGI)
      ENDIF
      ARGJ = ACJ*RIJ
      IF(ARGJ.GT.25.0D0)THEN
        EXPJ = 0.0D0
         DEXPJ = 0.0D0
      ELSE
        EXPJ = EXP(-ARGJ)
         DEXPJ = -ACJ*EXP(-ARGJ)
      ENDIF
      IF(IAMIN.EQ.1.AND.(IAMAX.EQ.7.OR.IAMAX.EQ.8))THEN
C
C       N-H OR O-H INTERACTIONS:
C
        IF(IAI.EQ.1)THEN
          HTERM = EXPI
           DHTERM = DEXPI
          XTERM = EXPJ
           DXTERM = DEXPJ
        ELSE
          HTERM = EXPJ
           DHTERM = DEXPJ
          XTERM = EXPI
           DXTERM = DEXPI
        ENDIF
        ENNIJ = EIJ*(1.0D0 + RIJ*XTERM + HTERM)
        DXENNIJ = DXEIJ*(1.0D0 + RIJ*XTERM + HTERM)
     &    + EIJ*(RIJ*DXTERM + XTERM + DHTERM)*(XI-XJ)/RIJ*Bohr2Ang
        DYENNIJ = DYEIJ*(1.0D0 + RIJ*XTERM + HTERM)
     &    + EIJ*(RIJ*DXTERM + XTERM + DHTERM)*(YI-YJ)/RIJ*Bohr2Ang
        DZENNIJ = DZEIJ*(1.0D0 + RIJ*XTERM + HTERM)
     &    + EIJ*(RIJ*DXTERM + XTERM + DHTERM)*(ZI-ZJ)/RIJ*Bohr2Ang
      ELSE
C
C       REGULAR INTERACTION:
C
        ENNIJ = EIJ*(1.0D0 + EXPI + EXPJ)
        DXENNIJ = DXEIJ*(1.0D0 + EXPI + EXPJ)
     &            + EIJ*(DEXPI + DEXPJ)*(XI-XJ)/RIJ*Bohr2Ang
        DYENNIJ = DYEIJ*(1.0D0 + EXPI + EXPJ)
     &            + EIJ*(DEXPI + DEXPJ)*(YI-YJ)/RIJ*Bohr2Ang
        DZENNIJ = DZEIJ*(1.0D0 + EXPI + EXPJ)
     &            + EIJ*(DEXPI + DEXPJ)*(ZI-ZJ)/RIJ*Bohr2Ang
      ENDIF
C
C     ADD GAUSSIANS FOR AM1 OR PM3 (NON-SPARKLE ATOMS ONLY).
C
      IF(AM1)THEN
        GTERM = 0.0D0
        DGTERM = 0.0D0
        IF(ORBSI)THEN
          DO 150 IGAUS=1,4
            ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2
            IF(ARGI.LT.25.0D0)THEN
              GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI)
         DGTERM = DGTERM - AGI(IGAUS)*2.0d0*BGI(IGAUS)*(RIJ-CGI(IGAUS))
     .                                                      *EXP(-ARGI)
            ENDIF
 150      CONTINUE
        ENDIF
        IF(ORBSJ)THEN
          DO 160 JGAUS=1,4
            ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2
            IF(ARGJ.LT.25.0D0)THEN
              GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ)
         DGTERM = DGTERM - AGJ(JGAUS)*2.0d0*BGJ(JGAUS)*(RIJ-CGJ(JGAUS))
     .                                                      *EXP(-ARGJ)
            ENDIF
 160      CONTINUE
        ENDIF
        ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
         DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
         DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
         DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
      ELSEIF(RM1)THEN
        GTERM = 0.0D0
        DGTERM = 0.0D0
        IF(ORBSI)THEN
          DO 151 IGAUS=1,4
            ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2
            IF(ARGI.LT.25.0D0)THEN
              GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI)
         DGTERM = DGTERM - AGI(IGAUS)*2.0d0*BGI(IGAUS)*(RIJ-CGI(IGAUS))
     .                                                      *EXP(-ARGI)
            ENDIF
 151      CONTINUE
        ENDIF
        IF(ORBSJ)THEN
          DO 161 JGAUS=1,4
            ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2
            IF(ARGJ.LT.25.0D0)THEN
              GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ)
         DGTERM = DGTERM - AGJ(JGAUS)*2.0d0*BGJ(JGAUS)*(RIJ-CGJ(JGAUS))
     .                                                      *EXP(-ARGJ)
            ENDIF
 161      CONTINUE
        ENDIF
        ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
         DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
         DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
         DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
      ELSEIF(AM1D)THEN
        GTERM = 0.0D0
        DGTERM = 0.0D0
        IF(ORBSI)THEN
          DO 152 IGAUS=1,4
            ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2
            IF(ARGI.LT.25.0D0)THEN
              GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI)
         DGTERM = DGTERM - AGI(IGAUS)*2.0d0*BGI(IGAUS)*(RIJ-CGI(IGAUS))
     .                                                      *EXP(-ARGI)
            ENDIF
 152      CONTINUE
        ENDIF
        IF(ORBSJ)THEN
          DO 162 JGAUS=1,4
            ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2
            IF(ARGJ.LT.25.0D0)THEN
              GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ)
         DGTERM = DGTERM - AGJ(JGAUS)*2.0d0*BGJ(JGAUS)*(RIJ-CGJ(JGAUS))
     .                                                      *EXP(-ARGJ)
            ENDIF
 162      CONTINUE
        ENDIF
        ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
         DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
         DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
         DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
      ELSEIF(PM3PDDG)THEN
        GTERM = 0.0D0
        DGTERM = 0.0D0
        IF(ORBSI)THEN
          DO IGAUS=1,4
            ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2
            IF(ARGI.LT.25.0D0)THEN
              GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI)
         DGTERM = DGTERM - AGI(IGAUS)*2.0d0*BGI(IGAUS)*(RIJ-CGI(IGAUS))
     .                                                      *EXP(-ARGI)
            ENDIF
          END DO
        ENDIF
        IF(ORBSJ)THEN
          DO JGAUS=1,4
            ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2
            IF(ARGJ.LT.25.0D0)THEN
              GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ)
         DGTERM = DGTERM - AGJ(JGAUS)*2.0d0*BGJ(JGAUS)*(RIJ-CGJ(JGAUS))
     .                                                      *EXP(-ARGJ)
            ENDIF
          END DO
        ENDIF
        ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
         DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
         DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
         DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
        DNAPNB = 1.0d0/(ZNUCI+ZNUCJ)
        DO IGAUS = 1,2
          DO JGAUS = 1,2
             TERM1 = ZNUCI*PA7(IGAUS,IAI) + ZNUCJ*PA7(JGAUS,IAJ)
             TERM2 = RIJ - DA7(IGAUS,IAI) - DA7(JGAUS,IAJ)
             TERM3 = DNAPNB*TERM1*EXP(-10.0d0*TERM2*TERM2)
             ENNIJ = ENNIJ + TERM3
             DXENNIJ = DXENNIJ - 20.0d0*TERM2*TERM3*(XI-XJ)/RIJ*Bohr2Ang
             DYENNIJ = DYENNIJ - 20.0d0*TERM2*TERM3*(YI-YJ)/RIJ*Bohr2Ang
             DZENNIJ = DZENNIJ - 20.0d0*TERM2*TERM3*(ZI-ZJ)/RIJ*Bohr2Ang
          END DO
        END DO
      ELSEIF(PM3)THEN
        GTERM = 0.0D0
         DGTERM = 0.0D0
        IF(ORBSI)THEN
          DO 180 IGAUS=1,2
            ARGI = BGI(IGAUS)*(RIJ-CGI(IGAUS))**2
            IF(ARGI.LT.25.0D0)THEN
              GTERM = GTERM + AGI(IGAUS)*EXP(-ARGI)
               DGTERM = DGTERM - 2.0d0*BGI(IGAUS)*(RIJ-CGI(IGAUS))
     .                                *AGI(IGAUS)*EXP(-ARGI)
            ENDIF
 180      CONTINUE
        ENDIF
        IF(ORBSJ)THEN
          DO 200 JGAUS=1,2
            ARGJ = BGJ(JGAUS)*(RIJ-CGJ(JGAUS))**2
            IF(ARGJ.LT.25.0D0)THEN
              GTERM = GTERM + AGJ(JGAUS)*EXP(-ARGJ)
               DGTERM = DGTERM - 2.0d0*BGJ(JGAUS)*(RIJ-CGJ(JGAUS))
     .                                *AGJ(JGAUS)*EXP(-ARGJ)
            ENDIF
 200      CONTINUE
        ENDIF
c
c *** Modification : GM - 01 fev 2000
c *** Change       : divcon.1.0.C011
c *** Object       : introduce PIF 
c                    (PM3 Parametrized Intermolecular Function for water)
c *** Modified     : GM - 13 mar 2010
c                    (additional PIF parameters)
c
       IF (PIF) THEN
        if (resinter.eq.0) then	
           ! intramolecular bond
           ! use standard PM3 term
           ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
           DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
           DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
           DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang

        elseif (resinter.eq.2) then
           call gwalid(iai,iaj,rij,ewalid,dewalid)
           ENNIJ = ENNIJ +  ewalid
           DXENNIJ = DXENNIJ + dewalid*(XI-XJ)/RIJ
           DYENNIJ = DYENNIJ + dewalid*(YI-YJ)/RIJ
           DZENNIJ = DZENNIJ + dewalid*(ZI-ZJ)/RIJ
        elseif (resinter.eq.1) then
            call gpif(iai,iaj,rij,epif,depif,lok)
            if (lok) then
              ENNIJ = ENNIJ +  epif
              DXENNIJ = DXENNIJ + depif*(XI-XJ)/RIJ
              DYENNIJ = DYENNIJ + depif*(YI-YJ)/RIJ
              DZENNIJ = DZENNIJ + depif*(ZI-ZJ)/RIJ
            else
              write(*,'("WARNING: Use standard interaction between",i5,
     &                 " and",i5)') iai, iaj
              ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
              DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
              DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
              DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
           endif
        else
            continue
        endif
       ELSEIF (mais) THEN        
            call gmais(iai,iaj,rij,emais,demais)
            ENNIJ = ENNIJ +  emais
            DXENNIJ = DXENNIJ + demais*(XI-XJ)/RIJ
            DYENNIJ = DYENNIJ + demais*(YI-YJ)/RIJ
            DZENNIJ = DZENNIJ + demais*(ZI-ZJ)/RIJ
       ELSE
           ENNIJ = ENNIJ + GTERM*ZIJ/RIJ
           DXENNIJ =  DXENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(XI-XJ)/RIJ*Bohr2Ang
           DYENNIJ =  DYENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(YI-YJ)/RIJ*Bohr2Ang
           DZENNIJ =  DZENNIJ + (DGTERM*ZIJ/RIJ 
     .                    - GTERM*ZIJ/RIJ**2)*(ZI-ZJ)/RIJ*Bohr2Ang
       ENDIF    
      ENDIF
      RETURN
      END