C
      SUBROUTINE SETUP(N3TM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (N3TMMN = 6)
C
      COMMON/PARAM/R(3426),VALEN,BETA,C,D2,ALPHA2,R2,D,ALPHA,RE
     $,SAB,SAS,SBS,DH,ALPHAH,RH,ACONST,RCUT1,RCUT2,DCUT,EZERO
      COMMON/IPARAM/NMETAL,NATOM
C
      DIMENSION X(1000),Y(1000),Z(1000)
C
      IPARA = 4
      ILAT  = 7       
C
C CHECK THE NUMBER OF COORDINATES SET BY THE CALLING PROGRAM
C
       IF (N3TM .LT. N3TMMN) THEN
           WRITE (6, 1800) N3TM
           STOP 'SETUP 1'
       ENDIF
C
C OPEN THE POTENTIAL DATA FILES
C
       OPEN (UNIT=IPARA, FILE='poth2niads.dat', STATUS='OLD', 
     *       FORM='FORMATTED', ERR=100)
       OPEN (UNIT=ILAT, FILE='poth2nilat.dat', STATUS='OLD', 
     *       FORM='FORMATTED', ERR=100)
C
C INPUT PARAMATERS FOR THE H2-METAL POTENTIAL.
C
      READ(IPARA,500) D2,ALPHA2,R2
      READ(IPARA,500) D,ALPHA,RE
      READ(IPARA,500) DH,ALPHAH,RH
      READ(IPARA,500) SAB,SAS,SBS
      READ(IPARA,500) VALEN,BETA,C
      READ(IPARA,500) ACONST,RCUT1,RCUT2
      READ(IPARA,600) EZERO,NMETAL,NATOM
C
      DCUT = RCUT1 - RCUT2 
C
C READ IN THE COORDINATES OF THE METAL ATOMS.
C                     
      L = 0
      DO 10 I= 1,NMETAL
        READ(ILAT,700) X(I),Y(I),Z(I)
        L = L + 1
        R(L) = X(I)
        L = L + 1
        R(L) = Y(I)
        L = L + 1
        R(L) = Z(I)
  10  CONTINUE 
      WRITE(6,1000)
      WRITE(6,1100) D2,ALPHA2,R2
      WRITE(6,1200) D,ALPHA,RE
      WRITE(6,1300) DH,ALPHAH,RH
      WRITE(6,1400) SAB,SAS,SBS
      WRITE(6,1500) VALEN,BETA,C
      WRITE(6,1600) ACONST,RCUT1,RCUT2
      WRITE(6,1700) EZERO,NMETAL,NATOM
C
C  CLOSE THE POTENTIAL DATA FILES
C
      CLOSE (UNIT = 10)
      CLOSE (UNIT = 11)
C
      RETURN 
C
  500 FORMAT(3F15.7)
  600 FORMAT(F15.10,2I5)
  700 FORMAT(3F15.8)
1000  FORMAT(/,2X,T5,'SETUP has been called for the DePristo H2-Metal ',
     *               'potential',//,2X,T5,'Potential parameters:')
1100  FORMAT(2X,T5,'For the H-Metal two-body potential: ',
     *       /,2X,T10,'D2 = ',F10.7,' eV','  Alpha2 = ',F10.7,
     *                ' bohr-1','  R2 = ',F10.7,' bohr')
1200  FORMAT(2X,T5,'For H-H interaction:',
     *       /,2X,T10,'D = ',F10.7,' eV','  Alpha = ',F10.7,
     *                ' bohr-1','  Re = ',F10.7,' bohr')
1300  FORMAT(2X,T5,'For the embedded energy:',
     *       /,2X,T10,'DH = ',F10.7,' eV','  Alpha = ',F10.7,
     *                ' bohr-1','  RH = ',F10.7,' bohr')
1400  FORMAT(2X,T5,'Sato parameters:',
     *       /,2X,T10,'SAB = ',F10.7,'  SAS = ',F10.7,
     *                '  SBS = ',F10.7)
1500  FORMAT(2X,T5,'For the density function:',
     *       /,2X,T10,'number of valence electrons = ',F6.3,
     *                '  Beta = ',F10.7,'  C = ',F10.7)
1600  FORMAT(2X,T5,'Lattice constant = ',F10.5,
     *       /,2X,T5,'Cutoff radiuses:',
     *       /,2X,T10,'RCUT1 = ',F10.7,' bohr',
     *                ' RCUT2 = ',F10.7,' bohr')
1700  FORMAT(2X,T5,'Zero of energy = ',F10.7,
     *       /,2X,T5,'Total number of metal atoms = ',I5,
     *       /,2X,T5,'Total number of gas atoms = ',I5)
1800  FORMAT(2X,T5,'Warning: N3TM is set equal to ',I3,
     *               ' but this potential routine',
     *       /,2X,T14,'requires N3TM be greater than or ',
     *                'equal to ',I3,/)
C
  100 WRITE(6,*)'Error opening potential data file'
      STOP 'SETUP 2'
C
      END
C
C
      SUBROUTINE SURF(V, X, DVDX, N3TM)
C
C     THIS SUBROUTINE CALCULATES THE VALUE OF THE POTENTIAL AND ITS GRADIENT
C     AS A FUNCTION OF THE CARTESIAN COORDINATES FOR H2 ON NI METAL SURFACE 
C     (FCC 100) USING LEE AND DEPRISTO FUNCTIONAL FORMS.
C
C     ***** NOTE - THE PARAMETERS USED IN CALCULATING THE POTENTIAL ARE IN
C                  UNITS OF BOHR AND EV. THE POTENTIAL IS CONVERTED FROM EV
C                  TO AU IN THE FINAL STEP.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION JAS,JAB,JBS
C
      COMMON/PARAM/R(3426),VALEN,BETA,C,D2,ALPHA2,R2,D,ALPHA,RE
     $,SAB,SAS,SBS,DH,ALPHAH,RH,ACONST,RCUT1,RCUT2,DCUT,EZERO
      COMMON/IPARAM/NMETAL,NATOM
C 
      DIMENSION X(N3TM), DVDX(N3TM)
C
      DIMENSION DRAB(6),DJAB(6),DQAB(6),DVAS(3),DVBS(3)
     $,DRHOA(3),DRHOB(3),DRSA(3),DRSB(3),DASUM(6),DBSUM(3),DV2A(3)
     $,DV2B(3),DADIF(6),DBDIF(3),DQAS(3),DQBS(3),DJAS(3),DJBS(3)
     $,DAEEXP(3),DBEEXP(3),DVHOMA(3),DVHOMB(3),TX(6)
C
      DATA AUCAL / 3.675024898D-2 /, PI / 3.141592654D0 /
      THRD=1.0D0/3.0D0
C
C
C     CALCULATE RAB=DISTANCE BETWEEN H ATOMS A AND B AND ITS DERIVATIVES DRAB
C     SKIP THIS IN CASE OF AN ATOM ON SURFACE.
C
      IF (NATOM .GT. 1) THEN
        DX=X(4)-X(1)
        DY=X(5)-X(2)
        DZ=X(6)-X(3)
        RAB=SQRT(DX*DX+DY*DY+DZ*DZ)
        DRAB(1)=-DX/RAB
        DRAB(2)=-DY/RAB
        DRAB(3)=-DZ/RAB
        DRAB(4)=-DRAB(1)
        DRAB(5)=-DRAB(2)
        DRAB(6)=-DRAB(3)
      ENDIF
      DO 1 I=1,6
    1 TX(I)=X(I)
C
C     CALL LATEXP TO PERFORM ALL THE SUMMATIONS OVER THE METAL LATTICE ATOMS
C
      CALL LATEXP(RHOA,RHOB,DRHOA,DRHOB,ASUM,BSUM,DASUM,DBSUM,ADIF,BDIF
     +,DADIF,DBDIF,X,N3TM)
C
C     CALCULATE THE TOTAL ELECTRON DENSITIES, RHOA AND RHOB, AND RSA AND RSB
C     ALONG WITH THE DERIVATIVES, DRHOA,DRHOB,DRSA,DRSB
C     RSA=(3/(4*RHOA*PI))**1/3
C
C FOR ATOM A:
C
      RHOA=RHOA*C*VALEN
      DO 110 J=1,3
      DRHOA(J)=DRHOA(J)*C*VALEN
  110 DRHOB(J)=DRHOB(J)*C*VALEN
      RSA=(0.75D0/(PI*RHOA))**THRD
      DO 120 J=1,3
120   DRSA(J)=-THRD*RSA*DRHOA(J)/RHOA 
C
C FOR ATOM B:
C             
      IF (NATOM .EQ. 1) GOTO 160
C
      RHOB=RHOB*C*VALEN
      RSB=(0.75D0/(PI*RHOB))**THRD
      DO 150 J=1,3
150   DRSB(J)=-THRD*RSB*DRHOB(J)/RHOB
C
C     CALCULATE V2A AND THE ANTI-MORSE ANALOGS OF V2A AND ITS DERIVATIVES
C
  160 CONTINUE
      V2A=ADIF*D2
      ASUM=ASUM*D2/2.0D0
      DO 200 J=1,3
      DV2A(J)=DADIF(J)*D2
      DASUM(J)=DASUM(J)*D2/2.0D0  
  200 CONTINUE
C
C     CALCULATE V2B AND THE ANTI-MORSE ANALOGS OF V2B AND ITS DERIVATIVE
C     SKIP IN CASE OF AN ATOM ON SURFACE.
C
      IF (NATOM .GT. 1) THEN
        V2B=BDIF*D2
        BSUM=BSUM*D2/2.0D0
        DO 210 J=1,3
         DV2B(J)=DBDIF(J)*D2
         DBSUM(J)=DBSUM(J)*D2/2.0D0 
  210   CONTINUE  
      ENDIF

C
C     CALCULATE VAS,VBS ; VAS=VHOMA+V2A
C     VHOMA=DH<EXP(-2ALPHAH(RSA-RH))-2EXP(-ALPHAH(RSA-RH))>
C    
C FOR ATOM A:
C
      AEXP=-ALPHAH*(RSA-RH)
      AEXP=EXP(AEXP)
      VHOMA=DH*(AEXP*AEXP-2.0D0*AEXP) 
      DO 220 J=1,3
  220 DVHOMA(J)=-ALPHAH*(VHOMA+AEXP*DH)*2.0D0*DRSA(J)
      VAS=VHOMA+V2A
      DO 225 J=1,3
  225 DVAS(J)=DVHOMA(J)+DV2A(J)
C
C FOR ATOM B: RETURN IN THE CASE OF AN ATOM ON SURFACE
C                                                     
      IF (NATOM .EQ. 1) THEN
         V = VAS*AUCAL - EZERO
         DO 600 J =1,3
  600    DVDX(J) = DVAS(J)*AUCAL 
         RETURN
      ELSE
         BEXP=-ALPHAH*(RSB-RH)
         BEXP=EXP(BEXP)
         VHOMB=DH*(BEXP*BEXP-2.0D0*BEXP)
         DO 230 J=1,3
  230    DVHOMB(J)=-ALPHAH*(VHOMB+BEXP*DH)*2.0D0*DRSB(J)
         VBS=VHOMB+V2B
         DO 235 J=1,3
  235    DVBS(J)=DVHOMB(J)+DV2B(J) 
      ENDIF
C
C     CALCULATE QAS,QBS AND JAS,JBS
C     QAS+JAS=VAS
C     QAS-JAS=((1-SAS)/(1+SAS))<(DH/2)[EXP(-2ALPHAH(RSA-RH))+
C               2EXP(-ALPHAH(RSA-RH))+SUMA
C
      AEXP=-ALPHAH*(RSA-RH)
      BEXP=-ALPHAH*(RSB-RH)
      AEXP=EXP(AEXP)
      BEXP=EXP(BEXP)
      AEEXP=AEXP*AEXP+2.0D0*AEXP
      BEEXP=BEXP*BEXP+2.0D0*BEXP
      DO 320 J=1,3
      DAEEXP(J)=-ALPHAH*(AEEXP-AEXP)*2.0D0*DRSA(J)
  320 DBEEXP(J)=-ALPHAH*(BEEXP-BEXP)*2.0D0*DRSB(J)
      ADIF=((1.0D0-SAS)/(1.0D0+SAS))*((DH/2.0D0)*AEEXP+ASUM)
      BDIF=((1.0D0-SBS)/(1.0D0+SBS))*((DH/2.0D0)*BEEXP+BSUM)
      DO 330 J=1,3
      DADIF(J)=((1.0D0-SAS)/(1.0D0+SAS))*((DH/2.0D0)*DAEEXP(J)+DASUM(J))
  330 DBDIF(J)=((1.0D0-SBS)/(1.0D0+SBS))*((DH/2.0D0)*DBEEXP(J)+DBSUM(J))
      QAS=(VAS+ADIF)/2.0D0
      QBS=(VBS+BDIF)/2.0D0
      JAS=(VAS-ADIF)/2.0D0
      JBS=(VBS-BDIF)/2.0D0 
      DO 340 J=1,3
      DQAS(J)=(DVAS(J)+DADIF(J))/2.0D0
      DQBS(J)=(DVBS(J)+DBDIF(J))/2.0D0
      DJAS(J)=(DVAS(J)-DADIF(J))/2.0D0
  340 DJBS(J)=(DVBS(J)-DBDIF(J))/2.0D0
C
C     CALCULATE QAB AND JAB
C     QAB+JAB=D[EXP(-2ALPHA(RAB-RE)-2EXP(-ALPHA(RAB-RE))]
C     QAB-JAB=((1-SAB)/(1+SAB))(D/2)[EXP(-2ALPHA(RAB-RE))+2EXP(-ALPHA(RAB-RE))]
C
      SEXP=EXP(-ALPHA*(RAB-RE))
      SUM=D*(SEXP*SEXP-2.0D0*SEXP)
      FACT=(1.0D0-SAB)/(1.0D0+SAB)
      DIF=FACT*(SEXP*SEXP+2.0D0*SEXP)*D*0.5D0
      DO 400 J=1,6
      DASUM(J)=-ALPHA*2.0D0*(SUM+D*SEXP)*DRAB(J)
 400  DADIF(J)=-ALPHA*2.0D0*(DIF-D*FACT*0.5D0*SEXP)*DRAB(J)
      QAB=0.5D0*(SUM+DIF)
      JAB=0.5D0*(SUM-DIF)
      DO 450 J=1,6
      DQAB(J)=0.5D0*(DASUM(J)+DADIF(J))
  450 DJAB(J)=0.5D0*(DASUM(J)-DADIF(J))
C
C     CALCULATE V=QAS+QBS+QAB-SQRT[JAB(JAB-JAS-JBS)+(JAS+JBS)**2]
C
      PART=SQRT(JAB*(JAB-JAS-JBS)+(JAS+JBS)*(JAS+JBS))
      V=QAS+QBS+QAB-PART
      V=(V+D)*AUCAL-EZERO
      DO 500 J=1,3
      DVDX(J)=DQAS(J)+DQAB(J)-0.5D0*(DJAB(J)*(JAB-JAS-JBS)+JAB*(
     $DJAB(J)-DJAS(J))+2.0D0*(JAS+JBS)*DJAS(J))/PART
      DVDX(J)=DVDX(J)*AUCAL
      DVDX(J+3)=DQBS(J)+DQAB(J+3)-0.5D0*(DJAB(J+3)*(JAB-JAS-JBS)
     $+JAB*(DJAB(J+3)-DJBS(J))+2.0D0*(JAS+JBS)*DJBS(J))/PART
      DVDX(J+3)=DVDX(J+ 3)*AUCAL
  500 CONTINUE
      RETURN
      END
                                       
C                    
      SUBROUTINE LATEXP(RHOA,RHOB,DRHOA,DRHOB,ASUM,BSUM,DASUM,DBSUM
     $,ADIF,BDIF,DADIF,DBDIF,X,N3TM)
C
C     THIS SUBROUTINE CALCULATES THE DISTANCES BETWEEN THE ADATOMS AND
C     METAL ATOMS OF THE LATTICE, RAS AND RBS, AND THEIR DERIVATIVES.
C     THE ADATOM COORDINATES ARE CONVERTED TO WITHIN THE UNIT CELL CENTERED
C     AT THE ORIGIN.  
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON/INTERN/RAS(1142),RBS(1142),DRAS(3,1142),DRBS(3,1142)
     $,INDEXA(800),INDEXB(800),NCOUNA,NCOUNB 
C
      COMMON/PARAM/R(3426),VALEN,BETA,C,D2,ALPHA2,R2,D,ALPHA,RE
     $,SAB,SAS,SBS,DH,ALPHAH,RH,ACONST,RCUT1,RCUT2,DCUT,EZERO
      COMMON/IPARAM/NMETAL,NATOM
C
      DIMENSION X(N3TM)
      DIMENSION DPRA1(3),DNRA1(3),DPRB1(3),DNRB1(3),DPDA1(3),DNDA1(3)
     +,DPDB1(3),DNDB1(3),DPSA1(3),DNSA1(3),DPSB1(3),DNSB1(3)
      DIMENSION DRHOA(3),DRHOB(3),DADIF(3),DBDIF(3),DASUM(3),DBSUM(3)
C
C
C     LOOP OVER ALL METAL ATOMS
C
      NCOUNA=0
      NCOUNB=0
      DO 100 I=1,NMETAL      
C
C     FIND THE DISTANCES BETWEEN THE ADATOMS AND THE METAL ATOM, RAS AND RBS
C     AND THEIR DERIVATIVES, DRAS AND DRBS. SKIP B'S IN CASE OF ATOM ON
C     SURFACE.
C
      JSTRT=1+(I-1)*3
      JEND=JSTRT+2
      RAS(I)=0.0D0
      RBS(I)=0.0D0
      K=0
      DO 70 J=JSTRT,JEND
       K=K+1
       DELAX=R(J)-X(K)
       RAS(I)=RAS(I)+DELAX*DELAX
       DRAS(K,I)=-DELAX
       IF (NATOM .EQ.1 ) GOTO 70
       DELBX=R(J)-X(K+3)
       RBS(I)=RBS(I)+DELBX*DELBX
       DRBS(K,I)=-DELBX
   70 CONTINUE 
      IF (RAS(I) .LE. RCUT2*RCUT2) THEN
         NCOUNA = NCOUNA + 1
         INDEXA(NCOUNA) = I
         RAS(I)=SQRT(RAS(I)) 
      ENDIF
      IF (RBS(I) .LE. RCUT2*RCUT2) THEN
         NCOUNB = NCOUNB + 1
         INDEXB(NCOUNB) = I
         RBS(I)=SQRT(RBS(I)) 
      ENDIF 
      DO 80 J=1,3
                       DRAS(J,I)=DRAS(J,I)/RAS(I)
      IF (NATOM .GT.1) DRBS(J,I)=DRBS(J,I)/RBS(I)   
   80 CONTINUE    
  100 CONTINUE
C
C     CALL LATSUM TO PERFORM THE SUMS OVER ALL METAL ATOMS
C
      CALL LATSUM(X,RA1,RB1,DPRA1,DPRB1,SA1,SB1,DPSA1,DPSB1,DA1
     +  ,DB1,DPDA1,DPDB1,DNRA1,DNRB1,DNSA1,DNSB1,DNDA1,DNDB1)
C
      RHOA=RA1
      RHOB=RB1
      ASUM=SA1
      BSUM=SB1
      ADIF=DA1
      BDIF=DB1
C
C     EVALUATE THE X,Y,Z DERIVATIVES
C
      DO 350 I=1,3
      TEMP1=DPRA1(I)
      TEMP2=DNRA1(I)
      DRHOA(I)=TEMP1+TEMP2
C
      TEMP1=DPSA1(I)
      TEMP2=DNSA1(I)
      DASUM(I)=TEMP1+TEMP2
C
      TEMP1=DPDA1(I)
      TEMP2=DNDA1(I)
      DADIF(I)=TEMP1+TEMP2
C    
      IF (NATOM .EQ. 1) GOTO 350
      TEMP1=DPRB1(I)
      TEMP2=DNRB1(I)
      DRHOB(I)=TEMP1+TEMP2
C
      TEMP1=DPSB1(I)
      TEMP2=DNSB1(I)
      DBSUM(I)=TEMP1+TEMP2
C
      TEMP1=DPDB1(I)
      TEMP2=DNDB1(I)
      DBDIF(I)=TEMP1+TEMP2
  350 CONTINUE
C
      RETURN         
      END
                                       
C
      SUBROUTINE LATSUM(X,RHOA,RHOB,DPRHOA,DPRHOB,ASUM,BSUM
     +,DPASUM,DPBSUM,ADIF,BDIF,DPADIF,DPBDIF,DNRHOA,DNRHOB,DNASUM
     +,DNBSUM,DNADIF,DNBDIF)
C
C     THIS SUBROUTINE CALCULATES THE VALUES WHICH INVOLVE SUMS OVER LATTICE 
C     METAL ATOMS: RHOA,RHOB,ASUM,BSUM,ADIF, AND BDIF, AND THEIR DERIVATIVES. 
C     THE DERIVATIVES ARE HANDLED AS SUMS OF THEIR NEGATIVE AND POSITIVE PARTS.
C     A "SOFTENED" CUT-OFF RADIUS IS USED.  THE EXPONENTIAL FUNCTIONS ARE 
C     MULTIPLIED BY 
C
C     ((R-R2)/(R1-R2)) - SIN[ PI*(2R-R1-R2)/(R2-R1)]/(2*PI)
C
C     WHEN RCUT1<RAS<RCUT2.
C     FOR RAS>RCUT2, THE CONTRIBUTION FROM THE METAL ATOM IS NOT INCLUDED IN THE
C     SUM.
C
C     RHOA=SUM(S) {EXP[-BETA*RAS]}
C     ASUM=SUM(S) {EXP[-2ALPHA2(RAS-R2)] + 2*EXP[-ALPHA2(RAS-R2)]}
C     ADIF=SUM(S) {EXP[-2ALPHA2(RAS-R2)] - 2*EXP[-ALPHA2(RAS-R2)]}
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C
      COMMON/INTERN/RAS(1142),RBS(1142),DRAS(3,1142),DRBS(3,1142)
     $,INDEXA(800),INDEXB(800),NCOUNA,NCOUNB
C
      COMMON/PARAM/R(3426),VALEN,BETA,C,D2,ALPHA2,R2,D,ALPHA,RE
     $,SAB,SAS,SBS,DH,ALPHAH,RH,ACONST,RCUT1,RCUT2,DCUT,EZERO
      COMMON/IPARAM/NMETAL,NATOM
C
      DIMENSION X(6),DPRHOA(3),DNRHOA(3),DPRHOB(3),DNRHOB(3),DPADIF(3)
     +,DNADIF(3),DPBDIF(3),DNBDIF(3),DPASUM(3),DNASUM(3),DPBSUM(3)
     +,DNBSUM(3)
C
      DATA PI / 3.141592654D0 /
C
      RHOA=0.0D0
      RHOB=0.0D0
      ASUM=0.0D0
      BSUM=0.0D0
      ADIF=0.0D0
      BDIF=0.0D0 
      DO 10 J=1,3
      DPRHOA(J)=0.0D0
      DPRHOB(J)=0.0D0
      DPASUM(J)=0.0D0
      DPBSUM(J)=0.0D0
      DPADIF(J)=0.0D0
      DPBDIF(J)=0.0D0
      DNRHOA(J)=0.0D0
      DNRHOB(J)=0.0D0
      DNASUM(J)=0.0D0
      DNBSUM(J)=0.0D0
      DNADIF(J)=0.0D0
      DNBDIF(J)=0.0D0
   10 CONTINUE
C
C     LOOP OVER ALL METAL ATOMS
C
C
C     FOR ADATOM A:
C
      DO 280 L=1,NCOUNA     
       I=INDEXA(L)                  
      IF(RAS(I).GT.RCUT1) THEN
        RR2=RAS(I)-RCUT2
        RR1=RAS(I)-RCUT1
        FACTA=(RR2/DCUT)+(SIN( PI*(RR1+RR2)/DCUT))/(2.0D0*PI)
        DFACTA=(1.0D0+COS( PI*(RR1+RR2)/DCUT))/DCUT
      ELSE IF(RAS(I).LE.RCUT1) THEN
        FACTA=1.0D0
        DFACTA=0.0D0
      END IF
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO RHOA 
C
      PARTA=EXP(-BETA*RAS(I))
      DO 40 J=1,3
      DPRTA=PARTA*DRAS(J,I)*(DFACTA-FACTA*BETA)
      IF(DPRTA.GT.0.0D0) DPRHOA(J)=DPRHOA(J)+DPRTA
      IF(DPRTA.LT.0.0D0) DNRHOA(J)=DNRHOA(J)+DPRTA
   40 CONTINUE
      RHOA=RHOA+PARTA*FACTA
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO DIFA AND DIFB
C
      AEXP=-ALPHA2*(RAS(I)-R2)
      AEXP=EXP(AEXP)
      PARTA=(AEXP*AEXP-2.0D0*AEXP)
      DO 180 J=1,3
      DPRTA=-(PARTA+AEXP)*DRAS(J,I)*ALPHA2*FACTA*2.0D0
      IF(DPRTA.GT.0.0D0) DPADIF(J)=DPADIF(J)+DPRTA
      IF(DPRTA.LT.0.0D0) DNADIF(J)=DNADIF(J)+DPRTA
      IF(DFACTA.EQ.0.0D0) GO TO 180
      DPRTA=PARTA*DRAS(J,I)*DFACTA
      IF(DPRTA.GT.0.0D0) DPADIF(J)=DPADIF(J)+DPRTA
      IF(DPRTA.LT.0.0D0) DNADIF(J)=DNADIF(J)+DPRTA
 180  CONTINUE
      ADIF=ADIF+PARTA*FACTA
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO ASUM 
C
      PARTA=(AEXP*AEXP+2.0D0*AEXP)
      DO 250 J=1,3
      DPRTA=-(PARTA-AEXP)*DRAS(J,I)*FACTA*ALPHA2*2.0D0
      IF(DPRTA.GT.0.0D0) DPASUM(J)=DPASUM(J)+DPRTA
      IF(DPRTA.LT.0.0D0) DNASUM(J)=DNASUM(J)+DPRTA
      IF(DFACTA.EQ.0.0D0) GO TO 250
      DPRTA=PARTA*DRAS(J,I)*DFACTA
      IF(DPRTA.GT.0.0D0) DPASUM(J)=DPASUM(J)+DPRTA
      IF(DPRTA.LT.0.0D0) DNASUM(J)=DNASUM(J)+DPRTA
  250 CONTINUE
      ASUM=ASUM+PARTA*FACTA   
  280 CONTINUE
C
      IF (NATOM .EQ. 1) RETURN  
C
C     FOR ADATOM B:
C
      DO 500 K =1,NCOUNB
        I=INDEXB(K)                
      IF(RBS(I).GT.RCUT1) THEN
        RR2=RBS(I)-RCUT2
        RR1=RBS(I)-RCUT1
        FACTB=(RR2/DCUT)+(SIN( PI*(RR1+RR2)/DCUT))/(2.0D0*PI)
        DFACTB=(1.0D0+COS( PI*(RR1+RR2)/DCUT))/DCUT
      ELSE IF(RBS(I).LE.RCUT1) THEN
        FACTB=1.0D0
        DFACTB=0.0D0
      END IF
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO RHOB
C
      PARTB=EXP(-BETA*RBS(I))
      DO 340 J=1,3
      DPRTB=PARTB*DRBS(J,I)*(DFACTB-FACTB*BETA)
      IF(DPRTB.GT.0.0D0) DPRHOB(J)=DPRHOB(J)+DPRTB
      IF(DPRTB.LT.0.0D0) DNRHOB(J)=DNRHOB(J)+DPRTB
  340 CONTINUE
      RHOB=RHOB+PARTB*FACTB
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO DIFB
C
      BEXP=-ALPHA2*(RBS(I)-R2)
      BEXP=EXP(BEXP)
      PARTB=(BEXP*BEXP-2.0D0*BEXP)
      DO 380 J=1,3
      DPRTB=-(PARTB+BEXP)*DRBS(J,I)*ALPHA2*FACTB*2.0D0
      IF(DPRTB.GT.0.0D0) DPBDIF(J)=DPBDIF(J)+DPRTB
      IF(DPRTB.LT.0.0D0) DNBDIF(J)=DNBDIF(J)+DPRTB
      IF(DFACTB.EQ.0.0D0) GO TO 380
      DPRTB=PARTB*DRBS(J,I)*DFACTB
      IF(DPRTB.GT.0.0D0) DPBDIF(J)=DPBDIF(J)+DPRTB
      IF(DPRTB.LT.0.0D0) DNBDIF(J)=DNBDIF(J)+DPRTB
 380  CONTINUE
      BDIF=BDIF+PARTB*FACTB
C
C     FIND THE METAL ATOM'S CONTRIBUTION TO BSUM
C
      PARTB=(BEXP*BEXP+2.0D0*BEXP)
      DO 450 J=1,3
      DPRTB=-(PARTB-BEXP)*DRBS(J,I)*FACTB*ALPHA2*2.0D0
      IF(DPRTB.GT.0.0D0) DPBSUM(J)=DPBSUM(J)+DPRTB
      IF(DPRTB.LT.0.0D0) DNBSUM(J)=DNBSUM(J)+DPRTB
      IF(DFACTB.EQ.0.0D0) GO TO 450
      DPRTB=PARTB*DRBS(J,I)*DFACTB
      IF(DPRTB.GT.0.0D0) DPBSUM(J)=DPBSUM(J)+DPRTB
      IF(DPRTB.LT.0.0D0) DNBSUM(J)=DNBSUM(J)+DPRTB
  450 CONTINUE
      BSUM=BSUM+PARTB*FACTB
  500 CONTINUE
C
      RETURN                    
      END