      subroutine generf
      use orbitals_module
      use grid_module
      implicit none
      integer icoord,i,mu
      double precision x, xa, xa2, arg, expo, type
      x=origin(1)-step(1)
      do i=1,np(1)
        x=x+step(1)
        do mu=1,ngto
           xa=x-rc(1,mu)
           xa2=xa*xa
           arg=-alfa(mu)*xa2
           expo=dexp(arg)
           type=ir(1,mu)
           if(type.eq.3) then
             fx(i,mu)=xa*xa2*expo
             dfx(i,mu)=(3.0d0+arg+arg)*xa2*expo
             dfx2(i,mu)=xa*expo*(6.0d0+14.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.2) then
             fx(i,mu)=xa2*expo
             dfx(i,mu)=(2.0d0+arg+arg)*xa*expo
             dfx2(i,mu)=expo*(2.0d0+10.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.1) then
             fx(i,mu)=xa*expo
             dfx(i,mu)=(1.0d0+arg+arg)*expo
             dfx2(i,mu)=-2*expo*xa*alfa(mu)*(3.0d0+2.0d0*arg)
           elseif(type.eq.0) then
             fx(i,mu)=expo
             dfx(i,mu)=-2.0d0*xa*alfa(mu)*expo
             dfx2(i,mu)=-2.0d0*expo*alfa(mu)*(2.0d0*arg+1.0d0)
           endif
        enddo
      enddo
      x=origin(2)-step(2)
      do i=1,np(2)
        x=x+step(2)
        do mu=1,ngto
           xa=x-rc(2,mu)
           xa2=xa*xa
           arg=-alfa(mu)*xa2
           expo=dexp(arg)
           type=ir(2,mu)
           if(type.eq.3) then
             fy(i,mu)=xa*xa2*expo
             dfy(i,mu)=(3.0d0+arg+arg)*xa2*expo
             dfy2(i,mu)=xa*expo*(6.0d0+14.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.2) then
             fy(i,mu)=xa2*expo
             dfy(i,mu)=(2.0d0+arg+arg)*xa*expo
             dfy2(i,mu)=expo*(2.0d0+10.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.1) then
             fy(i,mu)=xa*expo
             dfy(i,mu)=(1.0d0+arg+arg)*expo
             dfy2(i,mu)=-2*expo*xa*alfa(mu)*(3.0d0+2.0d0*arg)
           elseif(type.eq.0) then
             fy(i,mu)=expo
             dfy(i,mu)=-2.0d0*xa*alfa(mu)*expo
             dfy2(i,mu)=-2.0d0*expo*alfa(mu)*(2.0d0*arg+1.0d0)
           endif
        enddo
      enddo
      x=origin(3)-step(3)
      do i=1,np(3)
        x=x+step(3)
        do mu=1,ngto
           xa=x-rc(3,mu)
           xa2=xa*xa
           arg=-alfa(mu)*xa2
           expo=dexp(arg)
           type=ir(3,mu)
           if(type.eq.3) then
             fz(i,mu)=xa*xa2*expo
             dfz(i,mu)=(3.0d0+arg+arg)*xa2*expo
             dfz2(i,mu)=xa*expo*(6.0d0+14.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.2) then
             fz(i,mu)=xa2*expo
             dfz(i,mu)=(2.0d0+arg+arg)*xa*expo
             dfz2(i,mu)=expo*(2.0d0+10.0d0*arg+4.0d0*arg*arg)
           elseif(type.eq.1) then
             fz(i,mu)=xa*expo
             dfz(i,mu)=(1.0d0+arg+arg)*expo
             dfz2(i,mu)=-2*expo*xa*alfa(mu)*(3.0d0+2.0d0*arg)
           elseif(type.eq.0) then
             fz(i,mu)=expo
             dfz(i,mu)=-2.0d0*xa*alfa(mu)*expo
             dfz2(i,mu)=-2.0d0*expo*alfa(mu)*(2.0d0*arg+1.0d0)
           endif
        enddo
      enddo
      return
      end
