      subroutine elf_eval(rho,drho,d2rho,elfd,delf,d2elf,xp,yp,zp,level)
      use confs_module
      use orbitals_module
      use eval_module
      implicit none
      double precision rha,rhb,gradfa,gradfb,gradf,dgradf(3)
      double precision drha(3),drhb(3),d2rha(6),d2rhb(6),d2gradf(6)
      double precision dgradfa(3),dgradfb(3),d2gradfa(6),d2gradfb(6)
      double precision dratio(3),d2ratio(6),p163
      double precision ratio,x,y,z,r2
      double precision x0,x1,x2,x3,y0,y1,y2,y3,z0,z1,z2,z3
      double precision alfa2,alfax,alfax2,alfay,alfay2,alfaz,alfaz2
      double precision zp,yp,xp,ex,fermi,tol
      double precision bold(3),cold(3)
      double precision rho,elfd,drho(3),delf(3),d2rho(6),d2elf(6)
      integer i,j,k,l,mu,level
      integer ij,kl,imol,jmol,lmol,kmol
      common/threshold/ tol
      data fermi/0.121300565d0/
      p163=-16.0d0/3.0d0
      tol=2.871d-05
      if(ifci) then
         nmola=nmol
      endif
      do i=1,nmola
         indij(i)=(i-1)*(nmola+nmola-i)/2
      enddo
      do mu=1, ngto
         x=xp-rc(1,mu)
         y=yp-rc(2,mu)
         z=zp-rc(3,mu)
         r2=x*x+y*y+z*z
         ex=dexp(-alfa(mu)*r2)
         alfa2=-alfa(mu)*2.0d0
         alfax=alfa2*x
         alfax2=alfax*x
         alfay=alfa2*y
         alfay2=alfay*y
         alfaz=alfa2*z
         alfaz2=alfaz*z
         i=ir(1,mu)
         j=ir(2,mu)
         k=ir(3,mu)
         if(level.ge.0) then
           if(i.eq.0) then
             x0=1.0d0
             x1=alfax
             x2=alfa2*(1.0d0+alfax2)
           else if (i.eq.1) then
             x0=x
             x1=1.0d0+alfax2
             x2=alfax*(3.d0+alfax2)
           else if (i.eq.2) then
             x0=x*x
             x1=x*(2.d0+alfax2)
             x2=2.0d0+alfax2*(5.0d0+alfax2)
           else if (i.eq.3) then
             x0=x**3
             x1=x*x*(3.0d0+alfax2)
             x2=x*(6.0d0+alfax2*(7.d0+alfax2))
           else if (i.gt.3) then
             x0=x**i
             x1=x**(i-1)*(dfloat(i)+alfax2)
             x2=x**(i-2)*(dfloat(i)*dfloat(i-1)+alfax2*
     .         (dfloat(i+i+1)+alfax2))
           endif
           if(j.eq.0) then
             y0=1.0d0
             y1=alfay
             y2=alfa2*(1.0d0+alfay2)
           else if (j.eq.1) then
             y0=y
             y1=1.0d0+alfay2
             y2=alfay*(3.d0+alfay2)
           else if (j.eq.2) then
             y0=y*y
             y1=y*(2.d0+alfay2)
             y2=2.0d0+alfay2*(5.0d0+alfay2)
           else if (j.eq.3) then
             y0=y**3
             y1=y*y*(3.0d0+alfay2)
             y2=y*(6.0d0+alfay2*(7.d0+alfay2))
           else if (j.gt.3) then
             y0=y**j
             y1=y**(j-1)*(dfloat(j)+alfay2)
             y2=y**(j-2)*(dfloat(j)*dfloat(j-1)+alfay2*
     .         (dfloat(j+j+1)+alfay2))
           endif
           if(k.eq.0) then
             z0=1.0d0
             z1=alfaz
             z2=alfa2*(1.0d0+alfaz2)
           else if (k.eq.1) then
             z0=z
             z1=1.0d0+alfaz2
             z2=alfaz*(3.d0+alfaz2)
           else if (k.eq.2) then
             z0=z*z
             z1=z*(2.d0+alfaz2)
             z2=2.0d0+alfaz2*(5.0d0+alfaz2)
           else if (k.eq.3) then
             z0=z**3
             z1=z*z*(3.0d0+alfaz2)
             z2=z*(6.0d0+alfaz2*(7.d0+alfaz2))
           else if (k.gt.3) then
             z0=z**k
             z1=z**(k-1)*(dfloat(k)+alfaz2)
             z2=z**(k-2)*(dfloat(k)*dfloat(k-1)+alfaz2*
     .         (dfloat(k+k+1)+alfaz2))
           endif
         endif
         if(level.ge.1) then
           if(i.eq.0) then
             x3=alfa2*alfax*(3.0d0+alfax2)
           else if (i.eq.1) then
             x3=alfa2*(3.0d0+alfax2*(6.0d0+alfax2))
           else if (i.eq.2) then
             x3=alfax*(1.2d1+alfax2*(9.0d0+alfax2))
           else if (i.eq.3) then
             x3=6.0d0+alfax2*(2.7d1+alfax2*(1.2d1+alfax2))
           else if (i.gt.3) then
             x3=x**(i-3)*(dfloat(i*(i-1)*(i-2))+alfax2*
     .          (3.0d0*dfloat(i*i)+alfax2*(dfloat(3*i+3)+alfax2)))
           endif
           if(j.eq.0) then
             y3=alfa2*alfay*(3.0d0+alfay2)
           else if (j.eq.1) then
             y3=alfa2*(3.0d0+alfay2*(6.0d0+alfay2))
           else if (j.eq.2) then
             y3=alfay*(1.2d1+alfay2*(9.0d0+alfay2))
           else if (j.eq.3) then
             y3=6.0d0+alfay2*(2.7d1+alfay2*(1.2d1+alfay2))
           else if (j.gt.3) then
             y3=y**(j-3)*(dfloat(j*(j-1)*(j-2))+alfay2*
     .          (3.0d0*dfloat(j*j)+alfay2*(dfloat(3*j+3)+alfay2)))
           endif
           if(k.eq.0) then
             z3=alfa2*alfaz*(3.0d0+alfaz2)
           else if (k.eq.1) then
             z3=alfa2*(3.0d0+alfaz2*(6.0d0+alfaz2))
           else if (k.eq.2) then
             z3=alfaz*(1.2d1+alfaz2*(9.0d0+alfaz2))
           else if (k.eq.3) then
             z3=6.0d0+alfaz2*(2.7d1+alfaz2*(1.2d1+alfaz2))
           else if (k.gt.3) then
             z3=z**(k-3)*(dfloat(k*(k-1)*(k-2))+alfaz2*
     .          (3.0d0*dfloat(k*k)+alfaz2*(dfloat(3*k+3)+alfaz2)))
           endif
         endif
         if (dabs(ex).lt.1.0d-14) then
            extest(mu)=0
         else
            extest(mu)=1
         endif
         if(level.ge.0) then
           g(mu)=x0*y0*z0*ex
           gx(mu)=x1*y0*z0*ex
           gy(mu)=x0*y1*z0*ex
           gz(mu)=x0*y0*z1*ex
         endif
         if(level.ge.1) then
           fxx(mu)=x2*y0*z0*ex
           fxy(mu)=x1*y1*z0*ex
           fxz(mu)=x1*y0*z1*ex
           fyy(mu)=x0*y2*z0*ex
           fyz(mu)=x0*y1*z1*ex
           fzz(mu)=x0*y0*z2*ex
         endif
         if(level.ge.2) then
           fxxx(mu)=x3*y0*z0*ex
           fxxy(mu)=x2*y1*z0*ex
           fxxz(mu)=x2*y0*z1*ex
           fxyy(mu)=x1*y2*z0*ex
           fxyz(mu)=x1*y1*z1*ex
           fxzz(mu)=x1*y0*z2*ex
           fyyy(mu)=x0*y3*z0*ex
           fyyz(mu)=x0*y2*z1*ex
           fyzz(mu)=x0*y1*z2*ex
           fzzz(mu)=x0*y0*z3*ex
         endif
      enddo
