C 20 Aug 07 - KK,DGF - new module added
c
C*MODULE fmopbc  *DECK cellvol
      function cellvol(abclat,anglat)
      implicit double precision(a-h,o-z)
      parameter (eps=1.0d-06,one=1.0d+00,two=2.0d+00)
      dimension abclat(3),anglat(3)
c
      ca=cos(anglat(1))
      cb=cos(anglat(2))
      cg=cos(anglat(3))
      if(abs(ca).le.eps) ca=0 
      if(abs(cb).le.eps) cb=0 
      if(abs(cg).le.eps) cg=0 
c
      cellvol=abclat(1)*abclat(2)*abclat(3)*
     *        sqrt(one-ca*ca-cb*cb-cg*cg+two*ca*cb*cg)
c
      return
      end
c
C*MODULE fmopbc  *DECK chksymop
      subroutine chksymop(nat,cc,cct1,cct2)
      implicit double precision(a-h,o-z)
      parameter (rlarge=1.0d+12,rshort=0.25d+00)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      common /fmolat/ untang(3),untorg(3),respbc(4),abclat(3),anglat(3),
     *                symtra(3,24),symope(3,3,24),iatorg,nsymop,maxklms,
     *                ioporg(3),iopdir(3),iopabc(3),iopang(3)
      dimension cc(3,*),cct1(3,*),cct2(3,*)
      dimension trmat(3,3)
c     dimension xyzbox(3,8),u(3,3)
c     data      xyzbox/zero,zero,zero, one,zero,zero,
c    *                 zero,one,zero,  zero,zero,one,
c    *                 one,one,zero,   zero,one,one,
c    *                 one,zero,one,   one,one,one  /
c     check coordinates gererated by symmetry operations
      ismin=0
      jsmin=0
      iamin=0
      jamin=0
      do ig=1,nsymop
         call symolcc(0,0,0,ig,nat,cc,cct1,trmat)
c        if(idebug.eq.1) then
c           write(iw,9999) ik,il,im,ig
c           call prtmolcc(nat,ian,cct1)
c        endif
         do jg=ig,nsymop
            if(ig.ne.jg) then
               call symolcc(0,0,0,jg,nat,cc,cct2,trmat)
               rmin2=rlarge
               do i=1,nat
                  do j=1,nat
                     rij2=(cct1(1,i)-cct2(1,j))**2
     *                   +(cct1(2,i)-cct2(2,j))**2
     *                   +(cct1(3,i)-cct2(3,j))**2
                     if(rij2.lt.rmin2) then
                        rmin2=rij2
                        ismin=ig
                        jsmin=jg
                        iamin=i
                        jamin=j
                     endif
                  enddo
               enddo
            endif
         enddo
      enddo
c
c     if(idebug.ne.0) then
c        rij2=dsqrt(rmin2)/tobohr
c        write(iw,9998) rij2,iamin,jamin,ismin,jsmin
c        write(iw,9997)
c        call txctox(u,ierr)
c        do i=1,8
c           xf=xyzbox(1,i)
c           yf=xyzbox(2,i)
c           zf=xyzbox(3,i)
c           tx=u(1,1)*xf+u(1,2)*yf+u(1,3)*zf
c           ty=u(2,1)*xf+u(2,2)*yf+u(2,3)*zf
c           tz=u(3,1)*xf+u(3,2)*yf+u(3,3)*zf
c           write(iw,9996) i,tx/tobohr,ty/tobohr,tz/tobohr
c        enddo
c     endif
      if(rmin2.lt.rshort) then
         write(iw,9995) dsqrt(rmin2),iamin,jamin,ismin,jsmin
         call abrt
      endif
c
      return
c9999 format(' chksymop: atomic coordinates for (k,l,m,is)',4i4)
c9998 format(' chksymop: the shortest interatomic distance is, r=',f12.6
c    *,' A',/,10x,' between ia =',i3,' and ja =',i3,' generated by symme
c    *try operations,',/,10x,' is =',i2,' and js =',i2,
c    * ' , respectively, in (0,0,0) cell.')
c9997 format(' chksymop: cartesian coordinates of (0,0,0) box corners (A
c    *)')
c9996 format(i5,3f12.6)
 9995 format(' chksymop: too short interatomic diatance, r=',f12.6,
     * ' A',/,10x,' between ia =',i3,' and ja =',i3,' generated by symme
     *try operations,',/,10x,' is =',i2,' and js =',i2,
     * ' , respectively, in (0,0,0) cell.')
      end
