      subroutine cdgrxx(ncds,ncdspd,ilfunc,icfunc,cdfitc,alpha,coeffs,
     &                  coeffp,coeffd,coord,gradient)

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)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension g(6)
 
      dimension ilfunc(*),icfunc(*)

      dimension cdfitc(*),alpha(*),coeffs(*),coeffp(*),coeffd(*)
      dimension coord(3,*),gradient(3,*) 

      data one,two,three/1.0d0,2.0d0,3.0d0/
      data four,five,forty9/4.0d0,5.0d0,49.0d0/
      data zero,half,pi/0.0d0,0.5d0,3.1415926535898d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      sqrt493=sqrt(forty9/three)
      factor=(pi**(five/four)*sqrt(two))**2

      ncombos=(ncdspd*(ncdspd+1))/2

      nloop=ncombos/numprocs
      if((nloop*numprocs).lt.ncombos) nloop=nloop+1

      istart=myid*nloop+1
      iend=(myid+1)*nloop

      ij=0

      do 1001 i=1,ncdspd
      do 1001 j=1,i

      ij=ij+1

      if((ij.lt.istart).or.(ij.gt.iend)) goto 1001

      ilabel=ilfunc(ncds+i)
      jlabel=ilfunc(ncds+j)

      nc1=icfunc(ilabel)
      nc2=icfunc(jlabel)

      zetaa=alpha(ilabel)
      zetab=alpha(jlabel)

      axbx=coord(1,nc1)-coord(1,nc2)
      ayby=coord(2,nc1)-coord(2,nc2)
      azbz=coord(3,nc1)-coord(3,nc2)

      const=factor/(zetaa*zetab)

      zetaab=one/(zetaa+zetab)

      wxpx=(zetaa*coord(1,nc1)+zetab*coord(1,nc2))*zetaab-coord(1,nc1)
      wypy=(zetaa*coord(2,nc1)+zetab*coord(2,nc2))*zetaab-coord(2,nc1)
      wzpz=(zetaa*coord(3,nc1)+zetab*coord(3,nc2))*zetaab-coord(3,nc1)

      wxqx=(zetaa*coord(1,nc1)+zetab*coord(1,nc2))*zetaab-coord(1,nc2)
      wyqy=(zetaa*coord(2,nc1)+zetab*coord(2,nc2))*zetaab-coord(2,nc2)
      wzqz=(zetaa*coord(3,nc1)+zetab*coord(3,nc2))*zetaab-coord(3,nc2)

      u=zetaa*zetab*zetaab*(axbx*axbx+ayby*ayby+azbz*azbz)

      call augg(6,u,g)

      ss0=g(1)*const*sqrt(zetaab)
      ss1=g(2)*const*sqrt(zetaab)
      ss2=g(3)*const*sqrt(zetaab)
      ss3=g(4)*const*sqrt(zetaab)
      ss4=g(5)*const*sqrt(zetaab)
      ss5=g(6)*const*sqrt(zetaab)

      spx0=wxqx*ss1
      spy0=wyqy*ss1
      spz0=wzqz*ss1

      spx1=wxqx*ss2
      spy1=wyqy*ss2
      spz1=wzqz*ss2

      cfact=(ss0-zetaa*zetaab*ss1)/(two*zetab)

      sdxx0=wxqx*spx1+cfact
      sdxy0=wyqy*spx1
      sdxz0=wzqz*spx1
      sdyy0=wyqy*spy1+cfact
      sdyz0=wzqz*spy1
      sdzz0=wzqz*spz1+cfact

      pxs0=wxpx*ss1
      pys0=wypy*ss1
      pzs0=wzpz*ss1

      pxs1=wxpx*ss2
      pys1=wypy*ss2
      pzs1=wzpz*ss2

      pxs2=wxpx*ss3
      pys2=wypy*ss3
      pzs2=wzpz*ss3

      pxs3=wxpx*ss4
      pys3=wypy*ss4
      pzs3=wzpz*ss4

      pxs4=wxpx*ss5
      pys4=wypy*ss5
      pzs4=wzpz*ss5

      cfact=ss1*zetaab/two

      pxpx0=wxqx*pxs1+cfact
      pxpy0=wyqy*pxs1
      pxpz0=wzqz*pxs1
      pypx0=wxqx*pys1
      pypy0=wyqy*pys1+cfact
      pypz0=wzqz*pys1
      pzpx0=wxqx*pzs1
      pzpy0=wyqy*pzs1
      pzpz0=wzqz*pzs1+cfact

      cfact=ss2*zetaab/two

      pxpx1=wxqx*pxs2+cfact
      pxpy1=wyqy*pxs2
      pxpz1=wzqz*pxs2
      pypx1=wxqx*pys2
      pypy1=wyqy*pys2+cfact
      pypz1=wzqz*pys2
      pzpx1=wxqx*pzs2
      pzpy1=wyqy*pzs2
      pzpz1=wzqz*pzs2+cfact

      cfactx=(pxs0-zetaa*zetaab*pxs1)/(two*zetab)
      cfacty=(pys0-zetaa*zetaab*pys1)/(two*zetab)
      cfactz=(pzs0-zetaa*zetaab*pzs1)/(two*zetab)

      cfactpx=spx1*zetaab/two
      cfactpy=spy1*zetaab/two
      cfactpz=spz1*zetaab/two

      pxdxx0=wxqx*pxpx1+cfactx+cfactpx
      pxdxy0=wyqy*pxpx1
      pxdxz0=wzqz*pxpx1
      pxdyy0=wyqy*pxpy1+cfactx
      pxdyz0=wzqz*pxpy1
      pxdzz0=wzqz*pxpz1+cfactx
      pydxx0=wxqx*pypx1+cfacty
      pydxy0=wyqy*pypx1       +cfactpx
      pydxz0=wzqz*pypx1
      pydyy0=wyqy*pypy1+cfacty+cfactpy
      pydyz0=wzqz*pypy1
      pydzz0=wzqz*pypz1+cfacty
      pzdxx0=wxqx*pzpx1+cfactz
      pzdxy0=wyqy*pzpx1
      pzdxz0=wzqz*pzpx1       +cfactpx
      pzdyy0=wyqy*pzpy1+cfactz
      pzdyz0=wzqz*pzpy1       +cfactpy
      pzdzz0=wzqz*pzpz1+cfactz+cfactpz

      cfact=(ss0-zetab*zetaab*ss1)/(two*zetaa)

      dxxs0=wxpx*pxs1+cfact
      dxys0=wypy*pxs1
      dxzs0=wzpz*pxs1
      dyys0=wypy*pys1+cfact
      dyzs0=wzpz*pys1
      dzzs0=wzpz*pzs1+cfact

      cfact=(ss1-zetab*zetaab*ss2)/(two*zetaa)

      dxxs1=wxpx*pxs2+cfact
      dxys1=wypy*pxs2
      dxzs1=wzpz*pxs2
      dyys1=wypy*pys2+cfact
      dyzs1=wzpz*pys2
      dzzs1=wzpz*pzs2+cfact

      cfact=(ss2-zetab*zetaab*ss3)/(two*zetaa)

      dxxs2=wxpx*pxs3+cfact
      dxys2=wypy*pxs3
      dxzs2=wzpz*pxs3
      dyys2=wypy*pys3+cfact
      dyzs2=wzpz*pys3
      dzzs2=wzpz*pzs3+cfact

      cfact=(ss3-zetab*zetaab*ss4)/(two*zetaa)

      dxxs3=wxpx*pxs4+cfact
      dxys3=wypy*pxs4
      dxzs3=wzpz*pxs4
      dyys3=wypy*pys4+cfact
      dyzs3=wzpz*pys4
      dzzs3=wzpz*pzs4+cfact

      cfactx=pxs1*zetaab/two
      cfacty=pys1*zetaab/two
      cfactz=pzs1*zetaab/two

      dxxpx0=wxqx*dxxs1+cfactx+cfactx
      dxxpy0=wyqy*dxxs1
      dxxpz0=wzqz*dxxs1
      dxypx0=wxqx*dxys1+cfacty
      dxypy0=wyqy*dxys1+cfactx
      dxypz0=wzqz*dxys1
      dxzpx0=wxqx*dxzs1+cfactz
      dxzpy0=wyqy*dxzs1
      dxzpz0=wzqz*dxzs1+cfactx
      dyypx0=wxqx*dyys1
      dyypy0=wyqy*dyys1+cfacty+cfacty
      dyypz0=wzqz*dyys1
      dyzpx0=wxqx*dyzs1
      dyzpy0=wyqy*dyzs1+cfactz
      dyzpz0=wzqz*dyzs1+cfacty
      dzzpx0=wxqx*dzzs1
      dzzpy0=wyqy*dzzs1
      dzzpz0=wzqz*dzzs1+cfactz+cfactz

      cfactx=pxs2*zetaab/two
      cfacty=pys2*zetaab/two
      cfactz=pzs2*zetaab/two

      dxxpx1=wxqx*dxxs2+cfactx+cfactx
      dxxpy1=wyqy*dxxs2
      dxxpz1=wzqz*dxxs2
      dxypx1=wxqx*dxys2+cfacty
      dxypy1=wyqy*dxys2+cfactx
      dxypz1=wzqz*dxys2
      dxzpx1=wxqx*dxzs2+cfactz
      dxzpy1=wyqy*dxzs2
      dxzpz1=wzqz*dxzs2+cfactx
      dyypx1=wxqx*dyys2
      dyypy1=wyqy*dyys2+cfacty+cfacty
      dyypz1=wzqz*dyys2
      dyzpx1=wxqx*dyzs2
      dyzpy1=wyqy*dyzs2+cfactz
      dyzpz1=wzqz*dyzs2+cfacty
      dzzpx1=wxqx*dzzs2
      dzzpy1=wyqy*dzzs2
      dzzpz1=wzqz*dzzs2+cfactz+cfactz

      cfactxx=(dxxs0-zetaa*zetaab*dxxs1)/(two*zetab)
      cfactxy=(dxys0-zetaa*zetaab*dxys1)/(two*zetab)
      cfactxz=(dxzs0-zetaa*zetaab*dxzs1)/(two*zetab)
      cfactyy=(dyys0-zetaa*zetaab*dyys1)/(two*zetab)
      cfactyz=(dyzs0-zetaa*zetaab*dyzs1)/(two*zetab)
      cfactzz=(dzzs0-zetaa*zetaab*dzzs1)/(two*zetab)

      cfactxpx=pxpx1*zetaab/two
      cfactxpy=pxpy1*zetaab/two
      cfactxpz=pxpz1*zetaab/two
      cfactypx=pypx1*zetaab/two
      cfactypy=pypy1*zetaab/two
      cfactypz=pypz1*zetaab/two
      cfactzpx=pzpx1*zetaab/two
      cfactzpy=pzpy1*zetaab/two
      cfactzpz=pzpz1*zetaab/two

      dxxdxx0=wxqx*dxxpx1+cfactxx+cfactxpx+cfactxpx
      dxxdxy0=wyqy*dxxpx1
      dxxdxz0=wzqz*dxxpx1
      dxxdyy0=wyqy*dxxpy1+cfactxx
      dxxdyz0=wzqz*dxxpy1
      dxxdzz0=wzqz*dxxpz1+cfactxx
      dxydxx0=wxqx*dxypx1+cfactxy+cfactypx
      dxydxy0=wyqy*dxypx1        +cfactxpx
      dxydxz0=wzqz*dxypx1
      dxydyy0=wyqy*dxypy1+cfactxy+cfactxpy
      dxydyz0=wzqz*dxypy1
      dxydzz0=wzqz*dxypz1+cfactxy
      dxzdxx0=wxqx*dxzpx1+cfactxz+cfactzpx
      dxzdxy0=wyqy*dxzpx1
      dxzdxz0=wzqz*dxzpx1        +cfactxpx
      dxzdyy0=wyqy*dxzpy1+cfactxz
      dxzdyz0=wzqz*dxzpy1        +cfactxpy
      dxzdzz0=wzqz*dxzpz1+cfactxz+cfactxpz
      dyydxx0=wxqx*dyypx1+cfactyy
      dyydxy0=wyqy*dyypx1        +cfactypx+cfactypx
      dyydxz0=wzqz*dyypx1
      dyydyy0=wyqy*dyypy1+cfactyy+cfactypy+cfactypy
      dyydyz0=wzqz*dyypy1
      dyydzz0=wzqz*dyypz1+cfactyy
      dyzdxx0=wxqx*dyzpx1+cfactyz
      dyzdxy0=wyqy*dyzpx1        +cfactzpx
      dyzdxz0=wzqz*dyzpx1        +cfactypx
      dyzdyy0=wyqy*dyzpy1+cfactyz+cfactzpy
      dyzdyz0=wzqz*dyzpy1        +cfactypy
      dyzdzz0=wzqz*dyzpz1+cfactyz+cfactypz
      dzzdxx0=wxqx*dzzpx1+cfactzz
      dzzdxy0=wyqy*dzzpx1
      dzzdxz0=wzqz*dzzpx1        +cfactzpx+cfactzpx
      dzzdyy0=wyqy*dzzpy1+cfactzz
      dzzdyz0=wzqz*dzzpy1        +cfactzpy+cfactzpy
      dzzdzz0=wzqz*dzzpz1+cfactzz+cfactzpz+cfactzpz

      cfactx=(pxs0-zetab*zetaab*pxs1)/(two*zetaa)
      cfacty=(pys0-zetab*zetaab*pys1)/(two*zetaa)
      cfactz=(pzs0-zetab*zetaab*pzs1)/(two*zetaa)

      fxxxs0=wxpx*dxxs1+cfactx+cfactx
      fxxys0=wypy*dxxs1
      fxxzs0=wzpz*dxxs1
      fxyys0=wypy*dxys1+cfactx
      fxyzs0=wzpz*dxys1
      fxzzs0=wzpz*dxzs1+cfactx
      fyyys0=wypy*dyys1+cfacty+cfacty
      fyyzs0=wzpz*dyys1
      fyzzs0=wzpz*dyzs1+cfacty
      fzzzs0=wzpz*dzzs1+cfactz+cfactz

      cfactx=(pxs1-zetab*zetaab*pxs2)/(two*zetaa)
      cfacty=(pys1-zetab*zetaab*pys2)/(two*zetaa)
      cfactz=(pzs1-zetab*zetaab*pzs2)/(two*zetaa)

      fxxxs1=wxpx*dxxs2+cfactx+cfactx
      fxxys1=wypy*dxxs2
      fxxzs1=wzpz*dxxs2
      fxyys1=wypy*dxys2+cfactx
      fxyzs1=wzpz*dxys2
      fxzzs1=wzpz*dxzs2+cfactx
      fyyys1=wypy*dyys2+cfacty+cfacty
      fyyzs1=wzpz*dyys2
      fyzzs1=wzpz*dyzs2+cfacty
      fzzzs1=wzpz*dzzs2+cfactz+cfactz

      cfactx=(pxs2-zetab*zetaab*pxs3)/(two*zetaa)
      cfacty=(pys2-zetab*zetaab*pys3)/(two*zetaa)
      cfactz=(pzs2-zetab*zetaab*pzs3)/(two*zetaa)

      fxxxs2=wxpx*dxxs3+cfactx+cfactx
      fxxys2=wypy*dxxs3
      fxxzs2=wzpz*dxxs3
      fxyys2=wypy*dxys3+cfactx
      fxyzs2=wzpz*dxys3
      fxzzs2=wzpz*dxzs3+cfactx
      fyyys2=wypy*dyys3+cfacty+cfacty
      fyyzs2=wzpz*dyys3
      fyzzs2=wzpz*dyzs3+cfacty
      fzzzs2=wzpz*dzzs3+cfactz+cfactz

      cfactxx=dxxs1*zetaab/two
      cfactxy=dxys1*zetaab/two
      cfactxz=dxzs1*zetaab/two
      cfactyy=dyys1*zetaab/two
      cfactyz=dyzs1*zetaab/two
      cfactzz=dzzs1*zetaab/two

      fxxxpx0=wxqx*fxxxs1+cfactxx+cfactxx+cfactxx
      fxxxpy0=wyqy*fxxxs1
      fxxxpz0=wzqz*fxxxs1
      fxxypx0=wxqx*fxxys1+cfactxy+cfactxy
      fxxypy0=wyqy*fxxys1+cfactxx
      fxxypz0=wzqz*fxxys1
      fxxzpx0=wxqx*fxxzs1+cfactxz+cfactxz
      fxxzpy0=wyqy*fxxzs1
      fxxzpz0=wzqz*fxxzs1+cfactxx
      fxyypx0=wxqx*fxyys1+cfactyy
      fxyypy0=wyqy*fxyys1+cfactxy+cfactxy
      fxyypz0=wzqz*fxyys1
      fxyzpx0=wxqx*fxyzs1+cfactyz
      fxyzpy0=wyqy*fxyzs1+cfactxz
      fxyzpz0=wzqz*fxyzs1+cfactxy
      fxzzpx0=wxqx*fxzzs1+cfactzz
      fxzzpy0=wyqy*fxzzs1
      fxzzpz0=wzqz*fxzzs1+cfactxz+cfactxz
      fyyypx0=wxqx*fyyys1
      fyyypy0=wyqy*fyyys1+cfactyy+cfactyy+cfactyy
      fyyypz0=wzqz*fyyys1
      fyyzpx0=wxqx*fyyzs1
      fyyzpy0=wyqy*fyyzs1+cfactyz+cfactyz
      fyyzpz0=wzqz*fyyzs1+cfactyy
      fyzzpx0=wxqx*fyzzs1
      fyzzpy0=wyqy*fyzzs1+cfactzz
      fyzzpz0=wzqz*fyzzs1+cfactyz+cfactyz
      fzzzpx0=wxqx*fzzzs1
      fzzzpy0=wyqy*fzzzs1
      fzzzpz0=wzqz*fzzzs1+cfactzz+cfactzz+cfactzz

      cfactxx=dxxs2*zetaab/two
      cfactxy=dxys2*zetaab/two
      cfactxz=dxzs2*zetaab/two
      cfactyy=dyys2*zetaab/two
      cfactyz=dyzs2*zetaab/two
      cfactzz=dzzs2*zetaab/two

      fxxxpx1=wxqx*fxxxs2+cfactxx+cfactxx+cfactxx
      fxxxpy1=wyqy*fxxxs2
      fxxxpz1=wzqz*fxxxs2
      fxxypx1=wxqx*fxxys2+cfactxy+cfactxy
      fxxypy1=wyqy*fxxys2+cfactxx
      fxxypz1=wzqz*fxxys2
      fxxzpx1=wxqx*fxxzs2+cfactxz+cfactxz
      fxxzpy1=wyqy*fxxzs2
      fxxzpz1=wzqz*fxxzs2+cfactxx
      fxyypx1=wxqx*fxyys2+cfactyy
      fxyypy1=wyqy*fxyys2+cfactxy+cfactxy
      fxyypz1=wzqz*fxyys2
      fxyzpx1=wxqx*fxyzs2+cfactyz
      fxyzpy1=wyqy*fxyzs2+cfactxz
      fxyzpz1=wzqz*fxyzs2+cfactxy
      fxzzpx1=wxqx*fxzzs2+cfactzz
      fxzzpy1=wyqy*fxzzs2
      fxzzpz1=wzqz*fxzzs2+cfactxz+cfactxz
      fyyypx1=wxqx*fyyys2
      fyyypy1=wyqy*fyyys2+cfactyy+cfactyy+cfactyy
      fyyypz1=wzqz*fyyys2
      fyyzpx1=wxqx*fyyzs2
      fyyzpy1=wyqy*fyyzs2+cfactyz+cfactyz
      fyyzpz1=wzqz*fyyzs2+cfactyy
      fyzzpx1=wxqx*fyzzs2
      fyzzpy1=wyqy*fyzzs2+cfactzz
      fyzzpz1=wzqz*fyzzs2+cfactyz+cfactyz
      fzzzpx1=wxqx*fzzzs2
      fzzzpy1=wyqy*fzzzs2
      fzzzpz1=wzqz*fzzzs2+cfactzz+cfactzz+cfactzz

      cfactxxx=(fxxxs0-zetaa*zetaab*fxxxs1)/(two*zetab)
      cfactxxy=(fxxys0-zetaa*zetaab*fxxys1)/(two*zetab)
      cfactxxz=(fxxzs0-zetaa*zetaab*fxxzs1)/(two*zetab)
      cfactxyy=(fxyys0-zetaa*zetaab*fxyys1)/(two*zetab)
      cfactxyz=(fxyzs0-zetaa*zetaab*fxyzs1)/(two*zetab)
      cfactxzz=(fxzzs0-zetaa*zetaab*fxzzs1)/(two*zetab)
      cfactyyy=(fyyys0-zetaa*zetaab*fyyys1)/(two*zetab)
      cfactyyz=(fyyzs0-zetaa*zetaab*fyyzs1)/(two*zetab)
      cfactyzz=(fyzzs0-zetaa*zetaab*fyzzs1)/(two*zetab)
      cfactzzz=(fzzzs0-zetaa*zetaab*fzzzs1)/(two*zetab)

      cfactxxpx=dxxpx1*zetaab/two
      cfactxxpy=dxxpy1*zetaab/two
      cfactxxpz=dxxpz1*zetaab/two
      cfactxypx=dxypx1*zetaab/two
      cfactxypy=dxypy1*zetaab/two
      cfactxypz=dxypz1*zetaab/two
      cfactxzpx=dxzpx1*zetaab/two
      cfactxzpy=dxzpy1*zetaab/two
      cfactxzpz=dxzpz1*zetaab/two
      cfactyypx=dyypx1*zetaab/two
      cfactyypy=dyypy1*zetaab/two
      cfactyypz=dyypz1*zetaab/two
      cfactyzpx=dyzpx1*zetaab/two
      cfactyzpy=dyzpy1*zetaab/two
      cfactyzpz=dyzpz1*zetaab/two
      cfactzzpx=dzzpx1*zetaab/two
      cfactzzpy=dzzpy1*zetaab/two
      cfactzzpz=dzzpz1*zetaab/two

      fxxxdxx0=wxqx*fxxxpx1+cfactxxx+cfactxxpx+cfactxxpx+cfactxxpx
      fxxxdxy0=wyqy*fxxxpx1
      fxxxdxz0=wzqz*fxxxpx1
      fxxxdyy0=wyqy*fxxxpy1+cfactxxx
      fxxxdyz0=wzqz*fxxxpy1
      fxxxdzz0=wzqz*fxxxpz1+cfactxxx
      fxxydxx0=wxqx*fxxypx1+cfactxxy+cfactxypx+cfactxypx
      fxxydxy0=wyqy*fxxypx1         +cfactxxpx
      fxxydxz0=wzqz*fxxypx1
      fxxydyy0=wyqy*fxxypy1+cfactxxy+cfactxxpy
      fxxydyz0=wzqz*fxxypy1
      fxxydzz0=wzqz*fxxypz1+cfactxxy
      fxxzdxx0=wxqx*fxxzpx1+cfactxxz+cfactxzpx+cfactxzpx
      fxxzdxy0=wyqy*fxxzpx1
      fxxzdxz0=wzqz*fxxzpx1         +cfactxxpx
      fxxzdyy0=wyqy*fxxzpy1+cfactxxz
      fxxzdyz0=wzqz*fxxzpy1         +cfactxxpy
      fxxzdzz0=wzqz*fxxzpz1+cfactxxz+cfactxxpz
      fxyydxx0=wxqx*fxyypx1+cfactxyy+cfactyypx
      fxyydxy0=wyqy*fxyypx1         +cfactxypx+cfactxypx
      fxyydxz0=wzqz*fxyypx1
      fxyydyy0=wyqy*fxyypy1+cfactxyy+cfactxypy+cfactxypy
      fxyydyz0=wzqz*fxyypy1
      fxyydzz0=wzqz*fxyypz1+cfactxyy
      fxyzdxx0=wxqx*fxyzpx1+cfactxyz+cfactyzpx
      fxyzdxy0=wyqy*fxyzpx1         +cfactxzpx
      fxyzdxz0=wzqz*fxyzpx1         +cfactxypx
      fxyzdyy0=wyqy*fxyzpy1+cfactxyz+cfactxzpy
      fxyzdyz0=wzqz*fxyzpy1         +cfactxypy
      fxyzdzz0=wzqz*fxyzpz1+cfactxyz+cfactxypz
      fxzzdxx0=wxqx*fxzzpx1+cfactxzz+cfactzzpx
      fxzzdxy0=wyqy*fxzzpx1
      fxzzdxz0=wzqz*fxzzpx1         +cfactxzpx+cfactxzpx
      fxzzdyy0=wyqy*fxzzpy1+cfactxzz
      fxzzdyz0=wzqz*fxzzpy1         +cfactxzpy+cfactxzpy
      fxzzdzz0=wzqz*fxzzpz1+cfactxzz+cfactxzpz+cfactxzpz
      fyyydxx0=wxqx*fyyypx1+cfactyyy
      fyyydxy0=wyqy*fyyypx1         +cfactyypx+cfactyypx+cfactyypx
      fyyydxz0=wzqz*fyyypx1
      fyyydyy0=wyqy*fyyypy1+cfactyyy+cfactyypy+cfactyypy+cfactyypy
      fyyydyz0=wzqz*fyyypy1
      fyyydzz0=wzqz*fyyypz1+cfactyyy
      fyyzdxx0=wxqx*fyyzpx1+cfactyyz
      fyyzdxy0=wyqy*fyyzpx1         +cfactyzpx+cfactyzpx
      fyyzdxz0=wzqz*fyyzpx1         +cfactyypx
      fyyzdyy0=wyqy*fyyzpy1+cfactyyz+cfactyzpy+cfactyzpy
      fyyzdyz0=wzqz*fyyzpy1         +cfactyypy
      fyyzdzz0=wzqz*fyyzpz1+cfactyyz+cfactyypz
      fyzzdxx0=wxqx*fyzzpx1+cfactyzz
      fyzzdxy0=wyqy*fyzzpx1         +cfactzzpx
      fyzzdxz0=wzqz*fyzzpx1         +cfactyzpx+cfactyzpx
      fyzzdyy0=wyqy*fyzzpy1+cfactyzz+cfactzzpy
      fyzzdyz0=wzqz*fyzzpy1         +cfactyzpy+cfactyzpy
      fyzzdzz0=wzqz*fyzzpz1+cfactyzz+cfactyzpz+cfactyzpz
      fzzzdxx0=wxqx*fzzzpx1+cfactzzz
      fzzzdxy0=wyqy*fzzzpx1
      fzzzdxz0=wzqz*fzzzpx1         +cfactzzpx+cfactzzpx+cfactzzpx
      fzzzdyy0=wyqy*fzzzpy1+cfactzzz
      fzzzdyz0=wzqz*fzzzpy1         +cfactzzpy+cfactzzpy+cfactzzpy
      fzzzdzz0=wzqz*fzzzpz1+cfactzzz+cfactzzpz+cfactzzpz+cfactzzpz

      if(i.ne.j) then
                       scale=one
                 else
                       scale=half
                 endif

      n=ncds+(i-1)*10

      cfa01=cdfitc(n+01)*scale
      cfa02=cdfitc(n+02)*scale
      cfa03=cdfitc(n+03)*scale
      cfa04=cdfitc(n+04)*scale
      cfa05=cdfitc(n+05)*scale/sqrt493
      cfa06=cdfitc(n+06)*scale
      cfa07=cdfitc(n+07)*scale
      cfa08=cdfitc(n+08)*scale/sqrt493
      cfa09=cdfitc(n+09)*scale
      cfa10=cdfitc(n+10)*scale/sqrt493

      n=ncds+(j-1)*10

      cfb01=cdfitc(n+01)
      cfb02=cdfitc(n+02)
      cfb03=cdfitc(n+03)
      cfb04=cdfitc(n+04)
      cfb05=cdfitc(n+05)/sqrt493
      cfb06=cdfitc(n+06)
      cfb07=cdfitc(n+07)
      cfb08=cdfitc(n+08)/sqrt493
      cfb09=cdfitc(n+09)
      cfb10=cdfitc(n+10)/sqrt493

      css=cfa01*cfb01*coeffs(ilabel)*coeffs(jlabel)

      cpxs=cfa02*cfb01*coeffp(ilabel)*coeffs(jlabel)
      cpys=cfa03*cfb01*coeffp(ilabel)*coeffs(jlabel)
      cpzs=cfa04*cfb01*coeffp(ilabel)*coeffs(jlabel)

      cdxxs=cfa05*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cdxys=cfa06*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cdxzs=cfa07*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cdyys=cfa08*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cdyzs=cfa09*cfb01*coeffd(ilabel)*coeffs(jlabel)
      cdzzs=cfa10*cfb01*coeffd(ilabel)*coeffs(jlabel)

      cspx=cfa01*cfb02*coeffs(ilabel)*coeffp(jlabel)
      cspy=cfa01*cfb03*coeffs(ilabel)*coeffp(jlabel)
      cspz=cfa01*cfb04*coeffs(ilabel)*coeffp(jlabel)

      cpxpx=cfa02*cfb02*coeffp(ilabel)*coeffp(jlabel)
      cpxpy=cfa02*cfb03*coeffp(ilabel)*coeffp(jlabel)
      cpxpz=cfa02*cfb04*coeffp(ilabel)*coeffp(jlabel)
      cpypx=cfa03*cfb02*coeffp(ilabel)*coeffp(jlabel)
      cpypy=cfa03*cfb03*coeffp(ilabel)*coeffp(jlabel)
      cpypz=cfa03*cfb04*coeffp(ilabel)*coeffp(jlabel)
      cpzpx=cfa04*cfb02*coeffp(ilabel)*coeffp(jlabel)
      cpzpy=cfa04*cfb03*coeffp(ilabel)*coeffp(jlabel)
      cpzpz=cfa04*cfb04*coeffp(ilabel)*coeffp(jlabel)

      cdxxpx=cfa05*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdxxpy=cfa05*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdxxpz=cfa05*cfb04*coeffd(ilabel)*coeffp(jlabel)
      cdxypx=cfa06*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdxypy=cfa06*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdxypz=cfa06*cfb04*coeffd(ilabel)*coeffp(jlabel)
      cdxzpx=cfa07*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdxzpy=cfa07*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdxzpz=cfa07*cfb04*coeffd(ilabel)*coeffp(jlabel)
      cdyypx=cfa08*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdyypy=cfa08*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdyypz=cfa08*cfb04*coeffd(ilabel)*coeffp(jlabel)
      cdyzpx=cfa09*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdyzpy=cfa09*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdyzpz=cfa09*cfb04*coeffd(ilabel)*coeffp(jlabel)
      cdzzpx=cfa10*cfb02*coeffd(ilabel)*coeffp(jlabel)
      cdzzpy=cfa10*cfb03*coeffd(ilabel)*coeffp(jlabel)
      cdzzpz=cfa10*cfb04*coeffd(ilabel)*coeffp(jlabel)

      csdxx=cfa01*cfb05*coeffs(ilabel)*coeffd(jlabel)
      csdxy=cfa01*cfb06*coeffs(ilabel)*coeffd(jlabel)
      csdxz=cfa01*cfb07*coeffs(ilabel)*coeffd(jlabel)
      csdyy=cfa01*cfb08*coeffs(ilabel)*coeffd(jlabel)
      csdyz=cfa01*cfb09*coeffs(ilabel)*coeffd(jlabel)
      csdzz=cfa01*cfb10*coeffs(ilabel)*coeffd(jlabel)

      cpxdxx=cfa02*cfb05*coeffp(ilabel)*coeffd(jlabel)
      cpxdxy=cfa02*cfb06*coeffp(ilabel)*coeffd(jlabel)
      cpxdxz=cfa02*cfb07*coeffp(ilabel)*coeffd(jlabel)
      cpxdyy=cfa02*cfb08*coeffp(ilabel)*coeffd(jlabel)
      cpxdyz=cfa02*cfb09*coeffp(ilabel)*coeffd(jlabel)
      cpxdzz=cfa02*cfb10*coeffp(ilabel)*coeffd(jlabel)
      cpydxx=cfa03*cfb05*coeffp(ilabel)*coeffd(jlabel)
      cpydxy=cfa03*cfb06*coeffp(ilabel)*coeffd(jlabel)
      cpydxz=cfa03*cfb07*coeffp(ilabel)*coeffd(jlabel)
      cpydyy=cfa03*cfb08*coeffp(ilabel)*coeffd(jlabel)
      cpydyz=cfa03*cfb09*coeffp(ilabel)*coeffd(jlabel)
      cpydzz=cfa03*cfb10*coeffp(ilabel)*coeffd(jlabel)
      cpzdxx=cfa04*cfb05*coeffp(ilabel)*coeffd(jlabel)
      cpzdxy=cfa04*cfb06*coeffp(ilabel)*coeffd(jlabel)
      cpzdxz=cfa04*cfb07*coeffp(ilabel)*coeffd(jlabel)
      cpzdyy=cfa04*cfb08*coeffp(ilabel)*coeffd(jlabel)
      cpzdyz=cfa04*cfb09*coeffp(ilabel)*coeffd(jlabel)
      cpzdzz=cfa04*cfb10*coeffp(ilabel)*coeffd(jlabel)

      cdxxdxx=cfa05*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdxxdxy=cfa05*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdxxdxz=cfa05*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdxxdyy=cfa05*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdxxdyz=cfa05*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdxxdzz=cfa05*cfb10*coeffd(ilabel)*coeffd(jlabel)
      cdxydxx=cfa06*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdxydxy=cfa06*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdxydxz=cfa06*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdxydyy=cfa06*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdxydyz=cfa06*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdxydzz=cfa06*cfb10*coeffd(ilabel)*coeffd(jlabel)
      cdxzdxx=cfa07*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdxzdxy=cfa07*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdxzdxz=cfa07*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdxzdyy=cfa07*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdxzdyz=cfa07*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdxzdzz=cfa07*cfb10*coeffd(ilabel)*coeffd(jlabel)
      cdyydxx=cfa08*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdyydxy=cfa08*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdyydxz=cfa08*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdyydyy=cfa08*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdyydyz=cfa08*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdyydzz=cfa08*cfb10*coeffd(ilabel)*coeffd(jlabel)
      cdyzdxx=cfa09*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdyzdxy=cfa09*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdyzdxz=cfa09*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdyzdyy=cfa09*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdyzdyz=cfa09*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdyzdzz=cfa09*cfb10*coeffd(ilabel)*coeffd(jlabel)
      cdzzdxx=cfa10*cfb05*coeffd(ilabel)*coeffd(jlabel)
      cdzzdxy=cfa10*cfb06*coeffd(ilabel)*coeffd(jlabel)
      cdzzdxz=cfa10*cfb07*coeffd(ilabel)*coeffd(jlabel)
      cdzzdyy=cfa10*cfb08*coeffd(ilabel)*coeffd(jlabel)
      cdzzdyz=cfa10*cfb09*coeffd(ilabel)*coeffd(jlabel)
      cdzzdzz=cfa10*cfb10*coeffd(ilabel)*coeffd(jlabel)

      factorx=zero
      factory=zero
      factorz=zero

      zz=two*zetaa

      factorx=factorx-zz*css*pxs0
      factory=factory-zz*css*pys0
      factorz=factorz-zz*css*pzs0

      factorx=factorx-zz*(cspx*pxpx0+cspy*pxpy0+cspz*pxpz0)
      factory=factory-zz*(cspx*pypx0+cspy*pypy0+cspz*pypz0)
      factorz=factorz-zz*(cspx*pzpx0+cspy*pzpy0+cspz*pzpz0)

      factorx=factorx-zz*(csdxx*pxdxx0+csdxy*pxdxy0+csdxz*pxdxz0
     &                   +csdyy*pxdyy0+csdyz*pxdyz0+csdzz*pxdzz0)
      factory=factory-zz*(csdxx*pydxx0+csdxy*pydxy0+csdxz*pydxz0
     &                   +csdyy*pydyy0+csdyz*pydyz0+csdzz*pydzz0)
      factorz=factorz-zz*(csdxx*pzdxx0+csdxy*pzdxy0+csdxz*pzdxz0
     &                   +csdyy*pzdyy0+csdyz*pzdyz0+csdzz*pzdzz0)

      factorx=factorx+cpxs*ss0
      factory=factory+cpys*ss0
      factorz=factorz+cpzs*ss0

      factorx=factorx-zz*(cpxs*dxxs0+cpys*dxys0+cpzs*dxzs0)
      factory=factory-zz*(cpxs*dxys0+cpys*dyys0+cpzs*dyzs0)
      factorz=factorz-zz*(cpxs*dxzs0+cpys*dyzs0+cpzs*dzzs0)

      factorx=factorx+cdxxs*pxs0+cdxys*pys0+cdxzs*pzs0+cdxxs*pxs0
      factory=factory+cdxys*pxs0+cdyys*pys0+cdyzs*pzs0+cdyys*pys0
      factorz=factorz+cdxzs*pxs0+cdyzs*pys0+cdzzs*pzs0+cdzzs*pzs0

      factorx=factorx-zz*(cdxxs*fxxxs0+cdxys*fxxys0+cdxzs*fxxzs0
     &                   +cdyys*fxyys0+cdyzs*fxyzs0+cdzzs*fxzzs0)
      factory=factory-zz*(cdxxs*fxxys0+cdxys*fxyys0+cdxzs*fxyzs0
     &                   +cdyys*fyyys0+cdyzs*fyyzs0+cdzzs*fyzzs0)
      factorz=factorz-zz*(cdxxs*fxxzs0+cdxys*fxyzs0+cdxzs*fxzzs0
     &                   +cdyys*fyyzs0+cdyzs*fyzzs0+cdzzs*fzzzs0)

      factorx=factorx+cpxpx*spx0+cpxpy*spy0+cpxpz*spz0
      factory=factory+cpypx*spx0+cpypy*spy0+cpypz*spz0
      factorz=factorz+cpzpx*spx0+cpzpy*spy0+cpzpz*spz0

      factorx=factorx-zz*(cpxpx*dxxpx0+cpxpy*dxxpy0+cpxpz*dxxpz0
     &                   +cpypx*dxypx0+cpypy*dxypy0+cpypz*dxypz0
     &                   +cpzpx*dxzpx0+cpzpy*dxzpy0+cpzpz*dxzpz0)
      factory=factory-zz*(cpxpx*dxypx0+cpxpy*dxypy0+cpxpz*dxypz0
     &                   +cpypx*dyypx0+cpypy*dyypy0+cpypz*dyypz0
     &                   +cpzpx*dyzpx0+cpzpy*dyzpy0+cpzpz*dyzpz0)
      factorz=factorz-zz*(cpxpx*dxzpx0+cpxpy*dxzpy0+cpxpz*dxzpz0
     &                   +cpypx*dyzpx0+cpypy*dyzpy0+cpypz*dyzpz0
     &                   +cpzpx*dzzpx0+cpzpy*dzzpy0+cpzpz*dzzpz0)

      factorx=factorx+cdxxpx*pxpx0+cdxypx*pypx0+cdxzpx*pzpx0
     &               +cdxxpy*pxpy0+cdxypy*pypy0+cdxzpy*pzpy0
     &               +cdxxpz*pxpz0+cdxypz*pypz0+cdxzpz*pzpz0
     &               +cdxxpx*pxpx0+cdxxpy*pxpy0+cdxxpz*pxpz0
      factory=factory+cdxypx*pxpx0+cdyypx*pypx0+cdyzpx*pzpx0
     &               +cdxypy*pxpy0+cdyypy*pypy0+cdyzpy*pzpy0
     &               +cdxypz*pxpz0+cdyypz*pypz0+cdyzpz*pzpz0
     &               +cdyypx*pypx0+cdyypy*pypy0+cdyypz*pypz0
      factorz=factorz+cdxzpx*pxpx0+cdyzpx*pypx0+cdzzpx*pzpx0
     &               +cdxzpy*pxpy0+cdyzpy*pypy0+cdzzpy*pzpy0
     &               +cdxzpz*pxpz0+cdyzpz*pypz0+cdzzpz*pzpz0
     &               +cdzzpx*pzpx0+cdzzpy*pzpy0+cdzzpz*pzpz0

      factorx=factorx-zz*(cdxxpx*fxxxpx0+cdxxpy*fxxxpy0+cdxxpz*fxxxpz0
     &                   +cdxypx*fxxypx0+cdxypy*fxxypy0+cdxypz*fxxypz0
     &                   +cdxzpx*fxxzpx0+cdxzpy*fxxzpy0+cdxzpz*fxxzpz0
     &                   +cdyypx*fxyypx0+cdyypy*fxyypy0+cdyypz*fxyypz0
     &                   +cdyzpx*fxyzpx0+cdyzpy*fxyzpy0+cdyzpz*fxyzpz0
     &                   +cdzzpx*fxzzpx0+cdzzpy*fxzzpy0+cdzzpz*fxzzpz0)
      factory=factory-zz*(cdxxpx*fxxypx0+cdxxpy*fxxypy0+cdxxpz*fxxypz0
     &                   +cdxypx*fxyypx0+cdxypy*fxyypy0+cdxypz*fxyypz0
     &                   +cdxzpx*fxyzpx0+cdxzpy*fxyzpy0+cdxzpz*fxyzpz0
     &                   +cdyypx*fyyypx0+cdyypy*fyyypy0+cdyypz*fyyypz0
     &                   +cdyzpx*fyyzpx0+cdyzpy*fyyzpy0+cdyzpz*fyyzpz0
     &                   +cdzzpx*fyzzpx0+cdzzpy*fyzzpy0+cdzzpz*fyzzpz0)
      factorz=factorz-zz*(cdxxpx*fxxzpx0+cdxxpy*fxxzpy0+cdxxpz*fxxzpz0
     &                   +cdxypx*fxyzpx0+cdxypy*fxyzpy0+cdxypz*fxyzpz0
     &                   +cdxzpx*fxzzpx0+cdxzpy*fxzzpy0+cdxzpz*fxzzpz0
     &                   +cdyypx*fyyzpx0+cdyypy*fyyzpy0+cdyypz*fyyzpz0
     &                   +cdyzpx*fyzzpx0+cdyzpy*fyzzpy0+cdyzpz*fyzzpz0
     &                   +cdzzpx*fzzzpx0+cdzzpy*fzzzpy0+cdzzpz*fzzzpz0)

      factorx=factorx+cpxdxx*sdxx0+cpxdxy*sdxy0+cpxdxz*sdxz0
     &               +cpxdyy*sdyy0+cpxdyz*sdyz0+cpxdzz*sdzz0
      factory=factory+cpydxx*sdxx0+cpydxy*sdxy0+cpydxz*sdxz0
     &               +cpydyy*sdyy0+cpydyz*sdyz0+cpydzz*sdzz0
      factorz=factorz+cpzdxx*sdxx0+cpzdxy*sdxy0+cpzdxz*sdxz0
     &               +cpzdyy*sdyy0+cpzdyz*sdyz0+cpzdzz*sdzz0

      factorx=factorx-zz*(cpxdxx*dxxdxx0+cpxdxy*dxxdxy0+cpxdxz*dxxdxz0
     &                   +cpxdyy*dxxdyy0+cpxdyz*dxxdyz0+cpxdzz*dxxdzz0
     &                   +cpydxx*dxydxx0+cpydxy*dxydxy0+cpydxz*dxydxz0
     &                   +cpydyy*dxydyy0+cpydyz*dxydyz0+cpydzz*dxydzz0
     &                   +cpzdxx*dxzdxx0+cpzdxy*dxzdxy0+cpzdxz*dxzdxz0
     &                   +cpzdyy*dxzdyy0+cpzdyz*dxzdyz0+cpzdzz*dxzdzz0)
      factory=factory-zz*(cpxdxx*dxydxx0+cpxdxy*dxydxy0+cpxdxz*dxydxz0
     &                   +cpxdyy*dxydyy0+cpxdyz*dxydyz0+cpxdzz*dxydzz0
     &                   +cpydxx*dyydxx0+cpydxy*dyydxy0+cpydxz*dyydxz0
     &                   +cpydyy*dyydyy0+cpydyz*dyydyz0+cpydzz*dyydzz0
     &                   +cpzdxx*dyzdxx0+cpzdxy*dyzdxy0+cpzdxz*dyzdxz0
     &                   +cpzdyy*dyzdyy0+cpzdyz*dyzdyz0+cpzdzz*dyzdzz0)
      factorz=factorz-zz*(cpxdxx*dxzdxx0+cpxdxy*dxzdxy0+cpxdxz*dxzdxz0
     &                   +cpxdyy*dxzdyy0+cpxdyz*dxzdyz0+cpxdzz*dxzdzz0
     &                   +cpydxx*dyzdxx0+cpydxy*dyzdxy0+cpydxz*dyzdxz0
     &                   +cpydyy*dyzdyy0+cpydyz*dyzdyz0+cpydzz*dyzdzz0
     &                   +cpzdxx*dzzdxx0+cpzdxy*dzzdxy0+cpzdxz*dzzdxz0
     &                   +cpzdyy*dzzdyy0+cpzdyz*dzzdyz0+cpzdzz*dzzdzz0)

      factorx=factorx+cdxxdxx*pxdxx0+cdxxdxy*pxdxy0+cdxxdxz*pxdxz0
     &               +cdxxdyy*pxdyy0+cdxxdyz*pxdyz0+cdxxdzz*pxdzz0
     &               +cdxydxx*pydxx0+cdxydxy*pydxy0+cdxydxz*pydxz0
     &               +cdxydyy*pydyy0+cdxydyz*pydyz0+cdxydzz*pydzz0
     &               +cdxzdxx*pzdxx0+cdxzdxy*pzdxy0+cdxzdxz*pzdxz0
     &               +cdxzdyy*pzdyy0+cdxzdyz*pzdyz0+cdxzdzz*pzdzz0
     &               +cdxxdxx*pxdxx0+cdxxdxy*pxdxy0+cdxxdxz*pxdxz0
     &               +cdxxdyy*pxdyy0+cdxxdyz*pxdyz0+cdxxdzz*pxdzz0
      factory=factory+cdxydxx*pxdxx0+cdxydxy*pxdxy0+cdxydxz*pxdxz0
     &               +cdxydyy*pxdyy0+cdxydyz*pxdyz0+cdxydzz*pxdzz0
     &               +cdyydxx*pydxx0+cdyydxy*pydxy0+cdyydxz*pydxz0
     &               +cdyydyy*pydyy0+cdyydyz*pydyz0+cdyydzz*pydzz0
     &               +cdyzdxx*pzdxx0+cdyzdxy*pzdxy0+cdyzdxz*pzdxz0
     &               +cdyzdyy*pzdyy0+cdyzdyz*pzdyz0+cdyzdzz*pzdzz0
     &               +cdyydxx*pydxx0+cdyydxy*pydxy0+cdyydxz*pydxz0
     &               +cdyydyy*pydyy0+cdyydyz*pydyz0+cdyydzz*pydzz0
      factorz=factorz+cdxzdxx*pxdxx0+cdxzdxy*pxdxy0+cdxzdxz*pxdxz0
     &               +cdxzdyy*pxdyy0+cdxzdyz*pxdyz0+cdxzdzz*pxdzz0
     &               +cdyzdxx*pydxx0+cdyzdxy*pydxy0+cdyzdxz*pydxz0
     &               +cdyzdyy*pydyy0+cdyzdyz*pydyz0+cdyzdzz*pydzz0
     &               +cdzzdxx*pzdxx0+cdzzdxy*pzdxy0+cdzzdxz*pzdxz0
     &               +cdzzdyy*pzdyy0+cdzzdyz*pzdyz0+cdzzdzz*pzdzz0
     &               +cdzzdxx*pzdxx0+cdzzdxy*pzdxy0+cdzzdxz*pzdxz0
     &               +cdzzdyy*pzdyy0+cdzzdyz*pzdyz0+cdzzdzz*pzdzz0

      factorx=factorx-zz*(cdxxdxx*fxxxdxx0+cdxxdxy*fxxxdxy0
     &                   +cdxxdxz*fxxxdxz0+cdxxdyy*fxxxdyy0
     &                   +cdxxdyz*fxxxdyz0+cdxxdzz*fxxxdzz0
     &                   +cdxydxx*fxxydxx0+cdxydxy*fxxydxy0
     &                   +cdxydxz*fxxydxz0+cdxydyy*fxxydyy0
     &                   +cdxydyz*fxxydyz0+cdxydzz*fxxydzz0
     &                   +cdxzdxx*fxxzdxx0+cdxzdxy*fxxzdxy0
     &                   +cdxzdxz*fxxzdxz0+cdxzdyy*fxxzdyy0
     &                   +cdxzdyz*fxxzdyz0+cdxzdzz*fxxzdzz0
     &                   +cdyydxx*fxyydxx0+cdyydxy*fxyydxy0
     &                   +cdyydxz*fxyydxz0+cdyydyy*fxyydyy0
     &                   +cdyydyz*fxyydyz0+cdyydzz*fxyydzz0
     &                   +cdyzdxx*fxyzdxx0+cdyzdxy*fxyzdxy0
     &                   +cdyzdxz*fxyzdxz0+cdyzdyy*fxyzdyy0
     &                   +cdyzdyz*fxyzdyz0+cdyzdzz*fxyzdzz0
     &                   +cdzzdxx*fxzzdxx0+cdzzdxy*fxzzdxy0
     &                   +cdzzdxz*fxzzdxz0+cdzzdyy*fxzzdyy0
     &                   +cdzzdyz*fxzzdyz0+cdzzdzz*fxzzdzz0)
      factory=factory-zz*(cdxxdxx*fxxydxx0+cdxxdxy*fxxydxy0
     &                   +cdxxdxz*fxxydxz0+cdxxdyy*fxxydyy0
     &                   +cdxxdyz*fxxydyz0+cdxxdzz*fxxydzz0
     &                   +cdxydxx*fxyydxx0+cdxydxy*fxyydxy0
     &                   +cdxydxz*fxyydxz0+cdxydyy*fxyydyy0
     &                   +cdxydyz*fxyydyz0+cdxydzz*fxyydzz0
     &                   +cdxzdxx*fxyzdxx0+cdxzdxy*fxyzdxy0
     &                   +cdxzdxz*fxyzdxz0+cdxzdyy*fxyzdyy0
     &                   +cdxzdyz*fxyzdyz0+cdxzdzz*fxyzdzz0
     &                   +cdyydxx*fyyydxx0+cdyydxy*fyyydxy0
     &                   +cdyydxz*fyyydxz0+cdyydyy*fyyydyy0
     &                   +cdyydyz*fyyydyz0+cdyydzz*fyyydzz0
     &                   +cdyzdxx*fyyzdxx0+cdyzdxy*fyyzdxy0
     &                   +cdyzdxz*fyyzdxz0+cdyzdyy*fyyzdyy0
     &                   +cdyzdyz*fyyzdyz0+cdyzdzz*fyyzdzz0
     &                   +cdzzdxx*fyzzdxx0+cdzzdxy*fyzzdxy0
     &                   +cdzzdxz*fyzzdxz0+cdzzdyy*fyzzdyy0
     &                   +cdzzdyz*fyzzdyz0+cdzzdzz*fyzzdzz0)
      factorz=factorz-zz*(cdxxdxx*fxxzdxx0+cdxxdxy*fxxzdxy0
     &                   +cdxxdxz*fxxzdxz0+cdxxdyy*fxxzdyy0
     &                   +cdxxdyz*fxxzdyz0+cdxxdzz*fxxzdzz0
     &                   +cdxydxx*fxyzdxx0+cdxydxy*fxyzdxy0
     &                   +cdxydxz*fxyzdxz0+cdxydyy*fxyzdyy0
     &                   +cdxydyz*fxyzdyz0+cdxydzz*fxyzdzz0
     &                   +cdxzdxx*fxzzdxx0+cdxzdxy*fxzzdxy0
     &                   +cdxzdxz*fxzzdxz0+cdxzdyy*fxzzdyy0
     &                   +cdxzdyz*fxzzdyz0+cdxzdzz*fxzzdzz0
     &                   +cdyydxx*fyyzdxx0+cdyydxy*fyyzdxy0
     &                   +cdyydxz*fyyzdxz0+cdyydyy*fyyzdyy0
     &                   +cdyydyz*fyyzdyz0+cdyydzz*fyyzdzz0
     &                   +cdyzdxx*fyzzdxx0+cdyzdxy*fyzzdxy0
     &                   +cdyzdxz*fyzzdxz0+cdyzdyy*fyzzdyy0
     &                   +cdyzdyz*fyzzdyz0+cdyzdzz*fyzzdzz0
     &                   +cdzzdxx*fzzzdxx0+cdzzdxy*fzzzdxy0
     &                   +cdzzdxz*fzzzdxz0+cdzzdyy*fzzzdyy0
     &                   +cdzzdyz*fzzzdyz0+cdzzdzz*fzzzdzz0)

      gradient(1,nc1)=gradient(1,nc1)+factorx
      gradient(2,nc1)=gradient(2,nc1)+factory
      gradient(3,nc1)=gradient(3,nc1)+factorz

      gradient(1,nc2)=gradient(1,nc2)-factorx
      gradient(2,nc2)=gradient(2,nc2)-factory
      gradient(3,nc2)=gradient(3,nc2)-factorz

 1001 continue

      return
      end