*              ****
*              transform to molecular orbitals
*              ****
      do imol=1,nmol
         phi(imol)=0.0d0
         dphi(imol,1)=0.0d0
         dphi(imol,2)=0.0d0
         dphi(imol,3)=0.0d0
         do mu=1,ngto
            phi(imol)=phi(imol)+c(imol,mu)*g(mu)
            dphi(imol,1)=dphi(imol,1)+c(imol,mu)*gx(mu)
            dphi(imol,2)=dphi(imol,2)+c(imol,mu)*gy(mu)
            dphi(imol,3)=dphi(imol,3)+c(imol,mu)*gz(mu)
         enddo
         if(level.ge.1) then
            d2phi(imol,1)=0.0d0
            d2phi(imol,2)=0.0d0
            d2phi(imol,3)=0.0d0
            d2phi(imol,4)=0.0d0
            d2phi(imol,5)=0.0d0
            d2phi(imol,6)=0.0d0
            do mu=1,ngto
               d2phi(imol,1)=d2phi(imol,1)+c(imol,mu)*fxx(mu)
               d2phi(imol,2)=d2phi(imol,2)+c(imol,mu)*fxy(mu)
               d2phi(imol,3)=d2phi(imol,3)+c(imol,mu)*fxz(mu)
               d2phi(imol,4)=d2phi(imol,4)+c(imol,mu)*fyy(mu)
               d2phi(imol,5)=d2phi(imol,5)+c(imol,mu)*fyz(mu)
               d2phi(imol,6)=d2phi(imol,6)+c(imol,mu)*fzz(mu)
            enddo
         endif
         if (level.ge.2) then
            d3phi(imol,1)=0.0d0
            d3phi(imol,2)=0.0d0
            d3phi(imol,3)=0.0d0
            d3phi(imol,4)=0.0d0
            d3phi(imol,5)=0.0d0
            d3phi(imol,6)=0.0d0
            d3phi(imol,7)=0.0d0
            d3phi(imol,8)=0.0d0
            d3phi(imol,9)=0.0d0
            d3phi(imol,10)=0.0d0
            do mu=1,ngto
               d3phi(imol,1)=d3phi(imol,1)+c(imol,mu)*fxxx(mu)
               d3phi(imol,2)=d3phi(imol,2)+c(imol,mu)*fxxy(mu)
               d3phi(imol,3)=d3phi(imol,3)+c(imol,mu)*fxxz(mu)
               d3phi(imol,4)=d3phi(imol,4)+c(imol,mu)*fxyy(mu)
               d3phi(imol,5)=d3phi(imol,5)+c(imol,mu)*fxyz(mu)
               d3phi(imol,6)=d3phi(imol,6)+c(imol,mu)*fxzz(mu)
               d3phi(imol,7)=d3phi(imol,7)+c(imol,mu)*fyyy(mu)
               d3phi(imol,8)=d3phi(imol,8)+c(imol,mu)*fyyz(mu)
               d3phi(imol,9)=d3phi(imol,9)+c(imol,mu)*fyzz(mu)
               d3phi(imol,10)=d3phi(imol,10)+c(imol,mu)*fzzz(mu)
            enddo
         endif
      enddo
      do imol=1,nmola
         do jmol=imol,nmola
            lmol=indij(imol)+jmol
            rhorb(lmol)=phi(imol)*phi(jmol)
            dotorb(lmol,1)=phi(imol)*dphi(jmol,1)- 
     .                     phi(jmol)*dphi(imol,1)
            dotorb(lmol,2)=phi(imol)*dphi(jmol,2)- 
     .                     phi(jmol)*dphi(imol,2)
            dotorb(lmol,3)=phi(imol)*dphi(jmol,3)- 
     .                     phi(jmol)*dphi(imol,3)
            if(imol.ne.jmol) then
                rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
            endif
         enddo
      enddo
****
*              loop over molecular orbitals (one-int)
      rha=0.0d0
      gradfa=0.0d0
      do l=1,na
         lmol=ia(l)
         rha=rha+pa(l)*rhorb(lmol)
      enddo
*              compute D(r)
      do l=1,naa
         imol=iaa(l)
         jmol=jaa(l)
         kmol=kaa(l)
         lmol=laa(l)
         ij=indij(imol)+jmol
         kl=indij(kmol)+lmol
         gradfa=gradfa+paa(l)*(dotorb(ij,1)*dotorb(kl,1)+
     .          dotorb(ij,2)*dotorb(kl,2)+dotorb(ij,3)*dotorb(kl,3))
      enddo
      if(nb.eq.0) then
         rhb=rha
         gradfb=gradfa
      else