c
C*MODULE fmopbc  *DECK fmon123
      subroutine fmon123(numfrg,n1fmo,n2fmo,n3fmo)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension numfrg(*)
      common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
      n1fmo=0
      n2fmo=0
      n3fmo=0
      do ifg=1,nfg
        l1i=iand(numfrg(ifg),65535)
        n1fmo=n1fmo+l1i
        n2fmo=n2fmo+(l1i*l1i+l1i)/2
        n3fmo=n3fmo+l1i*l1i
      enddo
      return
      end
c
C*MODULE fmopbc  *DECK genklms
      subroutine genklms(k0,l0,m0,ig0,ngrint,ngresp,ngrptc,klms)
      implicit double precision(a-h,o-z)
      parameter (zero=0.0d+00)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      common /fmolat/ untang(3),untorg(3),respbc(4),abclat(3),anglat(3),
     *                symtra(3,24),symope(3,3,24),iatorg,nsymop,maxklms,
     *                ioporg(3),iopdir(3),iopabc(3),iopang(3)
      dimension klms(4,*)
      dimension mxklm(3),miklm(3),rrmax2(4)
      dimension u(3,3),w(3,3),xyz(3),xyz0(3),xyzi(3)
c
c     set groups which are within threshould distance.
c     note that distance is measured from (xorg,yorg,zorg).
c     rcutint < rcutesp < rptcint. rptcint is possibly zero (no use).
c
c     idebug=1
c
      ngrint=0
      ngresp=0
      ngrptc=0
      rcutint=respbc(1)
      rcutesp=respbc(2)
      rptcint=respbc(3)
c     if(igroup.le.0) then
         rmax=max(rcutint,rcutesp,rptcint)
c     else
c        rmax=rcutesp
c     endif 
      rmax2=rmax*rmax
      xyz(1)=zero
      xyz(2)=zero
      xyz(3)=zero
      call symolcc(k0,l0,m0,ig0,1,xyz,xyz0,w)
c
c     rough estimation of maximum k,l,m
      mxklm(1)=int(rmax/abclat(1))
      mxklm(2)=int(rmax/abclat(2))
      mxklm(3)=int(rmax/abclat(3))
      do i=1,3
         if(mxklm(i).gt.0) mxklm(i)=mxklm(i)-1
      enddo
c     estimation of maximum k,l,m
      call txctox(u,ierr)
      do i=1,3
   20    continue
         mxklm(i)=mxklm(i)+1
         xi=mxklm(i)*u(1,i)
         yi=mxklm(i)*u(2,i)
         zi=mxklm(i)*u(3,i)
         r00=xi**2+yi**2+zi**2
         if(r00.le.rmax2) go to 20
      enddo
c
      miklm(1)=-mxklm(1)+k0
      miklm(2)=-mxklm(2)+l0
      miklm(3)=-mxklm(3)+m0
      mxklm(1)= mxklm(1)+k0
      mxklm(2)= mxklm(2)+l0
      mxklm(3)= mxklm(3)+m0
c
c     if(idebug.eq.1) then
c        write(iw,9999) (mxklm(i),i=1,3)
c        write(iw,9998) (miklm(i),i=1,3)
c     endif
c     if(igroup.le.0) then
         nset=3
         if(rptcint.eq.zero) nset=2
         rrmax2(1)=zero
         rrmax2(2)=rcutint**2
         rrmax2(3)=rcutesp**2
         rrmax2(4)=rptcint**2
c     else
c        nset=1
c        rrmax2(1)=zero
c        rrmax2(2)=rcutesp**2
c     endif
c     start of cell loop
      ncount=0
      do ii=1,nset
         do ik=miklm(1),mxklm(1)
            do il=miklm(2),mxklm(2)
               do im=miklm(3),mxklm(3)
                  do ig=1,nsymop
                     if(ik.eq.k0.and.il.eq.l0.and.im.eq.m0.and.
     *                  ig.eq.ig0) go to 100
