       subroutine efips(m,n,icontr,inddyy,esp,efix,efiy,efiz,ccx,ccy,
     &                  ccz,g,density,zeta,px,py,pz,pxax,pyay,pzaz,
     &                  base,ss0,ss1,ps0,pscond,pxcx,pycy,pzcz,u,ss2,
     &                  saxs0,says0,sazs0,saxs1,says1,sazs1,paxs0,pays0,
     &                  pazs0,paxscond,payscond,pazscond)

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(*),ps0(n,*),
     &          pscond(m,*),pxcx(*),pycy(*),pzcz(*),u(*),ss2(*),
     &          saxs0(*),says0(*),sazs0(*),saxs1(*),says1(*),sazs1(*),
     &          paxs0(n,*),pays0(n,*),pazs0(n,*),paxscond(m,*),
     &          payscond(m,*),pazscond(m,*)

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

      twopi=two*pi

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

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

      do 1003 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)
 1003 ps0(i,3)=pzaz(i)*ss0(i)-pzcz(i)*ss1(i)
 
      do 1004 i=1,n
      saxs0(i)=two*zeta(i)*pxcx(i)*ss1(i)
      says0(i)=two*zeta(i)*pycy(i)*ss1(i)
 1004 sazs0(i)=two*zeta(i)*pzcz(i)*ss1(i)

      do 1005 i=1,n
      saxs1(i)=two*zeta(i)*pxcx(i)*ss2(i)
      says1(i)=two*zeta(i)*pycy(i)*ss2(i)
 1005 sazs1(i)=two*zeta(i)*pzcz(i)*ss2(i)

      do 1006 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)
 1006 pazs0(i,3)=pzaz(i)*sazs0(i)-pzcz(i)*sazs1(i)+ss1(i)

      do 1007 j=1,3
      do 1007 i=1,m
 1007 pscond(i,j)=zero

      do 1008 j=1,3
      do 1008 i=1,n
 1008 pscond(icontr(i),j)=pscond(icontr(i),j)+ps0(i,j)

      do 1009 j=1,3
      do 1009 i=1,m
      paxscond(i,j)=zero
      payscond(i,j)=zero
 1009 pazscond(i,j)=zero

      do 1010 j=1,3
      do 1010 i=1,n
      paxscond(icontr(i),j)=paxscond(icontr(i),j)+paxs0(i,j)
      payscond(icontr(i),j)=payscond(icontr(i),j)+pays0(i,j)
 1010 pazscond(icontr(i),j)=pazscond(icontr(i),j)+pazs0(i,j)

      do 1011 j=1,3
      do 1011 i=1,m
 1011 esp=esp-density(inddyy(i,j))*pscond(i,j)

      do 1012 j=1,3
      do 1012 i=1,m
      efix=efix-density(inddyy(i,j))*paxscond(i,j)
      efiy=efiy-density(inddyy(i,j))*payscond(i,j)
 1012 efiz=efiz-density(inddyy(i,j))*pazscond(i,j)

      return 
      end 