**    compute beta contribution (open shell case)
      do imol=nmola+1,nmol
         do jmol=imol,nmol
            lmol=indij(imol-nmola)+jmol-nmola
            rhorb(lmol)=phi(imol)*phi(jmol)
            dotorb(lmol,1)=phi(imol)*dphi(jmol,1)-phi(jmol)*dphi(imol,1)
            dotorb(lmol,2)=phi(imol)*dphi(jmol,2)-phi(jmol)*dphi(imol,2)
            dotorb(lmol,3)=phi(imol)*dphi(jmol,3)-phi(jmol)*dphi(imol,3)
            if(imol.ne.jmol) then
                rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
            endif
         enddo
      enddo
      rhb=0.0d0
      do l=1,nb
         lmol=ib(l)
         rhb=rhb+pb(l)*rhorb(lmol)
      enddo
      gradfb=0.0d0
      do l=1,nbb
         imol=ibb(l)
         jmol=jbb(l)
         kmol=kbb(l)
         lmol=lbb(l)
         ij=indij(imol)+jmol
         kl=indij(kmol)+lmol
         gradfb=gradfb+pbb(l)*(dotorb(ij,1)*dotorb(kl,1)+
     .          dotorb(ij,2)*dotorb(kl,2)+dotorb(ij,3)*dotorb(kl,3))
      enddo
      endif
      rho=rha+rhb+1.0d-16
      gradf=(gradfa+gradfb)+tol*rho
      ratio=gradf*gradf*fermi*rho**p163
      elfd=1.d0/(1.0d0+ratio)
      if(level.ge.1) then
         do imol=1,nmola
            do jmol=imol,nmola
               lmol=indij(imol)+jmol
               rhorb(lmol)=phi(imol)*phi(jmol)
               dotorb(lmol,1)=phi(imol)*dphi(jmol,1)- 
     .                     phi(jmol)*dphi(imol,1)
               dotorb(lmol,2)=phi(imol)*dphi(jmol,2)- 
     .                     phi(jmol)*dphi(imol,2)
               dotorb(lmol,3)=phi(imol)*dphi(jmol,3)- 
     .                     phi(jmol)*dphi(imol,3)
               drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                        phi(jmol)*dphi(imol,1)
               drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                        phi(jmol)*dphi(imol,2)
               drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                        phi(jmol)*dphi(imol,3)
               ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,1)
               ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,2)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,1)*dphi(jmol,2)
               ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,3)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,1)*dphi(jmol,3)
               ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,1)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,2)*dphi(jmol,1)
               ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)- 
     .                         phi(jmol)*d2phi(imol,4)
               ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,3)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,2)*dphi(jmol,3)
               ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,1)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,3)*dphi(jmol,1)
               ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,2)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,3)*dphi(jmol,2)
               ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)- 
     .                         phi(jmol)*d2phi(imol,6)
               if(imol.ne.jmol) then
                   rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                   drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                   drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                   drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
               endif
            enddo
         enddo
         drha(1)=0.0d0
         drha(2)=0.0d0
         drha(3)=0.0d0
         do l=1,na
            lmol=ia(l)
            drha(1)=drha(1)+pa(l)*drhorb(lmol,1)
            drha(2)=drha(2)+pa(l)*drhorb(lmol,2)
            drha(3)=drha(3)+pa(l)*drhorb(lmol,3)
         enddo
         dgradfa(1)=0.0d0
         dgradfa(2)=0.0d0
         dgradfa(3)=0.0d0
         do l=1,naa
            imol=iaa(l)
            jmol=jaa(l)
            kmol=kaa(l)
            lmol=laa(l)
            ij=indij(imol)+jmol
            kl=indij(kmol)+lmol
            dgradfa(1)=dgradfa(1)+paa(l)*(dotorb(ij,1)*ddotorb(kl,1)+
     .                 ddotorb(ij,1)*dotorb(kl,1)+
     .                 ddotorb(ij,4)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,4)+
     .                 ddotorb(ij,7)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,7))
            dgradfa(2)=dgradfa(2)+paa(l)*(dotorb(ij,1)*ddotorb(kl,2)+
     .                 ddotorb(ij,2)*dotorb(kl,1)+
     .                 ddotorb(ij,5)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,5)+
     .                 ddotorb(ij,8)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,8))
            dgradfa(3)=dgradfa(3)+paa(l)*(dotorb(ij,1)*ddotorb(kl,3)+
     .                 ddotorb(ij,3)*dotorb(kl,1)+
     .                 ddotorb(ij,6)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,6)+
     .                 ddotorb(ij,9)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,9))
         enddo
         if(nb.eq.0) then
            drhb(1)=drha(1)
            drhb(2)=drha(2)
            drhb(3)=drha(3)
            dgradfb(1)=dgradfa(1)
            dgradfb(2)=dgradfa(2)
            dgradfb(3)=dgradfa(3)
         else
**    compute beta contribution (open shell case)
         do imol=nmola+1,nmol
            do jmol=imol,nmol
               lmol=indij(imol-nmola)+jmol-nmola
               rhorb(lmol)=phi(imol)*phi(jmol)
               dotorb(lmol,1)=phi(imol)*dphi(jmol,1)-
     .                        phi(jmol)*dphi(imol,1)
               dotorb(lmol,2)=phi(imol)*dphi(jmol,2)-
     .                        phi(jmol)*dphi(imol,2)
               dotorb(lmol,3)=phi(imol)*dphi(jmol,3)-
     .                        phi(jmol)*dphi(imol,3)
               drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                        phi(jmol)*dphi(imol,1)
               drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                        phi(jmol)*dphi(imol,2)
               drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                        phi(jmol)*dphi(imol,3)
               ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,1)
               ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,2)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,1)*dphi(jmol,2)
               ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,3)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,1)*dphi(jmol,3)
               ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,1)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,2)*dphi(jmol,1)
               ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)- 
     .                         phi(jmol)*d2phi(imol,4)
               ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,3)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,2)*dphi(jmol,3)
               ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,1)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,3)*dphi(jmol,1)
               ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,2)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,3)*dphi(jmol,2)
               ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)- 
     .                         phi(jmol)*d2phi(imol,6)
               if(imol.ne.jmol) then
                   rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                   drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                   drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                   drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
               endif
            enddo
         enddo
         drhb(1)=0.0d0
         drhb(2)=0.0d0
         drhb(3)=0.0d0
         do l=1,nb
            lmol=ib(l)
            drhb(1)=drhb(1)+pb(l)*drhorb(lmol,1)
            drhb(2)=drhb(2)+pb(l)*drhorb(lmol,2)
            drhb(3)=drhb(3)+pb(l)*drhorb(lmol,3)
         enddo
         dgradfb(1)=0.0d0
         dgradfb(2)=0.0d0
         dgradfb(3)=0.0d0
         do l=1,nbb
            imol=ibb(l)
            jmol=jbb(l)
            kmol=kbb(l)
            lmol=lbb(l)
            ij=indij(imol)+jmol
            kl=indij(kmol)+lmol
            dgradfb(1)=dgradfb(1)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,1)+
     .                 ddotorb(ij,1)*dotorb(kl,1)+
     .                 ddotorb(ij,4)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,4)+
     .                 ddotorb(ij,7)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,7))
            dgradfb(2)=dgradfb(2)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,2)+
     .                 ddotorb(ij,2)*dotorb(kl,1)+
     .                 ddotorb(ij,5)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,5)+
     .                 ddotorb(ij,8)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,8))
            dgradfb(3)=dgradfb(3)+pbb(l)*(dotorb(ij,1)*ddotorb(kl,3)+
     .                 ddotorb(ij,3)*dotorb(kl,1)+
     .                 ddotorb(ij,6)*dotorb(kl,2)+
     .                 dotorb(ij,2)*ddotorb(kl,6)+
     .                 ddotorb(ij,9)*dotorb(kl,3)+
     .                 dotorb(ij,3)*ddotorb(kl,9))
         enddo
         endif
         drho(1)=drha(1)+drhb(1)
         drho(2)=drha(2)+drhb(2)
         drho(3)=drha(3)+drhb(3)
         dgradf(1)=dgradfa(1)+dgradfb(1)+tol*drho(1)
         dgradf(2)=dgradfa(2)+dgradfb(2)+tol*drho(2)
         dgradf(3)=dgradfa(3)+dgradfb(3)+tol*drho(3)
         dratio(1)=2.0d0*ratio*dgradf(1)/gradf+p163*ratio*drho(1)/rho
         dratio(2)=2.0d0*ratio*dgradf(2)/gradf+p163*ratio*drho(2)/rho
         dratio(3)=2.0d0*ratio*dgradf(3)/gradf+p163*ratio*drho(3)/rho
         delf(1)=-elfd*elfd*dratio(1)
         delf(2)=-elfd*elfd*dratio(2)
         delf(3)=-elfd*elfd*dratio(3)
      endif