c                    generate centers of (ik,il,im,ig) molecules
                     xyz(1)=zero
                     xyz(2)=zero
                     xyz(3)=zero
                     call symolcc(ik,il,im,ig,1,xyz,xyzi,w)
                     r0i=(xyz0(1)-xyzi(1))**2+(xyz0(2)-xyzi(2))**2
     *                                       +(xyz0(3)-xyzi(3))**2
                     if(r0i.gt.rrmax2(ii).and.r0i.le.rrmax2(ii+1)) then
                        ncount=ncount+1
                        if(ncount.gt.maxklms) call abrt
                        klms(1,ncount)=ik
                        klms(2,ncount)=il
                        klms(3,ncount)=im
                        klms(4,ncount)=ig
                     endif
  100                continue
                  enddo
               enddo
            enddo
         enddo
         if(ii.eq.1) ngrint=ncount
         if(ii.eq.2) ngresp=ncount
         if(ii.eq.3) ngrptc=ncount
      enddo
c     if(igroup.eq.1) then
c        ngrtot=ncount
c        if(ngrptc.eq.0) ngrptc=ngrtot
c     endif
c
c     if(igroup.le.0) then
         if(nset.eq.2) ngrptc=ngresp
         write(iw,9996) ngrint,ngresp,ngrptc
c     else
c        ngresp=ngrint
c        write(iw,9995) igroup,ngrint
c     endif
c     end cell loop
c
c     print lattice index list
c
      return
