      subroutine xcfunctional(nfunctional,ixcfit,igrad,
     &                        dnstya,dnstyb,dnstyt,
     &                        gradxa,gradya,gradza,
     &                        gradxb,gradyb,gradzb,
     &                        gradxt,gradyt,gradzt,
     &                        hssxxa,hssxya,hssxza,hssyya,hssyza,hsszza,
     &                        hssxxb,hssxyb,hssxzb,hssyyb,hssyzb,hsszzb,
     &                        hssxxt,hssxyt,hssxzt,hssyyt,hssyzt,hsszzt,
     &                        xcnrgy,xcpota,xcpotb)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      data small,expcut/1.0d-15,80.0d0/

      data zero,one,two,three,four,five,six,seven,dnine,eleven/
     &    0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,5.0d0,6.0d0,7.0d0,9.0d0,11.0d0/

      data third,threequarters,fourthirds,fivethirds/
     &    0.33333333333333d0,0.75d0,1.3333333333333d0,1.6666666666667d0/

      data pi,dnom,zetadnom,pialpha/
     &    3.1415926535898d0,5.1297633d0,0.519842099d0,1.636963801d0/

      data fperdew,cp1,cp2,cp3,cp4/
     &    0.11d0,0.023266d0,7.389d-06,8.723d0,0.472d0/

      data tenthousand,thirtythousand,cpa,cpb,cpc,cpd/
     &    10000.0d0,30000.0d0,0.001667d0,0.002568d0,1.745d0,0.004235d0/

      ex=zero
      ec=zero

      vxa=zero
      vxb=zero

      vca=zero
      vcb=zero

      if(nfunctional.gt.0) call drho(gradxa,gradya,gradza,
     &                               hssxxa,hssxya,hssxza,
     &                               hssyya,hssyza,hsszza,
     &                               gmoda,dlapa,xcvta)

      if(nfunctional.gt.0) call drho(gradxb,gradyb,gradzb,
     &                               hssxxb,hssxyb,hssxzb,
     &                               hssyyb,hssyzb,hsszzb,
     &                               gmodb,dlapb,xcvtb)

      if(nfunctional.gt.0) call drho(gradxt,gradyt,gradzt,
     &                               hssxxt,hssxyt,hssxzt,
     &                               hssyyt,hssyzt,hsszzt,
     &                               gmodt,dlapt,xcvtt)

      if(dnstyt.gt.small) then
                                rs=one/(pi*fourthirds*dnstyt)**third

                                zeta=(dnstya-dnstyb)/dnstyt

                                if(zeta.gt.+one) zeta=+one
                                if(zeta.lt.-one) zeta=-one
                          else
                                xcnrgy=zero

                                xcpota=zero
                                xcpotb=zero

                                return
                          endif

      ex=-(threequarters/(pialpha*rs))*((one+zeta)**fourthirds
     &                                 +(one-zeta)**fourthirds)/two

      factor=-threequarters*(three/pi)**third

      vxa=factor*((two*dnstya)**third)*fourthirds
      vxb=factor*((two*dnstyb)**third)*fourthirds

      if(nfunctional.gt.0) then

            if((ixcfit.eq.1).or.(igrad.eq.1)) then

                       call beckex(dnstya,dnstyb,dnstyt,gmoda,gmodb,
     &                             dlapa,dlapb,xcvta,xcvtb,ex,vxa,vxb)

                                              else

                       ex=zero

                       vxa=zero
                       vxb=zero

                                              endif

                           endif

      sqrtrs=sqrt(rs)

      aa=vwnf1(1,sqrtrs)/two
      bb=vwnf2(1,sqrtrs)/two

      ec=aa

      vca=aa-bb*sqrtrs/six
      vcb=aa-bb*sqrtrs/six

      cc=(vwnf1(2,sqrtrs)/two)-aa
      dd=(vwnf2(2,sqrtrs)/two)-bb
      ee=(vwnf1(3,sqrtrs)/two)*three/dnom

      sf1=((one+zeta)**third-(one-zeta)**third)*dnom/two
      sf2=(((one+zeta)**fourthirds+(one-zeta)**fourthirds)-two)/zetadnom

      ff=sf2*(cc*zeta**4+ee*(one-zeta**4)-(sqrtrs/six)
     &  *(dd*zeta**4+(one-zeta**4)*vwnf2(3,sqrtrs)*three/two/dnom))

      gg=sf1*(cc*zeta**4+ee*(one-zeta**4))+four*sf2*((cc-ee)*zeta**3)

      ec=ec+sf2*(cc*(zeta**4)+ee*(one-zeta**4))

      vca=vca+ff+(one-zeta)*gg
      vcb=vcb+ff-(one+zeta)*gg

      if(nfunctional.eq.0) goto 1001
      if((ixcfit.eq.0).and.(igrad.eq.0)) goto 1001

      d2minus1=one/((two**third)*sqrt(((one+zeta)/two)**fivethirds
     &                               +((one-zeta)/two)**fivethirds))

      cnperdew=cpa+(cpb+cp1*rs+cp2*rs*rs)
     &            /(one+cp3*rs+cp4*rs*rs+tenthousand*cp2*rs*rs*rs)

      phi=zero

      if(dnstyt.gt.small) phi=cpc*fperdew*(cpd/cnperdew)*gmodt
     &                       /dnstyt**(seven/six)

      expfac=zero

      if(phi.lt.expcut) expfac=exp(-phi)

      if(dnstyt.gt.small) ec=ec+d2minus1*expfac*cnperdew*gmodt*gmodt
     &                         /(dnstyt**(seven/three))

      term=one+cp3*rs+cp4*rs*rs+tenthousand*cp2*rs*rs*rs

      term1=((cp1+two*cp2*rs)/term)-((cpb+cp1*rs+cp2*rs*rs)
     &      *(cp3+two*cp4*rs+thirtythousand*cp2*rs*rs)/term**2)

      term2=zero

      if(dnstyt.gt.small) term2=-(four*pi/dnine)
     &                         *((four*pi*dnstyt/three)**(-fourthirds))

      term=term1*term2

      perdewc=zero

      aa=zero
      bb=zero
      cc=zero

      if(dnstyt.gt.small) then
                                perdewc=d2minus1*expfac*gmodt*gmodt
     &                                 *(phi*phi-phi-one)*term
     &                                 /(dnstyt**fourthirds)

                                aa=d2minus1*expfac*cnperdew
     &                                            /(dnstyt**third)

                                bb=(two-phi)*dlapt/dnstyt

                                cc=(fourthirds-(eleven*phi/three)
     &                            +((seven/six)*phi*phi))*gmodt*gmodt
     &                            /(dnstyt*dnstyt)
                          endif

      rhograd=dnstyt*gmodt

      if(rhograd.gt.small) then
                                 rhograd=phi*(phi-three)*xcvtt/rhograd
                           else
                                 rhograd=zero
                           endif

      rhot13=dnstyt**third

      rhoa23=dnstya**(two/three)
      rhob23=dnstyb**(two/three)

      dd=six*(dnstyt**4)/(d2minus1*d2minus1)

      dda=zero
      ddb=zero

      if(dd.gt.(small**4)) then
                                 dda=five*rhot13*(rhoa23-rhob23)/dd
                                 ddb=five*rhot13*(rhob23-rhoa23)/dd
                           endif

      ee=(two**(two/three))*(two-phi)*dnstyt

      eea=ee*(gradxb*gradxt+gradyb*gradyt+gradzb*gradzt)
      eeb=ee*(gradxa*gradxt+gradya*gradyt+gradza*gradzt)

      factor=(two**(two/three))*(one-phi)*gmodt*gmodt

      vca=vca+perdewc-aa*(bb-cc+rhograd-dda*(factor*dnstyb-eea))
      vcb=vcb+perdewc-aa*(bb-cc+rhograd-ddb*(factor*dnstya-eeb))

 1001 continue

      xcnrgy=ex+ec

      xcpota=vxa+vca
      xcpotb=vxb+vcb

      return
      end
