      subroutine generf2
      use grid_module
      use orbitals_module
      implicit none
c     double precision origin(3)
c     double precision step(3)
      double precision x, xa, xa2, arg, expo, tipo
      integer icoord,i,mu
      origin(1)=xmin
      origin(2)=ymin
      origin(3)=zmin
      step(1)=(xmax-xmin)/(np(1)-1)
      step(2)=(ymax-ymin)/(np(2)-1)
      step(3)=(zmax-zmin)/(np(3)-1)
      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)
           tipo=ir(1,mu)
           if(tipo.eq.3) then
             fx(i,mu)=xa*xa2*expo
             dfx(i,mu)=(3.0d0+arg+arg)*xa2*expo
             d2fx(i,mu)=2.0d0*(3.0d0+arg*(7.0d0+arg+arg))*xa*expo
           elseif(tipo.eq.2) then
             fx(i,mu)=xa2*expo
             dfx(i,mu)=(2.0d0+arg+arg)*xa*expo
             d2fx(i,mu)=2.0d0*(1.0d0+arg*(5.0d0+arg+arg))*expo
           elseif(tipo.eq.1) then
             fx(i,mu)=xa*expo
             dfx(i,mu)=(1.0d0+arg+arg)*expo
             d2fx(i,mu)=-2.0d0*alfa(mu)*xa*(arg+arg+3.0d0)*expo
           elseif(tipo.eq.0) then
             fx(i,mu)=expo
             dfx(i,mu)=-2.0d0*xa*alfa(mu)*expo
             d2fx(i,mu)=-2.0d0*alfa(mu)*(arg+arg+1.0d0)*expo
           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)
           tipo=ir(2,mu)
           if(tipo.eq.3) then
             fy(i,mu)=xa*xa2*expo
             dfy(i,mu)=(3.0d0+arg+arg)*xa2*expo
             d2fy(i,mu)=2.0d0*(3.0d0+arg*(7.0d0+arg+arg))*xa*expo
           elseif(tipo.eq.2) then
             fy(i,mu)=xa2*expo
             dfy(i,mu)=(2.0d0+arg+arg)*xa*expo
             d2fy(i,mu)=2.0d0*(1.0d0+arg*(5.0d0+arg+arg))*expo
           elseif(tipo.eq.1) then
             fy(i,mu)=xa*expo
             dfy(i,mu)=(1.0d0+arg+arg)*expo
             d2fy(i,mu)=-2.0d0*alfa(mu)*xa*(arg+arg+3.0d0)*expo
           elseif(tipo.eq.0) then
             fy(i,mu)=expo
             dfy(i,mu)=-2.0d0*xa*alfa(mu)*expo
             d2fy(i,mu)=-2.0d0*alfa(mu)*(arg+arg+1.0d0)*expo
           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)
           tipo=ir(3,mu)
           if(tipo.eq.3) then
             fz(i,mu)=xa*xa2*expo
             dfz(i,mu)=(3.0d0+arg+arg)*xa2*expo
             d2fz(i,mu)=2.0d0*(3.0d0+arg*(7.0d0+arg+arg))*xa*expo
           elseif(tipo.eq.2) then
             fz(i,mu)=xa2*expo
             dfz(i,mu)=(2.0d0+arg+arg)*xa*expo
             d2fz(i,mu)=2.0d0*(1.0d0+arg*(5.0d0+arg+arg))*expo
           elseif(tipo.eq.1) then
             fz(i,mu)=xa*expo
             dfz(i,mu)=(1.0d0+arg+arg)*expo
             d2fz(i,mu)=-2.0d0*alfa(mu)*xa*(arg+arg+3.0d0)*expo
           elseif(tipo.eq.0) then
             fz(i,mu)=expo
             dfz(i,mu)=-2.0d0*xa*alfa(mu)*expo
             d2fz(i,mu)=-2.0d0*alfa(mu)*(arg+arg+1.0d0)*expo
           endif
        enddo
      enddo
      return
      end