c9995 format(' igroup =',i4,/,' ngresp =',i10)
c9996 format(' igroup =',i4,/,' ngrint =',i10,/,' ngrint+ngresp =',
 9996 format(1x,' ngrint =',i10,/,' ngrint+ngresp =',
     *         i10,/,' ngrint+ngresp+ngrptc =',i10)
      end
c
C*MODULE fmopbc  *DECK genlatv
      subroutine genlatv(ifg,l1,l2,n2fmo,klms,untrot,espi,espg,r,wrk,
     *                   esplat)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK
      parameter (MXATM=2000,MXSH=5000)
c     COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                TT(432),INVT(48),NT
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      dimension klms(4,*),untrot(3,3,*),espi(*),espg(*),r(*),
     *          wrk(*),esplat(n2fmo,*)
c
c     Compute rotated monomer potentials VI for all group operations. 
c     Save them to ESPLAT only on masters to avoid double counting in 
c     the global sums later.
c
c     Assume equal in size fragments!
      maxl2=(maxl1*maxl1+maxl1)/2
      ind=(ifg-1)*maxl2+1
      if(maswrk) call dcopy(l2,espi,1,esplat(ind,1),1)
c     if(maswrk) write(iw,9000) 1
c     call prtri(espi,l1)
      do iu=1,nunesp
        ik=klms(1,iu)
        il=klms(2,iu)
        im=klms(3,iu)
        ig=klms(4,iu)
        if(ik.eq.0.and.il.eq.0.and.im.eq.0) then
c         if(maswrk) write(iw,9000) ig
          call dcopy(3*3,untrot(1,1,iu+1),1,tt,1)
          call TRPOSQ(tt,3)
c         Note we produce inverse rotation here, as needed below.
          call trmat
          call rotmom(r,l1,l1,0,0)
          call TFTRI(espg,espi,r,WRK,l1,l1,l1)
          if(maswrk) call dcopy(l2,espg,1,esplat(ind,ig),1)
c         call prtri(espg,l1)
        else
c         do ifg=1,nfg
c           iifg=ifg
c           call makemol(iifg,0,0,ilay,1,0,0,0,0,0,0,.true.)
c         enddo
        endif
      enddo
      call RUNITV(3,3,tt)
      call trmat
c
      return
c9000 format(1x,'Computing ESP for group operation = ',I2)
      end
c
C*MODULE fmopbc  *DECK pairvij
      subroutine pairvij(ilay,ifg,jfg,iu,nati,natfmob,nunptc,untxyz,
     *                   cfrg,espij,esp1i,esp1j,czi,czj,dolat)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical dolat 
      dimension untxyz(3,natfmob,0:nunptc+1),cfrg(*),espij(*),esp1i(*),
     *          esp1j(*),czi(*),czj(*)
      parameter (MXATM=2000,one=1.0D+00)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      common /fmolat/ untang(3),untorg(3),respbc(4),abclat(3),anglat(3),
     *                symtra(3,24),symope(3,3,24),iatorg,nsymop,maxklms,
     *                ioporg(3),iopdir(3),iopabc(3),iopang(3)
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
c
      nat1es=nat1e
      ncurshs=ncursh
      maxklmss=maxklms
      iskipesps=iskipesp
      nunesps=nunesp
c     ifmostps=ifmostp
c     modesps=modesp
      nunesp=0
      ifmostp=4
c     modesp=iand(modesp,inot(4))
      if(.not.dolat) maxklms=1
c     (to trigger fmoatfrg into following the lattice code).
      iskipesp=0
c
c     Generate JFG with its coordinates (unit cell iu)
c
      if(dolat) call dcopy(3*natfmob,untxyz(1,1,iu),1,cfrg,1)
      call CLOSDA('DELETE')
      CALL OPENDA(0)
      call makemol(jfg,0,0,ilay,0,0,0,0,0,0,0,.true.)
      natj=nat
      l2=(num*num+num)/2
c
c     Store ifg coordinates for the 2e ESP calculation 
c     (ifg is always from unit cell 0)
c
      if(dolat) call dcopy(3*natfmob,untxyz(1,1,0),1,cfrg,1)
c
c     Compute ESP VJ(I) (store in esp1j), J-block. 
c
      nat1e=natj+nati
      icurfg=jfg
      jcurfg=ifg
      ncursh=-1
c     A silly way to signal FMOESP. 
c
c     Store ifg coordinates for the 1e ESP calculation 
c
      call dcopy(nati,czi,1,zan(natj+1),1)
      call dcopy(nati*3,czi(nati+1),1,c(1,natj+1),1)
c
      call oneei
c     write(6,*) 'wwwV1j',jfg
c     call prtril(espij,num)
      call daxpy(l2,-one,espij,1,esp1j,1)
      call dcopy(natj,zan,1,czj,1)
      call dcopy(natj*3,c,1,czj(natj+1),1)
c
c     Generate IFG with its coordinates
c
      call CLOSDA('DELETE')
      CALL OPENDA(0)
      call makemol(ifg,0,0,ilay,0,0,0,0,0,0,0,.true.)
      l2=(num*num+num)/2
c
c     Store jfg coordinates for the ESP calculation
c
      if(dolat) call dcopy(3*natfmob,untxyz(1,1,iu),1,cfrg,1)
c
c     Compute ESP VI(J) (store in esp1i), I-block. 
c
c     nat1e is the same.
c
      icurfg=ifg
      jcurfg=jfg
      ncursh=-1-iu
      call dcopy(natj,czj,1,zan(nati+1),1)
      call dcopy(natj*3,czj(natj+1),1,c(1,nati+1),1)
c
      call oneei
c     write(6,*) 'wwwV1i',ifg
c     call prtril(espij,num)
      call daxpy(l2,-one,espij,1,esp1i,1)
c
      ncursh=ncurshs
      nat1e=nat1es
      nunesp=nunesps
c     modesp=modesps
c     ifmostp=ifmostps
      maxklms=maxklmss
      iskipesp=iskipesps
c
      if(dolat) call dcopy(3*natfmob,untxyz(1,1,iu),1,cfrg,1)
      return
      end
c
C*MODULE fmopbc  *DECK rotmo
      subroutine rotmo(v,ldv,nmo,l1i,ni,rotv)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (MXSH=5000,MXGTOT=20000)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      dimension v(ldv,nmo),rotv(ldv,nmo)
c
c     Rotate LCAO coefficients using precomputed AO transformation matrices
c     (set up in TRMAT, for a given rotation). In plain language, find the 
c     new MOs, after the atoms are rotated. 
c     MOs in V are a L1xL1 matrix, of which only the first NMO are rotated.
c     Rotated MOs in ROTV are a L1xL1 matrix, only the first NMO columns are set
c     (L1i is the argument to shift dimer KLOC into monomer; 
c     L1 itself is tacitly defined by NSHEL/basis set).
c     The rotation is easily set up by storing the rotation matrix in TT,
c     assigning NT=1 and calling TRMAT.
c     Ni gives the number of shells to skip in NSHEL (can be 0). 
c
c     call prsq(ptr,3,3,3)
c     call prsq(v,ldv,ldv,ldv)
      do ii=ni+1,nshell
        mini=kmin(ii)
        loci=kloc(ii)-l1i
        LIT = KTYPE(II)
c       write(6,*) 'wwwi',ii,mini,loci,LIT
c       S and the S part of L: spherically symmetric
        if(mini.eq.1) then
          call dcopy(nmo,v(loci,1),ldv,rotv(loci,1),ldv)
c         shift to P for L shells
          if(lit.eq.2) loci=loci+1
        endif
        if(lit.eq.2) then
          call MRARBR(ptr, 3, 3, 3,v(loci,1),ldv,nmo,rotv(loci,1),ldv)
        else if(lit.eq.3) then 
          call MRARBR(dtr, 6, 6, 6,v(loci,1),ldv,nmo,rotv(loci,1),ldv)
        else if(lit.eq.4) then
          call MRARBR(ftr,10,10,10,v(loci,1),ldv,nmo,rotv(loci,1),ldv)
        else if(lit.eq.5) then
          call MRARBR(gtr,15,15,15,v(loci,1),ldv,nmo,rotv(loci,1),ldv)
        endif
      enddo
c     call prsq(rotv,ldv,ldv,ldv)
      return
      END
c
C*MODULE fmopbc  *DECK rotmom
      subroutine rotmom(r,ldr,l1,ni,l1i)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (MXSH=5000,MXGTOT=20000,one=1.0D+00)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),ch(mxgtot),ci(mxgtot),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      dimension r(ldr,l1)
