      subroutine efids(m,n,icontr,inddyy,esp,efix,efiy,efiz,ccx,ccy,
     &                 ccz,g,density,zeta,px,py,pz,pxax,pyay,pzaz,
     &                 base,ss0,ss1,ss2,ps0,ps1,ds0,dscond,c,pxcx,
     &                 pycy,pzcz,u,twozeta,cx,cy,cz,ss3,saxs0,says0,
     &                 sazs0,saxs1,says1,sazs1,saxs2,says2,sazs2,paxs0,
     &                 pays0,pazs0,paxs1,pays1,pazs1,daxs0,days0,dazs0,
     &                 daxscond,dayscond,dazscond)

c  this subroutine written by roger gallant and alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

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

      dimension icontr(*),inddyy(m,*)

      dimension g(*),density(*),zeta(*),px(*),py(*),pz(*),pxax(*),
     &          pyay(*),pzaz(*),base(*),ss0(*),ss1(*),ss2(*),ps0(n,*),
     &          ps1(n,*),ds0(n,*),dscond(m,*),c(*),pxcx(*),pycy(*),
     &          pzcz(*),u(*),twozeta(*),ss3(*),saxs0(*),says0(*),
     &          sazs0(*),saxs1(*),says1(*),sazs1(*),saxs2(*),says2(*),
     &          sazs2(*),paxs0(n,*),pays0(n,*),pazs0(n,*),paxs1(n,*),
     &          pays1(n,*),pazs1(n,*),cx(*),cy(*),cz(*),daxs0(n,*),
     &          days0(n,*),dazs0(n,*),daxscond(m,*),dayscond(m,*),
     &          dazscond(m,*)

      data zero,half,one,two,pi/0.0d0,0.5d0,
     &                          1.0d0,2.0d0,3.1415926535898d0/

      do 1001 i=1,n
 1001 twozeta(i)=two*zeta(i)

      twopi=two*pi

      do 1002 i=1,n
      pxcx(i)=px(i)-ccx
      pycy(i)=py(i)-ccy
      pzcz(i)=pz(i)-ccz
 1002 u(i)=zeta(i)*(pxcx(i)*pxcx(i)+pycy(i)*pycy(i)+pzcz(i)*pzcz(i))

      do 1003 i=1,n
      call augg(4,u(i),g)
      ss0(i)=base(i)*g(1)*twopi/zeta(i)
      ss1(i)=base(i)*g(2)*twopi/zeta(i)
      ss2(i)=base(i)*g(3)*twopi/zeta(i)
 1003 ss3(i)=base(i)*g(4)*twopi/zeta(i)

      do 1004 i=1,n

      ps0(i,1)=pxax(i)*ss0(i)-pxcx(i)*ss1(i)
      ps0(i,2)=pyay(i)*ss0(i)-pycy(i)*ss1(i)
      ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)

      ps1(i,1)=pxax(i)*ss1(i)-pxcx(i)*ss2(i)
      ps1(i,2)=pyay(i)*ss1(i)-pycy(i)*ss2(i)
      ps1(i,3)=pzaz(i)*ss1(i)-pzcz(i)*ss2(i)

 1004 continue

      do 1005 i=1,n
 1005 c(i)=(ss0(i)-ss1(i))/twozeta(i)

      do 1006 i=1,n
      ds0(i,1)=pxax(i)*ps0(i,1)-pxcx(i)*ps1(i,1)+c(i)
      ds0(i,2)=pyay(i)*ps0(i,1)-pycy(i)*ps1(i,1)
      ds0(i,3)=pzaz(i)*ps0(i,1)-pzcz(i)*ps1(i,1)
      ds0(i,4)=pyay(i)*ps0(i,2)-pycy(i)*ps1(i,2)+c(i)
      ds0(i,5)=pzaz(i)*ps0(i,2)-pzcz(i)*ps1(i,2)
 1006 ds0(i,6)=pzaz(i)*ps0(i,3)-pzcz(i)*ps1(i,3)+c(i)

      do 1007 i=1,n

      saxs0(i)=twozeta(i)*pxcx(i)*ss1(i)
      says0(i)=twozeta(i)*pycy(i)*ss1(i)
      sazs0(i)=twozeta(i)*pzcz(i)*ss1(i)

      saxs1(i)=twozeta(i)*pxcx(i)*ss2(i)
      says1(i)=twozeta(i)*pycy(i)*ss2(i)
      sazs1(i)=twozeta(i)*pzcz(i)*ss2(i)

      saxs2(i)=twozeta(i)*pxcx(i)*ss3(i)
      says2(i)=twozeta(i)*pycy(i)*ss3(i)
      sazs2(i)=twozeta(i)*pzcz(i)*ss3(i)

 1007 continue

      do 1008 i=1,n

      paxs0(i,1)=pxax(i)*saxs0(i)-pxcx(i)*saxs1(i)+ss1(i)
      paxs0(i,2)=pyay(i)*saxs0(i)-pycy(i)*saxs1(i)
      paxs0(i,3)=pzaz(i)*saxs0(i)-pzcz(i)*saxs1(i)
      pays0(i,1)=pxax(i)*says0(i)-pxcx(i)*says1(i)
      pays0(i,2)=pyay(i)*says0(i)-pycy(i)*says1(i)+ss1(i)
      pays0(i,3)=pzaz(i)*says0(i)-pzcz(i)*says1(i)
      pazs0(i,1)=pxax(i)*sazs0(i)-pxcx(i)*sazs1(i)
      pazs0(i,2)=pyay(i)*sazs0(i)-pycy(i)*sazs1(i)
      pazs0(i,3)=pzaz(i)*sazs0(i)-pzcz(i)*sazs1(i)+ss1(i)

      paxs1(i,1)=pxax(i)*saxs1(i)-pxcx(i)*saxs2(i)+ss2(i)
      paxs1(i,2)=pyay(i)*saxs1(i)-pycy(i)*saxs2(i)
      paxs1(i,3)=pzaz(i)*saxs1(i)-pzcz(i)*saxs2(i)
      pays1(i,1)=pxax(i)*says1(i)-pxcx(i)*says2(i)
      pays1(i,2)=pyay(i)*says1(i)-pycy(i)*says2(i)+ss2(i)
      pays1(i,3)=pzaz(i)*says1(i)-pzcz(i)*says2(i)
      pazs1(i,1)=pxax(i)*sazs1(i)-pxcx(i)*sazs2(i)
      pazs1(i,2)=pyay(i)*sazs1(i)-pycy(i)*sazs2(i)
      pazs1(i,3)=pzaz(i)*sazs1(i)-pzcz(i)*sazs2(i)+ss2(i)

 1008 continue

      do 1009 i=1,n
      cx(i)=(saxs0(i)-saxs1(i))/twozeta(i)
      cy(i)=(says0(i)-says1(i))/twozeta(i)
 1009 cz(i)=(sazs0(i)-sazs1(i))/twozeta(i)

      do 1010 i=1,n
      daxs0(i,1)=pxax(i)*paxs0(i,1)-pxcx(i)*paxs1(i,1)+cx(i)+ps1(i,1)
      daxs0(i,2)=pyay(i)*paxs0(i,1)-pycy(i)*paxs1(i,1)
      daxs0(i,3)=pzaz(i)*paxs0(i,1)-pzcz(i)*paxs1(i,1)
      daxs0(i,4)=pyay(i)*paxs0(i,2)-pycy(i)*paxs1(i,2)+cx(i)
      daxs0(i,5)=pzaz(i)*paxs0(i,2)-pzcz(i)*paxs1(i,2)
      daxs0(i,6)=pzaz(i)*paxs0(i,3)-pzcz(i)*paxs1(i,3)+cx(i)
      days0(i,1)=pxax(i)*pays0(i,1)-pxcx(i)*pays1(i,1)+cy(i)
      days0(i,2)=pyay(i)*pays0(i,1)-pycy(i)*pays1(i,1)      +ps1(i,1)
      days0(i,3)=pzaz(i)*pays0(i,1)-pzcz(i)*pays1(i,1)
      days0(i,4)=pyay(i)*pays0(i,2)-pycy(i)*pays1(i,2)+cy(i)+ps1(i,2)
      days0(i,5)=pzaz(i)*pays0(i,2)-pzcz(i)*pays1(i,2)
      days0(i,6)=pzaz(i)*pays0(i,3)-pzcz(i)*pays1(i,3)+cy(i)
      dazs0(i,1)=pxax(i)*pazs0(i,1)-pxcx(i)*pazs1(i,1)+cz(i)
      dazs0(i,2)=pyay(i)*pazs0(i,1)-pycy(i)*pazs1(i,1)
      dazs0(i,3)=pzaz(i)*pazs0(i,1)-pzcz(i)*pazs1(i,1)      +ps1(i,1)
      dazs0(i,4)=pyay(i)*pazs0(i,2)-pycy(i)*pazs1(i,2)+cz(i)
      dazs0(i,5)=pzaz(i)*pazs0(i,2)-pzcz(i)*pazs1(i,2)      +ps1(i,2)
 1010 dazs0(i,6)=pzaz(i)*pazs0(i,3)-pzcz(i)*pazs1(i,3)+cz(i)+ps1(i,3)

      do 1011 j=1,6
      do 1011 i=1,m
 1011 dscond(i,j)=zero

      do 1012 j=1,6
      do 1012 i=1,n
 1012 dscond(icontr(i),j)=dscond(icontr(i),j)+ds0(i,j)

      do 1013 j=1,6
      do 1013 i=1,m
      daxscond(i,j)=zero
      dayscond(i,j)=zero
 1013 dazscond(i,j)=zero

      do 1014 j=1,6
      do 1014 i=1,n
      daxscond(icontr(i),j)=daxscond(icontr(i),j)+daxs0(i,j)
      dayscond(icontr(i),j)=dayscond(icontr(i),j)+days0(i,j)
 1014 dazscond(icontr(i),j)=dazscond(icontr(i),j)+dazs0(i,j)

      do 1015 j=1,6
      do 1015 i=1,m
 1015 esp=esp-density(inddyy(i,j))*dscond(i,j)

      do 1016 j=1,6
      do 1016 i=1,m
      efix=efix-density(inddyy(i,j))*daxscond(i,j)
      efiy=efiy-density(inddyy(i,j))*dayscond(i,j)
 1016 efiz=efiz-density(inddyy(i,j))*dazscond(i,j)

      return
      end