**** TMP
      if(level.ge.2) then
         do imol=1,nmola
            do jmol=imol,nmola
               lmol=indij(imol)+jmol
               rhorb(lmol)=phi(imol)*phi(jmol)
               dotorb(lmol,1)=phi(imol)*dphi(jmol,1)- 
     .                     phi(jmol)*dphi(imol,1)
               dotorb(lmol,2)=phi(imol)*dphi(jmol,2)- 
     .                     phi(jmol)*dphi(imol,2)
               dotorb(lmol,3)=phi(imol)*dphi(jmol,3)- 
     .                     phi(jmol)*dphi(imol,3)
               drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                        phi(jmol)*dphi(imol,1)
               drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                        phi(jmol)*dphi(imol,2)
               drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                        phi(jmol)*dphi(imol,3)
               ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,1)
               ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,2)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,1)*dphi(jmol,2)
               ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,3)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,1)*dphi(jmol,3)
               ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,1)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,2)*dphi(jmol,1)
               ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)- 
     .                         phi(jmol)*d2phi(imol,4)
               ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,3)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,2)*dphi(jmol,3)
               ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,1)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,3)*dphi(jmol,1)
               ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,2)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,3)*dphi(jmol,2)
               ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)- 
     .                         phi(jmol)*d2phi(imol,6)
               d2rhorb(lmol,1)=phi(imol)*d2phi(jmol,1)+
     .                      2.0d0*dphi(imol,1)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,1)
               d2rhorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                        dphi(imol,1)*dphi(jmol,2)+
     .                        dphi(imol,2)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,2)
               d2rhorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                        dphi(imol,1)*dphi(jmol,3)+
     .                        dphi(imol,3)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,3)
               d2rhorb(lmol,4)=phi(imol)*d2phi(jmol,4)+
     .                      2.0d0*dphi(imol,2)*dphi(jmol,2)+
     .                        phi(jmol)*d2phi(imol,4)
               d2rhorb(lmol,5)=phi(imol)*d2phi(jmol,5)+
     .                        dphi(imol,2)*dphi(jmol,3)+
     .                        dphi(imol,3)*dphi(jmol,2)+
     .                        phi(jmol)*d2phi(imol,5)
               d2rhorb(lmol,6)=phi(imol)*d2phi(jmol,6)+
     .                      2.0d0*dphi(imol,3)*dphi(jmol,3)+
     .                        phi(jmol)*d2phi(imol,6)
               d2dotorb(lmol,1)=phi(imol)*d3phi(jmol,1)-
     .                          phi(jmol)*d3phi(imol,1)+
     .                          dphi(imol,1)*d2phi(jmol,1)-
     .                          dphi(jmol,1)*d2phi(imol,1)
               d2dotorb(lmol,2)=phi(imol)*d3phi(jmol,2)-
     .                          phi(jmol)*d3phi(imol,2)+
     .                          dphi(imol,2)*d2phi(jmol,1)-
     .                          dphi(jmol,2)*d2phi(imol,1)
               d2dotorb(lmol,3)=phi(imol)*d3phi(jmol,3)-
     .                          phi(jmol)*d3phi(imol,3)+
     .                          dphi(imol,3)*d2phi(jmol,1)-
     .                          dphi(jmol,3)*d2phi(imol,1)
               d2dotorb(lmol,4)=phi(imol)*d3phi(jmol,4)-
     .                          phi(jmol)*d3phi(imol,4)+
     .                  2.0d0*(dphi(imol,2)*d2phi(jmol,2)-
     .                         dphi(jmol,2)*d2phi(imol,2))-
     .                         dphi(imol,1)*d2phi(jmol,4)+
     .                         dphi(jmol,1)*d2phi(imol,4)
               d2dotorb(lmol,5)=phi(imol)*d3phi(jmol,5)-
     .                          phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,1)*d2phi(jmol,5)+
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,2)*d2phi(jmol,3)-
     .                          dphi(jmol,2)*d2phi(imol,3)+
     .                          dphi(imol,3)*d2phi(jmol,2)-
     .                          dphi(jmol,3)*d2phi(imol,2)
               d2dotorb(lmol,6)=phi(imol)*d3phi(jmol,6)-
     .                          phi(jmol)*d3phi(imol,6)+
     .                  2.0d0*(dphi(imol,3)*d2phi(jmol,3)-
     .                         dphi(jmol,3)*d2phi(imol,3))-
     .                         dphi(imol,1)*d2phi(jmol,6)+
     .                         dphi(jmol,1)*d2phi(imol,6)
               d2dotorb(lmol,7)=phi(imol)*d3phi(jmol,2)-
     .                          phi(jmol)*d3phi(imol,2)+
     .                  2.0d0*(dphi(imol,1)*d2phi(jmol,2)-
     .                         dphi(jmol,1)*d2phi(imol,2))-
     .                         dphi(imol,2)*d2phi(jmol,1)+
     .                         dphi(jmol,2)*d2phi(imol,1)
               d2dotorb(lmol,8)=phi(imol)*d3phi(jmol,4)-
     .                          phi(jmol)*d3phi(imol,4)+
     .                          dphi(imol,1)*d2phi(jmol,4)-
     .                          dphi(jmol,1)*d2phi(imol,4)
               d2dotorb(lmol,9)=phi(imol)*d3phi(jmol,5)-
     .                          phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,2)*d2phi(jmol,3)+
     .                          dphi(jmol,2)*d2phi(imol,3)+
     .                          dphi(imol,1)*d2phi(jmol,5)-
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,3)*d2phi(jmol,2)-
     .                          dphi(jmol,3)*d2phi(imol,2)
               d2dotorb(lmol,10)=phi(imol)*d3phi(jmol,7)-
     .                          phi(jmol)*d3phi(imol,7)+
     .                          dphi(imol,2)*d2phi(jmol,4)-
     .                          dphi(jmol,2)*d2phi(imol,4)
               d2dotorb(lmol,11)=phi(imol)*d3phi(jmol,8)-
     .                          phi(jmol)*d3phi(imol,8)+
     .                          dphi(imol,3)*d2phi(jmol,4)-
     .                          dphi(jmol,3)*d2phi(imol,4)
               d2dotorb(lmol,12)=phi(imol)*d3phi(jmol,9)-
     .                          phi(jmol)*d3phi(imol,9)+
     .                  2.0d0*(dphi(imol,3)*d2phi(jmol,5)-
     .                         dphi(jmol,3)*d2phi(imol,5))-
     .                         dphi(imol,2)*d2phi(jmol,6)+
     .                         dphi(jmol,2)*d2phi(imol,6)
               d2dotorb(lmol,13)=phi(imol)*d3phi(jmol,3)-
     .                          phi(jmol)*d3phi(imol,3)+
     .                  2.0d0*(dphi(imol,1)*d2phi(jmol,3)-
     .                         dphi(jmol,1)*d2phi(imol,3))-
     .                         dphi(imol,3)*d2phi(jmol,1)+
     .                         dphi(jmol,3)*d2phi(imol,1)
               d2dotorb(lmol,15)=phi(imol)*d3phi(jmol,6)-
     .                          phi(jmol)*d3phi(imol,6)+
     .                          dphi(imol,1)*d2phi(jmol,6)-
     .                          dphi(jmol,1)*d2phi(imol,6)
               d2dotorb(lmol,14)=phi(imol)*d3phi(jmol,5)-
     .                           phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,3)*d2phi(jmol,2)+
     .                          dphi(jmol,3)*d2phi(imol,2)+
     .                          dphi(imol,1)*d2phi(jmol,5)-
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,2)*d2phi(jmol,3)-
     .                          dphi(jmol,2)*d2phi(imol,3)
               d2dotorb(lmol,16)=phi(imol)*d3phi(jmol,8)-
     .                          phi(jmol)*d3phi(imol,8)+
     .                  2.0d0*(dphi(imol,2)*d2phi(jmol,5)-
     .                         dphi(jmol,2)*d2phi(imol,5))-
     .                         dphi(imol,3)*d2phi(jmol,4)+
     .                         dphi(jmol,3)*d2phi(imol,4)
               d2dotorb(lmol,17)=phi(imol)*d3phi(jmol,9)-
     .                          phi(jmol)*d3phi(imol,9)+
     .                          dphi(imol,2)*d2phi(jmol,6)-
     .                          dphi(jmol,2)*d2phi(imol,6)
               d2dotorb(lmol,18)=phi(imol)*d3phi(jmol,10)-
     .                          phi(jmol)*d3phi(imol,10)+
     .                          dphi(imol,3)*d2phi(jmol,6)-
     .                          dphi(jmol,3)*d2phi(imol,6)
               if(imol.ne.jmol) then
                   rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                   drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                   drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                   drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
                   d2rhorb(lmol,1)=d2rhorb(lmol,1)+d2rhorb(lmol,1)
                   d2rhorb(lmol,2)=d2rhorb(lmol,2)+d2rhorb(lmol,2)
                   d2rhorb(lmol,3)=d2rhorb(lmol,3)+d2rhorb(lmol,3)
                   d2rhorb(lmol,4)=d2rhorb(lmol,4)+d2rhorb(lmol,4)
                   d2rhorb(lmol,5)=d2rhorb(lmol,5)+d2rhorb(lmol,5)
                   d2rhorb(lmol,6)=d2rhorb(lmol,6)+d2rhorb(lmol,6)
               endif
            enddo
         enddo
         d2rha(1)=0.0d0
         d2rha(2)=0.0d0
         d2rha(3)=0.0d0
         d2rha(4)=0.0d0
         d2rha(5)=0.0d0
         d2rha(6)=0.0d0
         do l=1,na
            lmol=ia(l)
            d2rha(1)=d2rha(1)+pa(l)*d2rhorb(lmol,1)
            d2rha(2)=d2rha(2)+pa(l)*d2rhorb(lmol,2)
            d2rha(3)=d2rha(3)+pa(l)*d2rhorb(lmol,3)
            d2rha(4)=d2rha(4)+pa(l)*d2rhorb(lmol,4)
            d2rha(5)=d2rha(5)+pa(l)*d2rhorb(lmol,5)
            d2rha(6)=d2rha(6)+pa(l)*d2rhorb(lmol,6)
         enddo
         d2gradfa(1)=0.0d0
         d2gradfa(2)=0.0d0
         d2gradfa(3)=0.0d0
         d2gradfa(4)=0.0d0
         d2gradfa(5)=0.0d0
         d2gradfa(6)=0.0d0
         do l=1,naa
            imol=iaa(l)
            jmol=jaa(l)
            kmol=kaa(l)
            lmol=laa(l)
            ij=indij(imol)+jmol
            kl=indij(kmol)+lmol
            d2gradfa(1)=d2gradfa(1)+paa(l)*
     .      (d2dotorb(ij,1)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,1)+
     .       ddotorb(ij,1)*ddotorb(kl,1)+ddotorb(ij,1)*ddotorb(kl,1)+
     .       d2dotorb(ij,7)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,7)+
     .       ddotorb(ij,4)*ddotorb(kl,4)+ddotorb(ij,4)*ddotorb(kl,4)+
     .       d2dotorb(ij,13)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,13)+
     .       ddotorb(ij,7)*ddotorb(kl,7)+ddotorb(ij,7)*ddotorb(kl,7))
            d2gradfa(2)=d2gradfa(2)+paa(l)*
     .      (d2dotorb(ij,2)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,2)+
     .       ddotorb(ij,1)*ddotorb(kl,2)+ddotorb(ij,2)*ddotorb(kl,1)+
     .       d2dotorb(ij,8)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,8)+
     .       ddotorb(ij,4)*ddotorb(kl,5)+ddotorb(ij,5)*ddotorb(kl,4)+
     .       d2dotorb(ij,14)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,14)+
     .       ddotorb(ij,7)*ddotorb(kl,8)+ddotorb(ij,8)*ddotorb(kl,7))
            d2gradfa(3)=d2gradfa(3)+paa(l)*
     .      (d2dotorb(ij,3)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,3)+
     .       ddotorb(ij,1)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,1)+
     .       d2dotorb(ij,9)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,9)+
     .       ddotorb(ij,4)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,4)+
     .       d2dotorb(ij,15)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,15)+
     .       ddotorb(ij,7)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,7))
            d2gradfa(4)=d2gradfa(4)+paa(l)*
     .      (d2dotorb(ij,4)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,4)+
     .       ddotorb(ij,2)*ddotorb(kl,2)+ddotorb(ij,2)*ddotorb(kl,2)+
     .       d2dotorb(ij,10)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,10)+
     .       ddotorb(ij,5)*ddotorb(kl,5)+ddotorb(ij,5)*ddotorb(kl,5)+
     .       d2dotorb(ij,16)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,16)+
     .       ddotorb(ij,8)*ddotorb(kl,8)+ddotorb(ij,8)*ddotorb(kl,8))
            d2gradfa(5)=d2gradfa(5)+paa(l)*
     .      (d2dotorb(ij,5)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,5)+
     .       ddotorb(ij,2)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,2)+
     .       d2dotorb(ij,11)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,11)+
     .       ddotorb(ij,5)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,5)+
     .       d2dotorb(ij,17)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,17)+
     .       ddotorb(ij,8)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,8))
            d2gradfa(6)=d2gradfa(6)+paa(l)*
     .      (d2dotorb(ij,6)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,6)+
     .       ddotorb(ij,3)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,3)+
     .       d2dotorb(ij,12)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,12)+
     .       ddotorb(ij,6)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,6)+
     .       d2dotorb(ij,18)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,18)+
     .       ddotorb(ij,9)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,9))
         enddo
         if(nb.eq.0) then
            d2rhb(1)=d2rha(1)
            d2rhb(2)=d2rha(2)
            d2rhb(3)=d2rha(3)
            d2rhb(4)=d2rha(4)
            d2rhb(5)=d2rha(5)
            d2rhb(6)=d2rha(6)
            d2gradfb(1)=d2gradfa(1)
            d2gradfb(2)=d2gradfa(2)
            d2gradfb(3)=d2gradfa(3)
            d2gradfb(4)=d2gradfa(4)
            d2gradfb(5)=d2gradfa(5)
            d2gradfb(6)=d2gradfa(6)
         else