c
c     Prepare the L1xL1 matrix to rotate MO coefficients (LCAO), for some
c     rotation (see rotmo).
c     ni,l1i gives the number of shells/AOs to skip (typically, 0). 
c
      call vclr(r,1,l1*l1)
      do ii=ni+1,nshell
        mini=kmin(ii)
        loci=kloc(ii)-l1i
        LIT = KTYPE(II)
c       write(6,*) 'wwwi',ii,mini,loci,LIT
c       S and the S part of L: spherically symmetric
        if(mini.eq.1) then
          r(loci,loci)=one
c         shift to P for L shells
          if(lit.eq.2) loci=loci+1
        endif
        if(lit.eq.2) then
          call blkcpy(ptr,3,3,r,l1,l1,loci,loci)
        else if(lit.eq.3) then 
          call blkcpy(dtr,6,6,r,l1,l1,loci,loci)
        else if(lit.eq.4) then
          call blkcpy(ftr,10,10,r,l1,l1,loci,loci)
        else if(lit.eq.5) then
          call blkcpy(gtr,15,15,r,l1,l1,loci,loci)
        endif
      enddo
c     call prsq(r,l1,ldr,l1)
      return
      END
c
C*MODULE fmopbc  *DECK symolcc
      subroutine symolcc(ik,il,im,ig,nat,ccinp,ccout,trmat)
      implicit double precision(a-h,o-z)
c     COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      common /fmolat/ untang(3),untorg(3),respbc(4),abclat(3),anglat(3),
     *                symtra(3,24),symope(3,3,24),iatorg,nsymop,maxklms,
     *                ioporg(3),iopdir(3),iopabc(3),iopang(3)
      dimension ccinp(3,*),ccout(3,*),trmat(3,3)
      dimension u(3,3),uinv(3,3),v(3,3)
