      use grid_module
      use confs_module
      use orbitals_module
      use eval_module
      use basins_module
      implicit none
      integer ier,level
      integer ksave, i,j,ilast,imol,mu
      double precision xp,yp,zp,rho,drho(3),d2rho(6),elfd,d2elf(6)
      double precision grad(3),grad_norm,xold,yold,zold,gelf(3)
      double precision delf
      double precision aux1(3,3), aux2(3,3), dx, dy,dz, valp,delt
      double precision eps,fret,gtol,xpt(3)
      double precision coeff(4),root(3),rooti(3)
      double precision, dimension(:,:), allocatable :: csave
      integer, dimension(:), allocatable :: orbsym
      character*1 blank
      character*40 fileelf,dumname,filein,filegam,filebas,molecule
      character*1 yes,ok
      character*3 func_name
      logical func_type,ifpi,rho_and_elf
      external func,func2
      common /files/ molecule,filein,filegam,fileelf,filebas,rho_and_elf
      common/ftype/func_type
      data blank/' '/, yes/'y'/, eps/1.0d-10/
      ksave=0
      call read_wfn
      allocate(indij(nmol*(nmol+1)/2))
      allocate(g(ngto),gx(ngto),gy(ngto),gz(ngto),extest(ngto))
      allocate(fxx(ngto),fxy(ngto),fxz(ngto),fyy(ngto),fyz(ngto))
      allocate(fzz(ngto),fxxx(ngto),fxxy(ngto),fxxz(ngto),fxyy(ngto))
      allocate(fxyz(ngto),fxzz(ngto),fyyy(ngto),fyyz(ngto),fyzz(ngto))
      allocate(fzzz(ngto),csave(nmol,ngto),orbsym(nmol))
      allocate(phi(nmol),dphi(nmol,3),d2phi(nmol,6),d3phi(nmol,10))
      allocate(rhorb(nmol*(nmol+1)/2),drhorb(nmol*(nmol+1)/2,3))
      allocate(d2rhorb(nmol*(nmol+1)/2,6),dotorb(nmol*(nmol+1)/2,3))
      allocate(ddotorb(nmol*(nmol+1)/2,9),d2dotorb(nmol*(nmol+1)/2,18))
      call read_elf2
      do i=1,nmol
         do j=1,ngto
            csave(i,j)=c(i,j)
         enddo
      enddo
      ifpi=.false.
      write(*,*) 'total elf only?'
      read(*,'(a)') ok
      if(ok.ne.'y') then
         write(*,*) 'input orbital symmetry : 0 sigma, 1 pi'
         read(*,*) (orbsym(i),i=1,nmol)
         ifpi=.true.
      endif
      do while (ksave.ge.0)
      write(*,*) 'function elf/rho '
      read(*,'(a3)') func_name
      if(func_name.eq.'elf') then
        func_type=.true.
      else
        func_type=.false.
      endif
      write(*,*) ' type of critical point 1 attractor, 0 otherwise'
     .           ,',  <0 fin '
      read(*,*) ksave
      if(ksave.ge.0) then
      write(*,*) ' coordinates'
      read(*,*) xold,yold,zold
      ilast=1
      if(ifpi) then
         ilast=3
      endif
      do i=1,ilast
      if(i.eq.1) then
         write(*,*) 'total wavefunction '
         do imol=1,nmol
            do mu=1,ngto
               c(imol,mu)=csave(imol,mu)
            enddo
         enddo
      elseif (i.eq.2) then
         write(*,*) 'sigma restricted wavefunction '
         do imol=1,nmol
            do mu=1,ngto
               c(imol,mu)=csave(imol,mu)*(1.0d0-dfloat(orbsym(imol)))
            enddo
         enddo
      elseif (i.eq.3) then
         write(*,*) 'pi restricted wavefunction '
         do imol=1,nmol
            do mu=1,ngto
               c(imol,mu)=csave(imol,mu)*dfloat(orbsym(imol))
            enddo
         enddo
      endif
      xp=xold
      yp=yold
      zp=zold
      level=0
      call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
      write(*,*) 'elf =',elfd, ' rho = ',rho
      gtol=1.0d-6
      xpt(1)=xp
      xpt(2)=yp
      xpt(3)=zp
      if(ksave.eq.0) call dfmfp(func2,xpt,fret,eps,ier)
      if(ksave.eq.1) call dfmfp(func,xpt,fret,eps,ier)
      delt=0.2d0
      delf=0.0d0
      xp=xpt(1)
      yp=xpt(2)
      zp=xpt(3)
      grad_norm=1.0d0
      do while(grad_norm.ge.1.0d-12)
         level=2
         call elf_eval(rho,drho,d2rho,elfd,gelf,d2elf,xp,yp,zp,level)
         if(func_type) then
            valp=elfd
            grad(1)=gelf(1)
            grad(2)=gelf(2)
            grad(3)=gelf(3)
            aux1(1,1)=d2elf(1)
            aux1(1,2)=d2elf(2)
            aux1(1,3)=d2elf(3)
            aux1(2,1)=d2elf(2)
            aux1(2,2)=d2elf(4)
            aux1(2,3)=d2elf(5)
            aux1(3,1)=d2elf(3)
            aux1(3,2)=d2elf(5)
            aux1(3,3)=d2elf(6)
         else
            valp=rho
            grad(1)=drho(1)
            grad(2)=drho(2)
            grad(3)=drho(3)
            aux1(1,1)=d2rho(1)
            aux1(1,2)=d2rho(2)
            aux1(1,3)=d2rho(3)
            aux1(2,1)=d2rho(2)
            aux1(2,2)=d2rho(4)
            aux1(2,3)=d2rho(5)
            aux1(3,1)=d2rho(3)
            aux1(3,2)=d2rho(5)
            aux1(3,3)=d2rho(6)
         endif
         grad_norm=grad(1)**2+grad(2)**2+grad(3)**2
         if(ksave.eq.1) then
         call minv3(aux1,aux2)
            dx=aux2(1,1)*grad(1)+aux2(1,2)*grad(2)+aux2(1,3)*grad(3)
            dy=aux2(2,1)*grad(1)+aux2(2,2)*grad(2)+aux2(2,3)*grad(3)
            dz=aux2(3,1)*grad(1)+aux2(3,2)*grad(2)+aux2(3,3)*grad(3)
            xp=xp-dx
            yp=yp-dy
            zp=zp-dz
         else
            grad_norm=0.0d0
         endif
      enddo
      coeff(1)=1.0d0
      coeff(2)=-(aux1(1,1)+aux1(2,2)+aux1(3,3))
      coeff(3)=aux1(1,1)*aux1(2,2)+aux1(1,1)*aux1(3,3)+
     .         aux1(2,2)*aux1(3,3)-aux1(1,2)*aux1(1,2)-
     .         aux1(1,3)*aux1(1,3)-aux1(2,3)*aux1(2,3)
      coeff(4)=aux1(1,1)*aux1(2,3)*aux1(2,3)+
     .         aux1(2,2)*aux1(1,3)*aux1(1,3)+
     .         aux1(3,3)*aux1(1,2)*aux1(1,2)-
     .         aux1(1,1)*aux1(2,2)*aux1(3,3)-
     .   2.0d0*aux1(1,2)*aux1(1,3)*aux1(2,3)
      call cubic(coeff,root,rooti)
      write(*,'(t12,"coordinates",t30,"elf",t40,"rho",t48,
     .              "norm of grad",t61,"nabla2 rho")')
      write(*,'(4f8.3,3f12.4)') xp,yp,zp,elfd,rho,grad_norm,
     .      d2rho(1)+d2rho(4)+d2rho(6)
      write(*,'(t8,"hessian matrix eigenvalues")')
      write(*,'(3d20.12)') root
      enddo
      endif
      enddo
      end
