      module rohf
        ! Restricted open-shell Hartree-Fock.
        use mtrx
        use uhf
        implicit none
        private
        integer::maxit_=20,nfock_=0
        real(8)::conv_=1d-4
        !---------------------------------------------------------------
        public::
     &      rohf_scf,
     &      rohf_eelc
        !---------------------------------------------------------------

      contains

        subroutine rohf_scf(nea,neb,sb,hb,tb,cm,em)
          implicit none
          integer,intent(in)::nea,neb
          real(8),intent(in)::sb(:),hb(:),tb(:)
          real(8),intent(inout)::cm(:,:)
          real(8),intent(out)::em(:)
          integer::i,nb,nb2
          real(8)::delta,ee
          real(8),allocatable::
     &      xb(:,:),pba(:),pbb(:),fb(:),pba0(:),fba(:),fbb(:)
          nb=SIZE(em)
          nb2=nb*(nb+1)/2
          allocate(xb(nb,nb))
          allocate(pba(nb2))
          allocate(pbb(nb2))
          allocate(fb(nb2))
          allocate(pba0(nb2))
          allocate(fba(nb2))
          allocate(fbb(nb2))
          call mtrx_xcanon(sb,xb)
          call uhf_dens(nea,cm,pba)
          call uhf_dens(neb,cm,pbb)
          do i=1,maxit_
            pba0(:)=pba(:)
            call uhf_fock(pba,pbb,hb,tb,fba)
            call uhf_fock(pbb,pba,hb,tb,fbb)
            call fock_(nea,neb,cm,sb,fba,fbb,fb)
            call mtrx_utau(xb,fb)
            call mtrx_eigsp(fb,em,cm)
            call mtrx_uc(xb,cm)
            call uhf_dens(nea,cm,pba)
            call uhf_dens(neb,cm,pbb)
            delta=SQRT((SUM(pba-pba0)**2)/4)
            if (delta<conv_) then
              exit
            end if
            if (i==maxit_) then
              call rohf_eelc(nea,neb,hb,tb,cm,ee)
              write(*,'("SCF did not converge")')
              write(*,'("i,ee,delta=",i4,2f13.6)') i,ee,delta
              stop(1)
            end if
          end do
        end subroutine

        subroutine rohf_eelc(nea,neb,hb,tb,cm,val)
          implicit none
          integer,intent(in)::nea,neb
          real(8),intent(in)::hb(:),tb(:),cm(:,:)
          real(8),intent(out)::val
          call uhf_eelc(nea,neb,hb,tb,cm,cm,val)
        end subroutine

        subroutine fockm_(nea,neb,cm,fba,fbb,val)
          ! Combined Fock matrix in the MO basis.
          implicit none
          integer,intent(in)::nea,neb
          real(8),intent(in)::cm(:,:),fba(:),fbb(:)
          real(8),intent(out)::val(:)
          integer::i,j,ij,nb,nb2
          real(8)::acc,bcc,aoo,boo,avv,bvv
          real(8),allocatable::fa(:),fb(:)
          call fockparam_(acc,bcc,aoo,boo,avv,bvv)
          nb2=SIZE(fba)
          call mtrx_iup2(nb2,i,nb)
          allocate(fa(nb2))
          allocate(fb(nb2))
          ! Fock matrix in the MO basis.
          fa(:)=fba(:)
          fb(:)=fbb(:)
          call mtrx_utau(cm,fa)
          call mtrx_utau(cm,fb)
          ! Form combined Fock matrix in the MO basis.
          val(:)=0d0
          do i=1,nb
            do j=1,i
              call mtrx_ip2(i,j,ij)
              if (i<=neb.and.j<=neb) then
                val(ij)=acc*fa(ij)+bcc*fb(ij)
              else if (neb<i.and.i<=nea.and.neb<j.and.j<=nea) then
                val(ij)=aoo*fa(ij)+boo*fb(ij)
              else if (i>nea.and.j>nea) then
                val(ij)=avv*fa(ij)+bvv*fb(ij)
              else if (neb<i.and.i<=nea.and.j<=neb) then
                val(ij)=fb(ij)
              else if (nea<i.and.j<=neb) then
                val(ij)=(fa(ij)+fb(ij))/2
              else if (nea<i.and.neb<j.and.j<=nea) then
                val(ij)=fa(ij)
              else
                write(*,'("error: out of range: fock()")')
              end if
            end do
          end do
        end subroutine

        subroutine fock_(nea,neb,cm,sb,fba,fbb,val)
          ! Combined Fock matrix in the AO basis.
          implicit none
          integer,intent(in)::nea,neb
          real(8),intent(in)::cm(:,:),sb(:),fba(:),fbb(:)
          real(8),intent(out)::val(:)
          integer::i,j,ij,nb,nb2
          real(8),allocatable::fa(:),fb(:),ct(:,:)
          nb2=SIZE(fba)
          call mtrx_iup2(nb2,i,nb)
          allocate(fa(nb2))
          allocate(fb(nb2))
          allocate(ct(nb,nb))
          ! The inverse of MO coefficients.
          do ij=1,SIZE(sb)
            call mtrx_iup2(ij,i,j)
            ct(i,j)=sb(ij)
            ct(j,i)=sb(ij)
          end do
          call mtrx_uc(TRANSPOSE(cm),ct)
          ! Fock matrix in the MO basis.
          fa(:)=fba(:)
          fb(:)=fbb(:)
          call mtrx_utau(cm,fa)
          call mtrx_utau(cm,fb)
          ! Form combined Fock matrix in the MO basis.
          call fockm_(nea,neb,cm,fba,fbb,val)
          ! Transform the Fock matrix to the AO basis.
          call mtrx_utau(ct,val)
        end subroutine

        subroutine fockparam_(acc,bcc,aoo,boo,avv,bvv)
          ! Canonicalization coefficients.
          implicit none
          real(8),intent(out)::acc,bcc,aoo,boo,avv,bvv
          ! Roothaan single matrix
          if (nfock_==0) then
            acc=-0.5d0
            bcc=1.5d0
            aoo=0.5d0
            boo=0.5d0
            avv=1.5d0
            bvv=-0.5d0
          else
            write(*,'("error: no such param: fockparam_()")')
            stop(1)
          end if
        end subroutine

      end module