c
c     generates molecules by symmetry operations, (kk,ll,mm,ig).
c     Form U * Tg * U**-1, 
c     U unit reorientation matrix (from UNTLAT)
c     Tg group symmetry operation. 
      call txctox(u,ierr)
      call matinv3(u,uinv,ierr)
      call MRARBR(symope(1,1,ig),3,3,3,uinv,3,3,v,3)
      call MRARBR(u,3,3,3,v,3,3,trmat,3)
c     call prsq(trmat,3,3,3)
c
c     translation in cell coordinates
      xf=untorg(1)+ik+symtra(1,ig)
      yf=untorg(2)+il+symtra(2,ig)
      zf=untorg(3)+im+symtra(3,ig)
c     translation in cartesian coordinates
      tx=u(1,1)*xf+u(1,2)*yf+u(1,3)*zf
      ty=u(2,1)*xf+u(2,2)*yf+u(2,3)*zf
      tz=u(3,1)*xf+u(3,2)*yf+u(3,3)*zf
c     generate atomic coordinates
      do i=1,nat
         xf=ccinp(1,i)
         yf=ccinp(2,i)
         zf=ccinp(3,i)
         ccout(1,i)=trmat(1,1)*xf+trmat(1,2)*yf+trmat(1,3)*zf+tx
         ccout(2,i)=trmat(2,1)*xf+trmat(2,2)*yf+trmat(2,3)*zf+ty
         ccout(3,i)=trmat(3,1)*xf+trmat(3,2)*yf+trmat(3,3)*zf+tz
c        write(6,*) 'xyz',ccout(1,i)-tx,ccout(2,i)-ty,ccout(3,i)-tz
      enddo
c     call prsq(ccinp,nat,3,3)
c     write(6,*) ' '
c     call prsq(ccout,nat,3,3)
c
c     if(idebug.eq.1) then
c        write(iw,9999) ik,il,im,ig
c        call prtmolcc(nat,ian,ccout)
c     endif
      return
c9999 format(' symolcc: generated atomic cartesian coordinates (A)',/,
c    *       '          for (k,l,m,is)=(',4i3,' )')
      end
c
C*MODULE fmopbc  *DECK txctox
      subroutine txctox(u,ierr)
      implicit double precision(a-h,o-z)
      parameter (zero=0.0d+00,one=1.0d+00,eps=1.0d-16)
      common /fmolat/ untang(3),untorg(3),respbc(4),abclat(3),anglat(3),
     *                symtra(3,24),symope(3,3,24),iatorg,nsymop,maxklms,
     *                ioporg(3),iopdir(3),iopabc(3),iopang(3)
      dimension u(3,3)
c     u ... transformation matrix from fractional coordinates
c           to cartessian coordinates, x(cc)=u*x(fc).
c     put b-axis in y-axis and a-axis in (x,y)-plane, and c-axis in z-axis.
c
      ierr=0
      a=abclat(1)  
      b=abclat(2)  
      c=abclat(3)  
      cos4=cos(anglat(1))
      cos5=cos(anglat(2))
      cos6=cos(anglat(3))
      sin6=sin(anglat(3))
      if(abs(sin6).lt.eps) then
         ierr=1
         return
      endif
      v = (cos5-cos4*cos6)
      vv=one-cos4*cos4-cos5*v-cos6*(cos6-cos4*cos5)
      vv=sqrt(vv)
      v = v/sin6
      u(1,1)=a*sin6
      u(1,2)=zero
      u(1,3)=c*v
      u(2,1)=a*cos6
      u(2,2)=b
      u(2,3)=c*cos4
      u(3,1)=zero
      u(3,2)=zero
      u(3,3)=c*vv/sin6
      return
      end
c
C*MODULE fmopbc  *DECK uniqpair
      subroutine uniqpair(fmozan,untxyz,unitnr,nsymeq,natfmob,nunptc,
     *                    iu1,iu2,nsu)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK
      parameter (tol=1.0D-10,np=24)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension fmozan(*),untxyz(3,natfmob,0:nunptc),unitnr(*),
     *          nsymeq(*),lp(np)