**    compute beta contribution (open shell case)
         do imol=nmola+1,nmol
            do jmol=imol,nmol
              lmol=indij(imol-nmola)+jmol-nmola
               rhorb(lmol)=phi(imol)*phi(jmol)
               dotorb(lmol,1)=phi(imol)*dphi(jmol,1)-
     .                        phi(jmol)*dphi(imol,1)
               dotorb(lmol,2)=phi(imol)*dphi(jmol,2)-
     .                        phi(jmol)*dphi(imol,2)
               dotorb(lmol,3)=phi(imol)*dphi(jmol,3)-
     .                        phi(jmol)*dphi(imol,3)
               drhorb(lmol,1)=phi(imol)*dphi(jmol,1)+
     .                        phi(jmol)*dphi(imol,1)
               drhorb(lmol,2)=phi(imol)*dphi(jmol,2)+
     .                        phi(jmol)*dphi(imol,2)
               drhorb(lmol,3)=phi(imol)*dphi(jmol,3)+
     .                        phi(jmol)*dphi(imol,3)
               ddotorb(lmol,1)=phi(imol)*d2phi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,1)
               ddotorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,2)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,1)*dphi(jmol,2)
               ddotorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,3)*dphi(jmol,1)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,1)*dphi(jmol,3)
               ddotorb(lmol,4)=phi(imol)*d2phi(jmol,2)+
     .                         dphi(imol,1)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,2)-
     .                         dphi(imol,2)*dphi(jmol,1)
               ddotorb(lmol,5)=phi(imol)*d2phi(jmol,4)- 
     .                         phi(jmol)*d2phi(imol,4)
               ddotorb(lmol,6)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,3)*dphi(jmol,2)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,2)*dphi(jmol,3)
               ddotorb(lmol,7)=phi(imol)*d2phi(jmol,3)+
     .                         dphi(imol,1)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,3)-
     .                         dphi(imol,3)*dphi(jmol,1)
               ddotorb(lmol,8)=phi(imol)*d2phi(jmol,5)+
     .                         dphi(imol,2)*dphi(jmol,3)-
     .                         phi(jmol)*d2phi(imol,5)-
     .                         dphi(imol,3)*dphi(jmol,2)
               ddotorb(lmol,9)=phi(imol)*d2phi(jmol,6)- 
     .                         phi(jmol)*d2phi(imol,6)
               d2rhorb(lmol,1)=phi(imol)*d2phi(jmol,1)+
     .                      2.0d0*dphi(imol,1)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,1)
               d2rhorb(lmol,2)=phi(imol)*d2phi(jmol,2)+
     .                        dphi(imol,1)*dphi(jmol,2)+
     .                        dphi(imol,2)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,2)
               d2rhorb(lmol,3)=phi(imol)*d2phi(jmol,3)+
     .                        dphi(imol,1)*dphi(jmol,3)+
     .                        dphi(imol,3)*dphi(jmol,1)+
     .                        phi(jmol)*d2phi(imol,3)
               d2rhorb(lmol,4)=phi(imol)*d2phi(jmol,4)+
     .                      2.0d0*dphi(imol,2)*dphi(jmol,2)+
     .                        phi(jmol)*d2phi(imol,4)
               d2rhorb(lmol,5)=phi(imol)*d2phi(jmol,5)+
     .                        dphi(imol,2)*dphi(jmol,3)+
     .                        dphi(imol,3)*dphi(jmol,2)+
     .                        phi(jmol)*d2phi(imol,5)
               d2rhorb(lmol,6)=phi(imol)*d2phi(jmol,6)+
     .                      2.0d0*dphi(imol,3)*dphi(jmol,3)+
     .                        phi(jmol)*d2phi(imol,6)
               d2dotorb(lmol,1)=phi(imol)*d3phi(jmol,1)-
     .                          phi(jmol)*d3phi(imol,1)+
     .                          dphi(imol,1)*d2phi(jmol,1)-
     .                          dphi(jmol,1)*d2phi(imol,1)
               d2dotorb(lmol,2)=phi(imol)*d3phi(jmol,2)-
     .                          phi(jmol)*d3phi(imol,2)+
     .                          dphi(imol,2)*d2phi(jmol,1)-
     .                          dphi(jmol,2)*d2phi(imol,1)
               d2dotorb(lmol,3)=phi(imol)*d3phi(jmol,3)-
     .                          phi(jmol)*d3phi(imol,3)+
     .                          dphi(imol,3)*d2phi(jmol,1)-
     .                          dphi(jmol,3)*d2phi(imol,1)
               d2dotorb(lmol,4)=phi(imol)*d3phi(jmol,4)-
     .                          phi(jmol)*d3phi(imol,4)+
     .                  2.0d0*(dphi(imol,2)*d2phi(jmol,2)-
     .                         dphi(jmol,2)*d2phi(imol,2))-
     .                         dphi(imol,1)*d2phi(jmol,4)+
     .                         dphi(jmol,1)*d2phi(imol,4)
               d2dotorb(lmol,5)=phi(imol)*d3phi(jmol,5)-
     .                          phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,1)*d2phi(jmol,5)+
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,2)*d2phi(jmol,3)-
     .                          dphi(jmol,2)*d2phi(imol,3)+
     .                          dphi(imol,3)*d2phi(jmol,2)-
     .                          dphi(jmol,3)*d2phi(imol,2)
               d2dotorb(lmol,6)=phi(imol)*d3phi(jmol,6)-
     .                          phi(jmol)*d3phi(imol,6)+
     .                  2.0d0*(dphi(imol,3)*d2phi(jmol,3)-
     .                         dphi(jmol,3)*d2phi(imol,3))-
     .                         dphi(imol,1)*d2phi(jmol,6)+
     .                         dphi(jmol,1)*d2phi(imol,6)
               d2dotorb(lmol,7)=phi(imol)*d3phi(jmol,2)-
     .                          phi(jmol)*d3phi(imol,2)+
     .                  2.0d0*(dphi(imol,1)*d2phi(jmol,2)-
     .                         dphi(jmol,1)*d2phi(imol,2))-
     .                         dphi(imol,2)*d2phi(jmol,1)+
     .                         dphi(jmol,2)*d2phi(imol,1)
               d2dotorb(lmol,8)=phi(imol)*d3phi(jmol,4)-
     .                          phi(jmol)*d3phi(imol,4)+
     .                          dphi(imol,1)*d2phi(jmol,4)-
     .                          dphi(jmol,1)*d2phi(imol,4)
               d2dotorb(lmol,9)=phi(imol)*d3phi(jmol,5)-
     .                          phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,2)*d2phi(jmol,3)+
     .                          dphi(jmol,2)*d2phi(imol,3)+
     .                          dphi(imol,1)*d2phi(jmol,5)-
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,3)*d2phi(jmol,2)-
     .                          dphi(jmol,3)*d2phi(imol,2)
               d2dotorb(lmol,10)=phi(imol)*d3phi(jmol,7)-
     .                          phi(jmol)*d3phi(imol,7)+
     .                          dphi(imol,2)*d2phi(jmol,4)-
     .                          dphi(jmol,2)*d2phi(imol,4)
               d2dotorb(lmol,11)=phi(imol)*d3phi(jmol,8)-
     .                          phi(jmol)*d3phi(imol,8)+
     .                          dphi(imol,3)*d2phi(jmol,4)-
     .                          dphi(jmol,3)*d2phi(imol,4)
               d2dotorb(lmol,12)=phi(imol)*d3phi(jmol,9)-
     .                          phi(jmol)*d3phi(imol,9)+
     .                  2.0d0*(dphi(imol,3)*d2phi(jmol,5)-
     .                         dphi(jmol,3)*d2phi(imol,5))-
     .                         dphi(imol,2)*d2phi(jmol,6)+
     .                         dphi(jmol,2)*d2phi(imol,6)
               d2dotorb(lmol,13)=phi(imol)*d3phi(jmol,3)-
     .                          phi(jmol)*d3phi(imol,3)+
     .                  2.0d0*(dphi(imol,1)*d2phi(jmol,3)-
     .                         dphi(jmol,1)*d2phi(imol,3))-
     .                         dphi(imol,3)*d2phi(jmol,1)+
     .                         dphi(jmol,3)*d2phi(imol,1)
               d2dotorb(lmol,15)=phi(imol)*d3phi(jmol,6)-
     .                          phi(jmol)*d3phi(imol,6)+
     .                          dphi(imol,1)*d2phi(jmol,6)-
     .                          dphi(jmol,1)*d2phi(imol,6)
               d2dotorb(lmol,14)=phi(imol)*d3phi(jmol,5)-
     .                           phi(jmol)*d3phi(imol,5)-
     .                          dphi(imol,3)*d2phi(jmol,2)+
     .                          dphi(jmol,3)*d2phi(imol,2)+
     .                          dphi(imol,1)*d2phi(jmol,5)-
     .                          dphi(jmol,1)*d2phi(imol,5)+
     .                          dphi(imol,2)*d2phi(jmol,3)-
     .                          dphi(jmol,2)*d2phi(imol,3)
               d2dotorb(lmol,16)=phi(imol)*d3phi(jmol,8)-
     .                          phi(jmol)*d3phi(imol,8)+
     .                  2.0d0*(dphi(imol,2)*d2phi(jmol,5)-
     .                         dphi(jmol,2)*d2phi(imol,5))-
     .                         dphi(imol,3)*d2phi(jmol,4)+
     .                         dphi(jmol,3)*d2phi(imol,4)
               d2dotorb(lmol,17)=phi(imol)*d3phi(jmol,9)-
     .                          phi(jmol)*d3phi(imol,9)+
     .                          dphi(imol,2)*d2phi(jmol,6)-
     .                          dphi(jmol,2)*d2phi(imol,6)
               d2dotorb(lmol,18)=phi(imol)*d3phi(jmol,10)-
     .                          phi(jmol)*d3phi(imol,10)+
     .                          dphi(imol,3)*d2phi(jmol,6)-
     .                          dphi(jmol,3)*d2phi(imol,6)
               if(imol.ne.jmol) then
                   rhorb(lmol)=rhorb(lmol)+rhorb(lmol)
                   drhorb(lmol,1)=drhorb(lmol,1)+drhorb(lmol,1)
                   drhorb(lmol,2)=drhorb(lmol,2)+drhorb(lmol,2)
                   drhorb(lmol,3)=drhorb(lmol,3)+drhorb(lmol,3)
                   d2rhorb(lmol,1)=d2rhorb(lmol,1)+d2rhorb(lmol,1)
                   d2rhorb(lmol,2)=d2rhorb(lmol,2)+d2rhorb(lmol,2)
                   d2rhorb(lmol,3)=d2rhorb(lmol,3)+d2rhorb(lmol,3)
                   d2rhorb(lmol,4)=d2rhorb(lmol,4)+d2rhorb(lmol,4)
                   d2rhorb(lmol,5)=d2rhorb(lmol,5)+d2rhorb(lmol,5)
                   d2rhorb(lmol,6)=d2rhorb(lmol,6)+d2rhorb(lmol,6)
              endif
            enddo
         enddo
         d2rhb(1)=0.0d0
         d2rhb(2)=0.0d0
         d2rhb(3)=0.0d0
         d2rhb(4)=0.0d0
         d2rhb(5)=0.0d0
         d2rhb(6)=0.0d0
         do l=1,nb
            lmol=ib(l)
            d2rhb(1)=d2rhb(1)+pb(l)*d2rhorb(lmol,1)
            d2rhb(2)=d2rhb(2)+pb(l)*d2rhorb(lmol,2)
            d2rhb(3)=d2rhb(3)+pb(l)*d2rhorb(lmol,3)
            d2rhb(4)=d2rhb(4)+pb(l)*d2rhorb(lmol,4)
            d2rhb(5)=d2rhb(5)+pb(l)*d2rhorb(lmol,5)
            d2rhb(6)=d2rhb(6)+pb(l)*d2rhorb(lmol,6)
         enddo
         d2gradfb(1)=0.0d0
         d2gradfb(2)=0.0d0
         d2gradfb(3)=0.0d0
         d2gradfb(4)=0.0d0
         d2gradfb(5)=0.0d0
         d2gradfb(6)=0.0d0
         do l=1,nbb
            imol=ibb(l)
            jmol=jbb(l)
            kmol=kbb(l)
            lmol=lbb(l)
            ij=indij(imol)+jmol
            kl=indij(kmol)+lmol
            d2gradfb(1)=d2gradfb(1)+pbb(l)*
     .      (d2dotorb(ij,1)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,1)+
     .       ddotorb(ij,1)*ddotorb(kl,1)+ddotorb(ij,1)*ddotorb(kl,1)+
     .       d2dotorb(ij,7)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,7)+
     .       ddotorb(ij,4)*ddotorb(kl,4)+ddotorb(ij,4)*ddotorb(kl,4)+
     .       d2dotorb(ij,13)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,13)+
     .       ddotorb(ij,7)*ddotorb(kl,7)+ddotorb(ij,7)*ddotorb(kl,7))
            d2gradfb(2)=d2gradfb(2)+pbb(l)*
     .      (d2dotorb(ij,2)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,2)+
     .       ddotorb(ij,1)*ddotorb(kl,2)+ddotorb(ij,2)*ddotorb(kl,1)+
     .       d2dotorb(ij,8)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,8)+
     .       ddotorb(ij,4)*ddotorb(kl,5)+ddotorb(ij,5)*ddotorb(kl,4)+
     .       d2dotorb(ij,14)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,14)+
     .       ddotorb(ij,7)*ddotorb(kl,8)+ddotorb(ij,8)*ddotorb(kl,7))
            d2gradfb(3)=d2gradfb(3)+pbb(l)*
     .      (d2dotorb(ij,3)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,3)+
     .       ddotorb(ij,1)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,1)+
     .       d2dotorb(ij,9)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,9)+
     .       ddotorb(ij,4)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,4)+
     .       d2dotorb(ij,15)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,15)+
     .       ddotorb(ij,7)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,7))
            d2gradfb(4)=d2gradfb(4)+pbb(l)*
     .      (d2dotorb(ij,4)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,4)+
     .       ddotorb(ij,2)*ddotorb(kl,2)+ddotorb(ij,2)*ddotorb(kl,2)+
     .       d2dotorb(ij,10)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,10)+
     .       ddotorb(ij,5)*ddotorb(kl,5)+ddotorb(ij,5)*ddotorb(kl,5)+
     .       d2dotorb(ij,16)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,16)+
     .       ddotorb(ij,8)*ddotorb(kl,8)+ddotorb(ij,8)*ddotorb(kl,8))
            d2gradfb(5)=d2gradfb(5)+pbb(l)*
     .      (d2dotorb(ij,5)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,5)+
     .       ddotorb(ij,2)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,2)+
     .       d2dotorb(ij,11)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,11)+
     .       ddotorb(ij,5)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,5)+
     .       d2dotorb(ij,17)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,17)+
     .       ddotorb(ij,8)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,8))
            d2gradfb(6)=d2gradfb(6)+pbb(l)*
     .      (d2dotorb(ij,6)*dotorb(kl,1)+dotorb(ij,1)*d2dotorb(kl,6)+
     .       ddotorb(ij,3)*ddotorb(kl,3)+ddotorb(ij,3)*ddotorb(kl,3)+
     .       d2dotorb(ij,12)*dotorb(kl,2)+dotorb(ij,2)*d2dotorb(kl,12)+
     .       ddotorb(ij,6)*ddotorb(kl,6)+ddotorb(ij,6)*ddotorb(kl,6)+
     .       d2dotorb(ij,18)*dotorb(kl,3)+dotorb(ij,3)*d2dotorb(kl,18)+
     .       ddotorb(ij,9)*ddotorb(kl,9)+ddotorb(ij,9)*ddotorb(kl,9))
         enddo
         endif
         d2rho(1)=d2rha(1)+d2rhb(1)
         d2rho(2)=d2rha(2)+d2rhb(2)
         d2rho(3)=d2rha(3)+d2rhb(3)
         d2rho(4)=d2rha(4)+d2rhb(4)
         d2rho(5)=d2rha(5)+d2rhb(5)
         d2rho(6)=d2rha(6)+d2rhb(6)
         d2gradf(1)=d2gradfa(1)+d2gradfb(1)+tol*d2rho(1)
         d2gradf(2)=d2gradfa(2)+d2gradfb(2)+tol*d2rho(2)
         d2gradf(3)=d2gradfa(3)+d2gradfb(3)+tol*d2rho(3)
         d2gradf(4)=d2gradfa(4)+d2gradfb(4)+tol*d2rho(4)
         d2gradf(5)=d2gradfa(5)+d2gradfb(5)+tol*d2rho(5)
         d2gradf(6)=d2gradfa(6)+d2gradfb(6)+tol*d2rho(6)
         d2ratio(1)=dratio(1)*(2.0d0*dgradf(1)/gradf+p163*drho(1)/rho)+
     .       2.0d0*ratio*(d2gradf(1)-dgradf(1)*dgradf(1)/gradf)/gradf+
     .       p163*ratio*(d2rho(1)-drho(1)*drho(1)/rho)/rho
         d2ratio(2)=dratio(2)*(2.0d0*dgradf(1)/gradf+p163*drho(1)/rho)+
     .       2.0d0*ratio*(d2gradf(2)-dgradf(1)*dgradf(2)/gradf)/gradf+
     .       p163*ratio*(d2rho(2)-drho(1)*drho(2)/rho)/rho
         d2ratio(3)=dratio(3)*(2.0d0*dgradf(1)/gradf+p163*drho(1)/rho)+
     .       2.0d0*ratio*(d2gradf(3)-dgradf(1)*dgradf(3)/gradf)/gradf+
     .       p163*ratio*(d2rho(3)-drho(1)*drho(3)/rho)/rho
         d2ratio(4)=dratio(2)*(2.0d0*dgradf(2)/gradf+p163*drho(2)/rho)+
     .       2.0d0*ratio*(d2gradf(4)-dgradf(2)*dgradf(2)/gradf)/gradf+
     .       p163*ratio*(d2rho(4)-drho(2)*drho(2)/rho)/rho
         d2ratio(5)=dratio(3)*(2.0d0*dgradf(2)/gradf+p163*drho(2)/rho)+
     .       2.0d0*ratio*(d2gradf(5)-dgradf(2)*dgradf(3)/gradf)/gradf+
     .       p163*ratio*(d2rho(5)-drho(2)*drho(3)/rho)/rho
         d2ratio(6)=dratio(3)*(2.0d0*dgradf(3)/gradf+p163*drho(3)/rho)+
     .       2.0d0*ratio*(d2gradf(6)-dgradf(3)*dgradf(3)/gradf)/gradf+
     .       p163*ratio*(d2rho(6)-drho(3)*drho(3)/rho)/rho
      d2elf(1)=elfd*elfd*(2.0d0*elfd*dratio(1)*dratio(1)-d2ratio(1))
      d2elf(2)=elfd*elfd*(2.0d0*elfd*dratio(1)*dratio(2)-d2ratio(2))
      d2elf(3)=elfd*elfd*(2.0d0*elfd*dratio(1)*dratio(3)-d2ratio(3))
      d2elf(4)=elfd*elfd*(2.0d0*elfd*dratio(2)*dratio(2)-d2ratio(4))
      d2elf(5)=elfd*elfd*(2.0d0*elfd*dratio(2)*dratio(3)-d2ratio(5))
      d2elf(6)=elfd*elfd*(2.0d0*elfd*dratio(3)*dratio(3)-d2ratio(6))
      endif
***** TMP2
      return
      end