c
c     24 - roughly corresponds to the group operations 
c
      ntot=0
      nsu=0
      do iu=iu1,iu2
        en=0
        do iat=1,natfmo
          cx=untxyz(1,iat,0)
          cy=untxyz(2,iat,0)
          cz=untxyz(3,iat,0)
          za=fmozan(iat)
          do jat=1,natfmo
            en=en+za*fmozan(jat)/sqrt((untxyz(1,jat,iu)-cx)**2+
     *               (untxyz(2,jat,iu)-cy)**2+(untxyz(3,jat,iu)-cz)**2)
          enddo
        enddo
        unitnr(iu)=en
c
        nsymeq(iu)=1
        do ju=iu1,iu-1
          if(abs(en-unitnr(ju)).lt.tol) then
            nsymeq(ju)=nsymeq(ju)+1
            nsymeq(iu)=-ju
            unitnr(iu)=0
            goto 100
          endif
        enddo
        nsu=nsu+1
  100   continue
        ntot=ntot+1
      enddo
      if(maswrk) then
        write(iw,9900)
        do iu=iu1,iu2
          if(nsymeq(iu).gt.0) then
            ii=0 
            do ju=iu+1,iu2
              if(nsymeq(ju).eq.-iu.and.ii.lt.np) then
                ii=ii+1 
                lp(ii)=ju
              endif
            enddo
            write(iw,9910) iu,nsymeq(iu),(lp(i),i=1,ii)
          endif
        enddo
        write(iw,9920) ntot,nsu
      endif
      return
 9900 format(1x,'Equivalent groups are:',/1x,'  i    Ni  equiv groups')
 9910 format(i5,i5,9999(i5))
 9920 format(1x,'Total number of groups=',I5,', unique=',I5,/)
      end
c
C*MODULE fmopbc  *DECK pbcvij
      subroutine addv2(l1,l2,mapi,mapj,vi,vj,vij)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension mapi(*),mapj(*),vi(*),vj(*),vij(*)
c
      call vclr(vij,1,l2)
      loop=0
      do i=1,l1
        ii=mapi(i)
        ji=mapj(i)
        do j=1,i
          loop=loop+1
          ij=mapi(j)
          jj=mapj(j)
c         fill the I-block, as V(I)-VI(J)
          if(ii.ne.0.and.ij.ne.0) then
            ii1=max(ii,ij)
            ij1=min(ii,ij)
            loopi=(ii1*ii1-ii1)/2+ij1 
            vij(loop)=vi(loopi)
          endif
c         fill the J-block, as V(J)-VJ(I)
          if(ji.ne.0.and.jj.ne.0) then
            ji1=max(ji,jj)
            jj1=min(ji,jj)
            loopj=(ji1*ji1-ji1)/2+jj1
c           vij(loop)=vij(loop)+vj(loopj)
            vij(loop)=vj(loopj)
c           The overlapping piece is computed twice (in I and J),
c           with exactly the same values, since it corresponds
c           to the effect of environment upon the "AO" junction of
c           two connected fragments. The values shold NOT be added. 
          endif
        enddo
      enddo
c     call prtril(vij,l1)
c     call vclr(vij,1,l2)  
      return
      end
C*MODULE fmopbc  *DECK latrij
      SUBROUTINE latrij(ifg,indat,fmozan,natfmob,untxyz,iabdfg,jabdfg,
     *                  rij,rwrk,vdwrad,nunptc,nsymeq)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      dimension indat(*),fmozan(*),untxyz(3,natfmob,*),iabdfg(*),
     *          jabdfg(*),vdwrad(*),rij(*),rwrk(*),nsymeq(*)
      common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
c
      mfg2=(nfg*nfg-nfg)/2
      noff=mfg2+(ifg-1)*nfg+1
      do iu=1,nunptc
c       skip symmetry equivalent groups.
        if(nsymeq(iu).gt.0) then
          call monor(.false.,ifg,indat,fmozan,untxyz(1,1,iu+1),iabdfg,
     *               jabdfg,rij(noff),rwrk,vdwrad)
          noff=noff+nfg*nfg
        endif
      enddo
c     write(6,*) 'wwwrr',(rij(mfg2+i),i=1,nunptc)
      RETURN
      END
