C 11 Apr 08 - MWS - synchronize DFGRID common block
C  4 Mar 08 - DGF - FMOMINP: fix the RITRIM check
C  7 Dec 07 - DGF - change version number printing.
C 28 Aug 07 - DGF - small printing changes
C 20 Aug 07 - DGF - various changes for FMO 3.1 release
c 12 Jul 07 - MC  - add FMO-TDDFT arguments
C 24 Mar 07 - MWS - pad FRGINF common block
C 22 Dec 06 - TN  - synchronise EFPFMO common
C  8 Nov 06 - DGF - various changes for FMO 3.0 release
C  6 Nov 06 - MWS - adjust wavefunction and GDDI common block
C 22 Feb 06 - TN  - FMOPROP: include EFP/FMO model
C 21 Nov 05 - DGF - various changes for FMO 2.1 release
C 19 Sep 05 - IA  - synchronize FRGINF common
C  6 Jul 05 - DGF - FMOCCP: update print routine for CCSD(TQ)
C  1 Jun 05 - DGF - fixes for the 2nd release
c 31 Jan 05 - DGF,KK - add I/O module for the FMO method
c
C*MODULE fmoio   *DECK fmominp
      SUBROUTINE fmominp(itask,nder,ifgfmo0,ichfg,mulfg,frgnam,layfrg,
     *                   indat,scffrg,fmoscf,fmoci,fmodft,fmocc,mpnfmo,
     *                   fmotd,nacut,modmol,molfrg,nprfrg,orshft,lbody,
     *                   gcorrel,fmoq,iexcit,exfid)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 tstring
      logical GOPARR,DSKWRK,MASWRK,gcorrel,fmoq
      PARAMETER (nnam=37,MaxLay=5,ten=1.0D+01)
      DIMENSION QNAM(NNAM),KQNAM(NNAM),ichfg(*),mulfg(*),frgnam(*),
     *          layfrg(*),indat(*),scffrg(*),fmoscf(*),fmoci(*),
     *          fmodft(*),fmocc(*),mpnfmo(*),fmotd(*),molfrg(*),
     *          nprfrg(*),lbody(*),iexcit(5)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      EQUIVALENCE (tstring,dstring)
      DATA KQNAM /1,1,1,1,1, 5,1,23,23,3, 43,1,1,1,1, 1,5,5,5,5,
     *            1,1,3,1,1, 1,1,5,3,1,    1,1,1,1,51, 5,3/
      DATA QNAM/8HNFRAG   ,8HNLAYER  ,8HICHARG  ,8HMULT    ,8HINDAT   ,
     *          8HFRGNAM  ,8HLAYER   ,8HRESPAP  ,8HRESPPC  ,8HRESDIM  ,
     *          8HRITRIM  ,8HMAXKND  ,8HMAXCAO  ,8HNBSSE   ,8HMPLEVL  ,
     *          8HIXESP   ,8HDFTTYP  ,8HSCFFRG  ,8HCITYP   ,8HCCTYP   ,
     *          8HNACUT   ,8HNFRND   ,8HRCORSD  ,8HMOLFRG  ,8HNBODY   ,
     *          8HNOPFRG  ,8HMODMOL  ,8HSCFTYP  ,8HORSHFT  ,8HMODGRD  ,
     *          8HRESPCT  ,8HMAXBND  ,8HMODESP  ,8HIVMUL   ,8HIEXCIT  ,
     *          8HTDDFT   ,8HEXFID   /
      data fmogrp/8HFMO     /,rnone/8HNONE    /
c
c     this subroutine reads FMO input.
c     itask=1 light version: only read dimensions and options
c     itask=2 heavy version: read everything
c     Keywords in QNAM must be arranged appropriately.
c
c     atomic population approximation 
      RESPAP(1)=0
      RESPAP(2)=0
c     point charge approximation 
      RESPPC(1)=2.0D+00
      RESPPC(2)=0
c     separated dimer approximation (energy without SCF) 
      RESDIM=2.0D+00
c     restri(1,2,3) thresholds to ignore SCF trimer contribution
c     restri(4) a threshold to ignore correlated trimer contribution
c     The condition to NOT ignore trimer SCF contribution is:
c     rmin.le.restri(1).or.rmax.le.restri(2)).and.(rmax.le.
c    *            restri(3).or.restri(3).eq.0)).or.restri(1).eq.0
c     where rmin and rmax are returned by fmodist3(ifg,jfg,kfg,rmin,rmax)
c     The condition to NOT ignore trimer correlated contribution is:
c     rmax.le.restri(4).or.restri(4).eq.0
      RESPCT=0.0D+00
c
      defval=2.0D+00
      restri(1)=defval
      restri(2)=defval
      restri(3)=defval
      restri(4)=defval
c     separated dimer approximation (energy without correlation)
      rcorsd=2.0D+00
      nbsse=0
      IXESP=0
c     bit additive
c     1 add exchange terms to ESP, usually not needed.
c     2 remove delta-D terms(2nd/3rd order energy correction).
c       (That is, terms of the form Tr(delta-Dij*Vij) or Tr(delta-Dijk*Vijk).
c       Useful to make energy gradient analytic by changing the 
c       energy, not the gradient!
c     4 add derivative delta-D terms (not implemented)
c     32 do not add separated dimer derivative contributions. 
c 
      NACUT=0
c     slice the molecule into fragments each having NACUT atoms,
c     if NACUT is set by the user, INDAT is ignored. All other options must be
c     set by the user.
      nfmopal=0
c     "friend option" - provide compatibility with other implementations of FMO
c     nfmopal=1 GAUSSIAN-94, round off interfragment distances to 0.1 waals. 
c
c     the B parameter in the projection operators
      orshft=1.0d+06
c
c     n-body expansion of FMO 
      nbody=2
      if(ifgfmo0.ne.0) nbody=1
      call icopy(MaxLay,-1,0,lbody,1)
c
c     Selected monomer option (molfrg), bit additive.
c     1 skip SCF mode. If set, do not perform SCF calculations of non-selected
c       n-mers otherwise, perform correlated calculations of selected and SCF
c       of all n-mers (n>1).  
c     2 intra/intermolecular (set/unset) choice. Intramolecular means
c       do n-mers within the selected group, otherwise between selected and 
c       the rest.
c     4 save memory by not allocating space for n-mer arrays (active only
c       with intramolecular mode (2) and, at present, with n=3 only).
c     "Skip SCF" is not useful with nbody=3 and the current code. Perhaps in 
c     this case one should compute all dimers and then selected trimers with at
c     most one monomer from molfrg? 
c
      modmol=0
      modgrd=0
c     modgrd is bit additive.
c     1 do not subtract ESP from Lagrangian. This is needed to restore the
c       old published MP2 gradient.
c     2 add ESP derivatives (T. Nagata).
c     Options for ESP
c     modesp=0 the old uniform distances in dimers and trimers 
c     modesp=1 the n-body consistent distances except for connected n-mers 
c     modesp=2 the n-body consistent distances throughout 
c     modesp=4 additive dimer(trimer) ESP approximation (unconnected only)
c     modesp=8 use the modesp=4 approximation for connected cases.
c     modesp=16 FDIFF for ESPs
      modesp=1
c
c     The length of the ESP multipolar expansion.
      ivmul=1
c
c     only two dimension parameters are user accessible. 
c
c     at most 10 sets of LMOs
      MAXKND=10
c     at most 5 MOs in LMOs. 5 is for C (1s + four sp3). 
      maxcao=5
c     maxbas=1
      maxbnd=-1
c
      icurfg=0 
      jcurfg=0 
      kcurfg=0 
      ncursh=0
c     internal flags to alter computations of integrals 
c
c
c     initialize fmo-tddft options
      call viclr(iexcit,1,5) 
      if(TDDFTYP.ne.rnone) iexcit(2)=1
      iexcit(4)=-1
      exfid=0
c     NOPFRG:
c     print-out level separately for each fragment (bit additive). 
c     0 no additional print-out 
c     1 set the value of NPRINT=7 (print orbitals).
c     2 set MVOQ to +6 to obtain better virtual orbitals.
c     4 generate cube file, the grid is chosen automatically. 
c     8 remove one electron, donor fragment (meaningful only for nbody=1).
c     16 add one electron, acceptor fragment (meaningful only for nbody=1).
c     32 use 10 times tighter icut for DIRSCF jobs involving this fragment. 
c  
      if(itask.ne.1.and.itask.ne.2) then
        write(6,*) 'Internal Task error',itask
        call abrt
      endif
      KQNAM(25)=MaxLay*10+1
      if(itask.gt.1) then
        KQNAM(3)=nfg*10+1
        KQNAM(4)=nfg*10+1
        KQNAM(5)=(natfmo+nfg+1)*10+1
        KQNAM(6)=nfg*10+5
        KQNAM(7)=nfg*10+1
        KQNAM(15)=nlayer*10+1
        KQNAM(17)=nlayer*10+5
        KQNAM(18)=nfg*10+5
        KQNAM(19)=nlayer*10+5
        KQNAM(20)=nlayer*10+5
        KQNAM(24)=nfg*10+1
        KQNAM(26)=nfg*10+1
        KQNAM(28)=nlayer*10+5
        KQNAM(36)=nlayer*10+5
        do i=1,nfg
          ichfg(i)=0
          mulfg(i)=1
          layfrg(i)=1
          WRITE(UNIT=tstring,FMT='(A3,I5.5)') 'frg',i
          frgnam(i)=dstring 
          scffrg(i)=scftyp
          molfrg(i)=0
          nprfrg(i)=0
        enddo
        do i=1,nlayer
          mpnfmo(i)=mplevl
          fmoscf(i)=scftyp
          fmodft(i)=dftype
c         Although no check is done here, only grid-based DFT is supported. 
          fmocc(i)=cityp
          fmoci(i)=cityp
          fmotd(i)=TDDFTYP
        enddo
        do i=1,natfmo+nfg+1
          indat(1)=1 
        enddo
      else
        nfg=1
        nlayer=1
        KQNAM(3)=10*1+9
        KQNAM(4)=10*1+9
        KQNAM(5)=10*1+9
        KQNAM(6)=10*5+9
        KQNAM(7)=10*1+9
        KQNAM(15)=10*1+9
        KQNAM(17)=10*5+9
        KQNAM(18)=10*5+9
        KQNAM(19)=10*5+9
        KQNAM(20)=10*5+9
        KQNAM(24)=10*1+9
        KQNAM(26)=10*1+9
        KQNAM(28)=10*5+9
        KQNAM(36)=10*5+9
      endif
      CALL NAMEIO(IR,JRET,FMOGRP,NNAM,QNAM,KQNAM,
     *            nfg,nlayer,ichfg,mulfg,indat,frgnam,layfrg,rESPAP,
     *            RESPPC,RESDIM,restri,MAXKND,maxcao,nbsse,mpnfmo,IXESP,
     *            fmodft,scffrg,fmoci,fmocc,NACUT,nfmopal,rcorsd,molfrg,
     *            lbody,nprfrg,modmol,fmoscf,orshft,modgrd,respct,
     *            maxbnd,modesp,ivmul,iexcit, fmotd,exfid,0,0,0, 
     *            0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,  0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,*) 'Error reading $FMO'
         CALL ABRT
      END IF
      if(lbody(1).eq.-2) then
        lbody(1)=abs(lbody(1))
        fmoq=.true. 
      else
        fmoq=.false. 
      endif
c
      if(lbody(2).lt.0) call icopy(nlayer-1,lbody(1),0,lbody(2),1)
      if(lbody(1).lt.0) call icopy(nlayer,nbody,0,lbody,1)
      nbody=lbody(inamax(MaxLay,lbody,1))
      if(nbody.lt.0) nbody=0 
c     bizarre artefact: should use max value, not max abs. value in inamax. 
c
      if(nbody.gt.3) then
        if(maswrk) write(iw,9050) nbody 
        call abrt
      endif
      if(iexcit(2).gt.nbody) then
        if(maswrk) write(iw,9060)
        call abrt
      endif
      if(TDDFTYP.ne.rnone.and.iexcit(1).eq.0) then 
        if(maswrk) write(iw,*) 'TDDFT fragment undefined in iexcit.'
        call abrt
      endif 
      if(iexcit(4).lt.0) then
        iexcit(4)=0
        if(iexcit(2).ge.2) iexcit(4)=2
      endif
c     FMO-TDDFT gradient requires iexcit(3)=1. 
      if(TDDFTYP.ne.rnone.and.nder.gt.0) iexcit(3)=1 
c
      if(maxbnd.le.0) maxbnd=nfg*2+1
c
      if(nbody.le.2) call vclr(restri,1,4)
c
c     note that the 2nd elements are not used in Gaussian(?)
      if(iand(nfmopal,1).ne.0) then
        itmp=int(respap(1)*ten) 
        respap(1)=itmp/ten
        itmp=int(resppc(1)*ten) 
        resppc(1)=itmp/ten
        itmp=int(resdim*ten) 
        resdim=itmp/ten
      endif
      gcorrel=.false.  
      if(itask.gt.1) then
        do ilay=1,nlayer
c         gcorrel=gcorrel.or.mpnfmo(ilay).ne.0.or.fmoci(ilay).ne.rnone
c    *                   .or.fmocc(ilay).ne.rnone
c         Only one of these may be set. Check now if that is so.
          ncorme=0
          if(mpnfmo(ilay).ne.0) ncorme=ncorme+1
          if(fmoci(ilay).ne.rnone) ncorme=ncorme+1
          if(fmocc(ilay).ne.rnone) ncorme=ncorme+1
          if(fmotd(ilay).ne.rnone) ncorme=ncorme+1
          if(ncorme.gt.1) then
            write(6,*) 'Check MP,CI,CC and TD in layer ',ilay
            call abrt
          endif
          gcorrel=gcorrel.or.ncorme.ne.0
        enddo
      endif 
      if(.not.gcorrel) then
        rcorsd=0.0D+00
        restri(4)=0.0D+00
      endif
      if(gcorrel.and.rcorsd.eq.0.0D+00) then
        write(iw,*) 'Please set rcorsd to a nonzero value.'
        call abrt
      endif
      if(restri(4).gt.restri(3)) then
        write(iw,9040) restri(3),restri(4)
        call abrt
      endif
      if(iand(ixesp,1).ne.0.and.nder.gt.0) then
        write(iw,*) 'No gradient for exchange terms.' 
        call abrt
      endif
      if(iand(modgrd,2).ne.0) then
        if(respap(1).ne.0) then
          write(iw,*) 'ESP derivative (modgrd=2) requires respap=0.'
          call abrt
        endif
        if(nbody.lt.2) then
          write(iw,*) 'ESP derivative (modgrd=2) requires nbody>1.'
          call abrt
        endif
        if(iand(modesp,3).ne.0) then
          write(iw,*) 'ESP derivative (modgrd=2) requires modesp=0.'
          call abrt
        endif
      endif
c
c     FMO/BD always needs full allocation, since LMOs are regenerated.
c
      if(itask.gt.1.and.iand(ixesp,128).ne.0.and.maxcao.ne.maxcbs) then
        write(iw,*) 'Please set maxcao in $FMO to',maxcbs
        call abrt
      endif
c
      if(resdim.ne.0.0D+00.and.nbody.gt.2.and.nder.ne.0) then
        write(iw,*) '3-body gradient is not available with RESDIM.ne.0.'
        call abrt
c       While this seems to be easy, it is not so.
c       To implement one has to do this:
c       1. change grad accumulation to separate 1,2,3-body, 
c          rather than 1,1+2, 1+2+3 (easy).
c       2. Add RESDIM derivative terms to the 3-body gradient in TVDER, with
c          precomputed weights (not hard). 
c       3. Add RESDIM derivative terms to the 3-body gradient in FMOESP, with
c          precomputed weights (not hard). 
c       4. Accumulate different potentials VI-J (with precomputed weight
c          factors), to be used in a separate call to ESDGRD (messy). 
      endif
      if(gcorrel.and.rcorsd.ne.0.0D+00.and.nbody.gt.2.and.nder.ne.0.and.
     *   iand(ixesp,64).eq.0) then
        write(iw,*) '3-body gradient is not available with RCORSD.ne.0.'
c       Implementing this would require adding a loop 1,2 (=RHF,corr),
c       similar to monoscf. That is not very hard. 
c       Actually, the only way to run a 3-body correlated gradient is
c       to turn off ALL approximations!
        call abrt
      endif
c     if(ixesp.eq.-1) then
c       ixesp=0
c       if(nder.gt.0.and.resdim.gt.0) ixesp=32
c     endif
c     nfg is some estimate of the cross-cuts (extra bond cuts due to X-shaped
c     crossovers in proteins; S-S links, mostly).
c     mxatmf=mxatm*nfg
      maxvec=100
      maxabd=4
c     maxcbs will be set elsewhere
      return
c
 9040 format(1x,'Invalid trimer option ritrim(3)=',f10.5/
     *       1x,'      is smaller than ritrim(4)=',f10.5)
 9050 format(1x,'The only meaningful values for nbody are: 0,1,2,3',
     *          'and you have nbody=',I3/)
 9060 format(1x,'nbody should not be larger than iexcit(2).',/)
      end
c
C*MODULE fmoio   *DECK fmopinp
      SUBROUTINE fmopinp(itask,mannod,mcnvfmo,mconvex,mconfg,ijvec,
     *                   ngrfmo,prtdst,ncvscf,natdata,numdata,ibfconv,
     *                   mapconv,ngrmax,savgrd,loadbf,loadgr,l1dir,ngm,
     *                   vdwrad,grdpad,mpcmit,convpcm,pcmoff,imect,n0bda
     *                  ,ne0bda,r0bda,e0bda,efmo0,nepl0,epl0ds,eint0,
     *                   rappri,irestl)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK,ISGDDI,PAROUT,INITGDDI,savgrd,SG1
      PARAMETER (nnam=40,MAXNZ=137,nstfmo=10,maxpst=10,ONE=1.0D+00,
     *           MXGRID=10)
      DIMENSION QNAM(NNAM),KQNAM(NNAM),mannod(*),mcnvfmo(*),mconfg(*),
     *          ijvec(2,*),ngrfmo(maxpst,*),prtdst(3),ncvscf(2),
     *          ibfconv(*),mapconv(*),loadbf(maxpst,*),loadgr(maxpst,*),
     *          ngm(3),vdwrad(MAXNZ),vdwr0(MAXNZ),r0bda(*),
     *          e0bda(4,ne0bda,*),efmo0(*),epl0ds(*),eint0(4),rappri(3)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT,INITGDDI
      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
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      DATA KQNAM /1,1,3,1,1, 1,1,21,1,1, 93,1,33,1,3, 1,3,1,1,1, 
     *            1,1,1,1,31, 3,3,1,3,3, 1,1,1,3,3, 3,3,43,33,1/
      DATA QNAM/8HNGUESS  ,8HMAXIT   ,8HCONV    ,8HMODORB  ,8HMODPAR  ,
     *          8HMCONV   ,8HIREST   ,8HNCVSCF  ,8HIJVEC   ,8HNPRINT  ,
     *          8HESPSCA  ,8HNGRFMO  ,8HPRTDST  ,8HMANNOD  ,8HCNVDMP  ,
     *          8HMCONFG  ,8HCOROFF  ,8HIBFCON  ,8HMAPCON  ,8HMAXAOC  ,
     *          8HMODPRP  ,8HLOADBF  ,8HLOADGR  ,8HNAODIR  ,8HNGRID   ,
     *          8HVDWRAD  ,8HGRDPAD  ,8HNPCMIT  ,8HCNVPCM  ,8HPCMOFF  ,
     *          8HIMECT   ,8HIPIEDA  ,8HN0BDA   ,8HR0BDA   ,8HE0BDA   ,
     *          8HEFMO0   ,8HEPL0DS  ,8HEINT0   ,8HRAPPRI  ,8HIRESTL  /
      data fmogrp/8HFMOPRP  /
c     note: unknown parameters are set to 2.5A.
      data vdwr0/1.20d+00,1.20d+00,1.37d+00,1.45d+00,1.45d+00,1.50d+00,
     *          1.50d+00,1.40d+00,1.35d+00,1.30d+00,1.57d+00,1.36d+00,
     *          1.24d+00,1.17d+00,1.80d+00,1.75d+00,1.70d+00,12*2.5d+00,
     *          2.5d+00,4*2.5d+00,2.3d+00,102*2.5d+00/
c
c     this subroutine reads FMO property input.
c     care should be exercised about e0bda: it is used both as
c     e0bda(4,n0bda) and e0bda(4,nbdfg,nlayer)
c
      nguess=2
c     bit additive
c        1 use free molecule guess for monomers
c        2 use monomer density for dimer guess (if 0, use Huckel guess)
c        4 insert HMO projection operator in the Huckel guess 
c        8 apply dimer cut orbital projection to the dimer initial guess
c       16 do RHF for each dimer before MCSCF and match orbitals
c       32 a modifier of 16: copy dimer RHF core oorbitals(but use MCSCF active)
c       64 apply projection operators to dimer initial orbitals in the hope it
c          will remove the ill effect of the fractioned bonds upon convergence.
c      128 do not do restart during geometry optimisations (that is, always
c          use the same (Huckel) guess for monomers).
c      256 the "desperate" option: if SCF does not converge, try the
c          alternative converger (flip between SOSCF and DIIS).
c      512 manual orbital reordering, i.e. reorder initial orbitals using
c          $GUESS options (IORDER), useful for MCSCF.
c     1024 orthogonalise MCSCF dimer initial orbitals
c     2048 a modifier of 256: when going "desperate" do not read orbitals
c          of the previous (diverged) SCF.
c     Options 4,8,16,32,64 are effective only if some bonds are cut.
c     Options 64,128,512,1024 are obscure and almost never used.
c     Options 256 and 2048 may be useful in hard cases. 
      mxitfg=30
c     monomer SCF convergence (energy)
      convfg=1.0D-07
c     a threshold to force no static correlation (MCSCF/DFT) until energy 
c     converges up to coroff. If 0, this option is ignored.
c     This option is supposed to be useful mainly for DFT. Sensible values
c     are about 1e-2 ... 1e-4. Beside helping convergence, this option can 
c     speed the run up. 
c     coroff is not set if DFTTYP is not defined in $CONTRL.
c     This can be very confusing.  
      coroff=0.0d+00
      if(ndftfg.ne.0) coroff=1.0d-03
c
      prtdst(1)=0.0D+00
c     all interfragment distances less then prtdst(1) will be printed.
      prtdst(2)=0.5D+00
c     nonzero interfragment distances less then prtdst(2) will be marked.
c     fragments closer than prtdst(3) will be obliged to have a fractioned
c     bond defined in $FMOBND. -1.0D+00 effectively turns off such check.
c     0.6 is a good value to check C-C bonds.
c     The cost of checking will grow rapidly with prtdst(3) (not parallelised).
      prtdst(3)=0.6D+00
c 
c     bit additive
c     MODORB=0 exchange density between fragments
c     MODORB=1 exchange MOs between fragments
c     MODORB=2 orbital energies between fragments
c     For DFT MODORB=1 must be set (i.e., MODORB=1 or MODORB=3).
c     DFT will run with either of the two, but MODORB=1 does not properly
c     reorder the initial dimer orbitals, so MODORB=3 is the only proper choice.
c     MCSCF requires both orbitals and energies (MODORB=3 only).
      MODORB=0
      if(NDFTFG.ne.0) MODORB=3
c     if(NDFTFG.ne.0.or.SOSCF) MODORB=modorb+1
c     if(SOSCF)                modorb=modorb+2
c
c     MODPAR controls parallel options (bit-additive)
c     MODPAR=  1 turns on/off heavy job first strategy (aimed to reduce
c                waiting on remaining jobs at barrier points)
c     MODPAR=  2 changes ESP parallisation strategy:
c                if not set, parallelise loops over shells in each fragment,
c                otherwise parallelise loop over fragments.
c                Note that the latter way is seldom useful (i.e., many short-
c                contracted shells on slow network and few nodes). The other 
c                option is preferred, although it generates more parallel 
c                traffic to perform dynamic load balancing.
c     modpar=  4 broadcast all fragments done by a group at once rather
c                than fragment by fragment. 
c     modpar=  8 alters the behaviour of the fragment
c                initialisation: if set, fragments are always done in the
c                reverse order (nfg, nfg-1, ...1) because distance calculation 
c                costs decrease in the same order and they usually prevail over
c                making Huckel orbitals or running free monomer SCF. Note that
c                during SCC (monomer SCF) iterations the order in which monomers
c                are done is determined by MODPAR=1. 
c     modpar= 16 if set, LMO projectors (FMOHOP) will not be parallelised 
c     modpar= 32 if set, separated dimers will be done one per node, within
c                the group. This is often useful for MP2 when dimer group
c                size is large and there are many separated dimers. This option
c                requires MODPAR=1 be set and resdim be nonzero to be active.
c                The load balancing becomes bizarre: each new dimer is got with
c                either DLB or SLB at the group level but is assigned statically
c                to the group members irregardless of $CONTRL BALTYP. 
c                Nota bene! For restart jobs with irststp>=4 the data in $FMOEND
c                (if any) should be generated by a run with the same options
c                (that is, modpar=32), otherwise parallel deadlocks can occur. 
c     modpar= 64 Distribute file F40 for restarts irest>1. The file should be
c                put only on the grand master node and it will broadcast it to
c                all other nodes. Moreover, if the file F40 is copied to slaves
c                and modpar=64, then file errors will result due to 
c                improper indexing of the existing file.
c     modpar=128 flip fragment indices in separated dimers. It improves
c                parallelisation for massively parallel runs.
c                
      MODPAR=13
c
      call viclr(ngrfmo,1,maxpst*nlayer)
c     Set the number of GDDI groups:
c     ngrfmo(1) SCF monomers
c     ngrfmo(2) SCF dimers
c     ngrfmo(3) all trimers 
c     ngrfmo(4) correlated monomers 
c     ngrfmo(5) separated dimers
c     ngrfmo(6) SCF monomers after MCSCF monomers
c     ngrfmo(7) SCF dimers after MCSCF dimers 
c            8,9,10 reserved
c     These are defined inidivually for each layer.
c     The default is all 0s, using whatever is in $GDDI.
c  
c     mannod in $FMO similarly to $GDDI defines manual node division into 
c     groups separately for each ngrfmo(i) above.  
c     -1 defines automatic division. Note that to use manual division NGRFMO
c     must also be set.
      nprfmo=1
c     bits 1-2
c     0 normal output
c     1 reduced output
c     2 minimum output (not implemented yet)
c     4 print interfragment distances
c     8 print Mulliken charges
c
c     the ESP from the Fock matrix and rediagonalise.
      call dacopy(9,one,espsca,1)
      espscf=one
c     scale ESP during the first three iterations by multiplying by ESPSCA.
      cnvdmp=0
c
c     Set a converger for each step in FMO
      do mfmostp=1,nstfmo
        mcnvfmo(mfmostp)=-1
      enddo
      call viclr(loadbf,1,maxpst*nlayer)
      call viclr(loadgr,1,maxpst*nlayer)
c     loadbf and loadgr define semidynamic load balancing.
c     n-mers with more AOs than loadbf are done statically on the first 
c     loadgr groups, e.g. loadbf(1) and loadgr(1) are used for monomer SCF
c     (see ngrfmo(1:maxpst) description for explanation of steps). 
c     Other groups will use dynamic load balancing. When all static tasks
c     are done, dynamic load balancing will be used by the selected loadgr
c     groups as well. The following conditions must be satisfied in order to use
c     semidynamic load balancing: 
c     a) GDDI in use 
c     b) nonzero loadbf(ifmostp)+loadgr(ifmostp), 
c     c) $GDDI BALTYP=NXTVAL 
c     d) iand(modpar,1).ne.0 . 
c     It appears that ngrfmo(6) should be 0 as well (no MCSCF selected monomer
c     group division). 
c
c     rappri is in the order: rmin,rmax,rstep
      call vclr(rappri,1,3)
      if(itask.ne.1.and.itask.ne.2) then
        write(6,*) 'Internal Task error',itask
        call abrt
      endif
      KQNAM(6)=nstfmo*10+1
      KQNAM(12)=maxpst*nlayer*10+1
      KQNAM(22)=maxpst*nlayer*10+1
      KQNAM(23)=maxpst*nlayer*10+1
      KQNAM(26)=MAXNZ*10+1
      if(itask.gt.1) then
        KQNAM(9)=2*maxvec*10+1
        KQNAM(14)=ngrmax*maxpst*10+1
        KQNAM(16)=(nfg+mconvex*4)*10+1
        KQNAM(18)=2*natdata*10+1
        KQNAM(19)=natdata*maxl1c*10+1
        KQNAM(34)=10*ne0bda+3
        KQNAM(35)=10*4*ne0bda*nlayer+3
        KQNAM(36)=10*ne0bda*nlayer+3
        KQNAM(37)=10*nepl0+3
        do i=1,nfg+mconvex*4
          mconfg(i)=-1
        enddo
        do i=1,ngrmax*maxpst*nlayer
          mannod(i)=-1
        enddo
        call viclr(ijvec,1,2*maxvec)
        do i=1,2*natdata
          ibfconv(i)=0
        enddo
        do i=1,natdata*numdata
          mapconv(i)=-1
        enddo
        if(nbdfg.ne.0.and.ipieda.ne.0) then
          call vclr(r0bda,1,ne0bda)
          call vclr(e0bda,1,4*ne0bda*nlayer)
        endif
        if(ne0bda.ne.0) call vclr(efmo0,1,ne0bda*nlayer)
        if(nepl0.ne.0) call vclr(epl0ds,1,nepl0)
      else
        KQNAM(9)=10*1+9
        KQNAM(14)=10*1+9
        KQNAM(16)=10*1+9
        KQNAM(18)=10*1+9
        KQNAM(19)=10*1+9
        KQNAM(34)=10*3+9
        KQNAM(35)=10*3+9
        KQNAM(36)=10*3+9
        KQNAM(37)=10*3+9
      endif
c
c     restart options: step and where is step
c     at present only layer 1 steps can be restarted.
c     Supported options: 
c       irststp=2 resume monomer SCF 
c       irststp=4 resume dimer SCF 
c       irststp=5 add BSSE corrections
c     n-multiple of 100 can be added to skip n layers.
c     e.g. 702 means skip 7 layers (i.e. jump to layer 8), step 2.
c     all restart jobs require F40 with monomer densities.
c
      irststp=0
      irestl=0
c     switch the SCF converger after ncvscf(1) monomer SCF iterations (that is,
c     switch between SOSCF <-> FULLNR). Enforce the converger in MCONV(2) after
c     ncvscf(2), i.e. overwrite convertors in MCONFG. 
      ncvscf(1)=9999
      ncvscf(2)=9999
c  
      modprp=0
c     FMO properties:
c      1 electron density (AO-based matrix, in-core).
c      2 electron density (AO-based matrix, disk-based) - reserved.
c      4 electron density (on a grid, in-core); produces a Gaussian cube file. 
c      8 electron density (on a grid, in-core); produces a sparse cube file. 
c     16 automatically generate grid for modprp = 4 or 8.
c
      l1dir=0
c     parameter enforcing direct SCF if the current basis size is small.
c     l1dir is the maximum basis set size per 100 mln integrals.
c     The actual condition to enforce DIRSCF is:
c     (L1/l1dir)**4 * 10d+8 * k .gt. nintic*nproc
c     where L1 is the basis set size, k is 1.5 for LABSIZ=1,NWDVAR=2 or else 2.
c     nintic is the memory for the incore integrals buffer. 
c     The recommended value is 180 (based on experience).
      call viclr(ngm,1,3)
c     NGM defines the number of grid points for monomer cube files.
c     The choice of monomers for which cube files are to be computed 
c     is done in NPRFRG.
      call dcopy(MAXNZ,vdwr0,1,vdwrad,1)
c
c     Grid padding in a.u. Contributions to density on grid will be ignored
c     beyond the box surrounding an n-mer padded with grdpad from all sides.
c
      grdpad=2.0D+00
c
      mpcmit=30
      convpcm=1.0D-07
      pcmoff=0
c
c     charge transfer calculation method. 0,1,2,3 are supported.
c
      imect=4
c
      ipieda=0
c
      n0bda=0
      call vclr(eint0,1,4)
c
      CALL NAMEIO(IR,JRET,FMOGRP,NNAM,QNAM,KQNAM,
     *            NGUESS,mxitfg,convfg,modorb,modpar,mcnvfmo,irststp,
     *            ncvscf,ijvec,nprfmo,espsca,ngrfmo,prtdst,mannod,
     *            cnvdmp,mconfg,coroff,ibfconv,mapconv,maxl1c,modprp,
     *            loadbf,loadgr,l1dir,ngm,vdwrad,grdpad,mpcmit,convpcm,
     *            pcmoff,imect,ipieda,n0bda,r0bda,e0bda,efmo0,epl0ds,
     *            eint0,rappri,irestl, 0,0,0,0,0, 0,0,0,0,0,
     *            0,0,0,0,0, 0,0,0,0,0, 0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,*) 'Error reading $FMOPRP'
         CALL ABRT
      END IF
c
c     replace level default values by the general default 
      do mfmostp=1,nstfmo
c       -1 in fmoconv resets to defaults
       if(mcnvfmo(mfmostp).eq.-1) call fmoconv(mcnvfmo(mfmostp),-1,0,dum
     *                           ,.true.,.true.,.false.,.false.,.false.)
      enddo
      savgrd=iand(irststp,1024).ne.0
      if(irststp.ge.1024) irststp=irststp-1024
c     An option to save gradient after each dimer calculation for the purpose
c     of the following RUNTYP=GRADIENT restarts.
      if(irststp.gt.nstfmo.or.irestl.gt.nlayer) then
        write(iw,*) 'Invalid restart options',irststp,irestl
        call abrt
      endif
      if(iand(modpar,2).ne.0) then
        if(maswrk) write(iw,*) 'MODPAR=2 option has an unfixed bug.'
c       The bug has something to do with GOPARR set to false in parallel
c       execution; it shows up if NPROC>NFG-1??
c       call abrt
      endif 
c     if(.not.parout.and.n1.ne.0.and.(n1.ne.n2.or.n1.ne.ngroups)) then
c       if(maswrk) write(iw,9020) 
c     endif 
      if(itask.gt.1) then
        loop=0
        do k=1,nlayer
          do i=1,maxpst
            nnodi=0 
            do j=1,ngrfmo(i,k)
              loop=loop+1
              nnodi=nnodi+mannod(loop) 
            enddo
c           nnodi is negative if nodes are not divided manually
            if(nnodi.gt.0.and.nnodi.ne.nnglob.or.ngrfmo(i,k)
     *         .gt.nnglob) then
              if(maswrk) write(iw,9020) i,k,ngrfmo(i,k),nnodi,nnglob
              call abrt
            endif 
          enddo 
        enddo 
      endif
      if(ipieda.gt.0.and.iand(modorb,1).eq.0) then
        if(maswrk) write(iw,*) 'Setting modorb=1 flag...'
        modorb=modorb+1
      endif
c
c     if(ncvscf.eq.0) ncvscf=9999
c     if(ixesp.eq.-1) then
c       ixesp=0
c       if(nder.gt.0.and.resdim.gt.0) ixesp=32
c     endif
      return
 9020 format(/1x,'Wrong group division (ngrfmo/mannod) for step',I2,
     *           ', layer',I1,':',3I6,/)
      end
c
C*MODULE fmoio   *DECK fmocinp
      SUBROUTINE fmocinp(fmozan,fmoc,fmomas,dolat)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical GOPARR,DSKWRK,MASWRK,dolat
      PARAMETER (nnam=14)
      DIMENSION QNAM(NNAM),KQNAM(NNAM),u(3,3),ctr(3),
     *          fmozan(*),fmoc(3,*),fmomas(*)
      COMMON /FMCOM / X(1)
      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
      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)
      DATA KQNAM /1,33,33,43,33, 33,1,3,3,31, 31,31,31,1/
      DATA QNAM/8HIATORG  ,8HUNTANG  ,8HUNTORG  ,8HRESPBC  ,8HABCLAT  ,
     *          8HANGLAT  ,8HNSYMOP  ,8HSYMTRA  ,8HSYMOPE  ,8HIOPORG  ,
     *          8HIOPDIR  ,8HIOPABC  ,8HIOPANG  ,8HMAXKLM  /
      data fmogrp/8HFMOLAT  /
c
c     this subroutine reads lattice FMO input.
c
      tobohr=1.0D+00/0.52917724924D+00
      torad=acos(0.0D+00)/90.0D+00
c     iatorg: where to put the origin.
c              <0 at (0,0,0) 
c              =0 center of mass (default),
c              >0 at iatorg-th atom.
      KQNAM(8)=3*48*10+3
      KQNAM(9)=3*3*48*10+3
      iatorg=0
      nsymop=0
      do i=1,3
        untang(i)=0
        untorg(i)=0
        abclat(i)=0
        anglat(i)=0
        ioporg(i)=0
        iopdir(i)=0
        iopabc(i)=0
        iopang(i)=0
      enddo
      call vclr(respbc,1,4)
      call vclr(symtra,1,3*24)
      call vclr(symope,1,3*3*24)
      maxklms=10000
c
      CALL NAMEIO(IR,JRET,FMOGRP,NNAM,QNAM,KQNAM,
     *            iatorg,untang,untorg,respbc,abclat,anglat,nsymop,
     *            symtra,symope,ioporg,iopdir,iopabc,iopang,maxklms,0,
     *            0,0,0,0,0,  0,0,0,0,0, 0,0,0,0,0,
     *            0,0,0,0,0,  0,0,0,0,0, 0,0,0,0,0,
     *            0,0,0,0,0,  0,0,0,0,0, 0,0,0,0,0, 0,0,0,0)
      IF(JRET.EQ.2) THEN
         IF (MASWRK) WRITE(IW,*) 'Error reading $FMO'
         CALL ABRT
      END IF
      dolat=JRET.EQ.0
      if(.not.dolat) then
        maxklms=0
        return
      endif
c
      call dscal(3,torad,untang,1)
      call dscal(3,torad,anglat,1)
c     if(iuntrd.eq.1) call dscal(3,tobohr,abclat,1)
      call dscal(3,tobohr,abclat,1)
      call dscal(4,tobohr,respbc,1)
c
      vol=cellvol(abclat,anglat)
      if(maswrk) write(iw,9976) vol/tobohr**3
c     orientation of molecule
      call euleru(untang(1),untang(2),untang(3),u)
c     shift center-of-mass to origin
      iwhere=0
      if(iatorg.lt.0) then
        call vclr(ctr,1,3)
      else if(iatorg.eq.0) then
        call fndcntr(fmozan,fmoc,fmomas,iwhere,ctr(1),ctr(2),ctr(3))
      else
        call dcopy(3,fmoc(1,iatorg),1,ctr,1)
      endif
      do iat=1,natfmo
         xx=fmoc(1,i)-ctr(1)
         yy=fmoc(2,i)-ctr(2)
         zz=fmoc(3,i)-ctr(3)
         fmoc(1,i)=xx*u(1,1)+yy*u(1,2)+zz*u(1,3)
         fmoc(2,i)=xx*u(2,1)+yy*u(2,2)+zz*u(2,3)
         fmoc(3,i)=xx*u(3,1)+yy*u(3,2)+zz*u(3,3)
      enddo
c     check symmetry operations 
      CALL VALFM(LOADFM)
      lct=LOADFM+1
      last=lct+3*natfmo*2
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
      call chksymop(natfmo,fmoc,x(lct),x(lct+3*natfmo))
      CALL RETFM(NEED)
      return
 9976 format(' cell volume (A**3) = ',f20.4)
      end
c
C*MODULE fmoio   *DECK fmoxyz
      SUBROUTINE fmoxyz(gprnam,iscan,natfmo,fmozan,fmoc,izbas)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Parameter (MXATM=2000)
      character*99 STR,STRF77
      character*4 atm,c4dum
      dimension iargs(5)
      character*8 atname,ATOMNM,gprnam
      dimension fmoZAN(*),fmoC(3,*),izbas(*)
      logical GOPARR,DSKWRK,MASWRK,dolend
      COMMON /INFO  / Cdum(MXATM,3),IAN(MXATM),NATOMS,IUNTRD,ATM(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      data iargs/1,1,0,0,0/,c4dum/' '/
c
c     this subroutine reads FMO atomic coordinates.
c     parstat: GroupMasterBcast
c     iscan=0 normal reading and storing
c     iscan=1 scan to find NATFMO/maxbas, do not store fmozan, fmoc and izbas. 
c             In this case fmozan, fmoc and izbas should be declared in
c             the calling arguments as dum, dum3(3) and idum. 
c
      UNITS=1.0D+00/0.52917724924D+00
c     if(iand(nfmopal,1).ne.0) UNITS=1.889725989D+00
      if(iand(nfmopal,1).ne.0) UNITS=1.0D+00/0.529177249D+00
      CALL SEQREW(IR)
      CALL FNDGRP(IR,gprnam,IEOF)
      if(ieof.ne.0) then
        write(iw,*) gprnam,' not found' 
        call abrt
      endif
      natfmo=0
      if(maswrk) then
  100   continue
          natfmo=natfmo+1
          natfmo0=natfmo
          if(iscan.ne.0) natfmo0=1
c         if(natfmo.gt.MXATMF) then
c           write(iw,*) 'Too many atoms',natfmo,MXATMF 
c           call abrt
c         endif
         atname='        '
         READ(IR,9000,END=600,ERR=600) STR
         if(dolend(str)) goto 120
         call convsf77(80,str,strf77,5,iargs)
         READ(strf77,*,END=600,ERR=600) atname,ATOMNM,
     *                                  (fmoC(i,natfmo0),i=1,3)
c        READ(IR,*,END=120,ERR=600) atname,ATOMNM,(fmoC(i,natfmo),i=1,3)
          call UPRCAS(atname,8)
          call UPRCAS(ATOMNM,8)
c         if(atname.eq.'$END') goto 120
          if(iuntrd.eq.1) then
            do i=1,3
              fmoC(i,natfmo0)=fmoC(i,natfmo0)*units
            enddo
          endif
          imode=0
c         write(6,*) 'Reading',natfmo
c         write(6,9999) atname,ATOMNM 
c9999     format(2A8)
          call zsymnum(ATOMNM,c4dum,imode)
          fmoZAN(natfmo0)=imode
c         write(6,*) 'wwwatm',imode,(fmoC(i,natfmo0),i=1,3)
          ipos=ifndchr(atname,8,'.')
          if(ipos.ne.0) then
            read(UNIT=atname(ipos+1:8),FMT='(I3)') ibas
            if(ibas.le.0) then
              write(6,*) 'Invalid basis number in $FMOXYZ',ibas
              call abrt
            endif
          else
            ibas=1
          endif
          izbas(natfmo0)=ibas
          if(iscan.ne.0) maxbas=max(maxbas,ibas)
c         write(6,*) 'readbas',ibas
        goto 100
      endif
  120 continue
      IF (GOPARR) CALL DDI_BCAST(2400,'I',natfmo,1,MASTER)
      natfmo=natfmo-1
      if(natfmo.le.0) then
        write(iw,*) 'No atomic coordinates found' 
        call abrt
      endif
      if(iscan.eq.0.and.GOPARR) then
        CALL DDI_BCAST(2401,'F',fmoC,3*natfmo,MASTER)
        CALL DDI_BCAST(2402,'F',fmoZAN,natfmo,MASTER)
        CALL DDI_BCAST(2403,'I',izbas,natfmo,MASTER)
      endif
      if(iscan.eq.1.and.GOPARR) then
        CALL DDI_BCAST(2401,'I',maxbas,1,MASTER)
      endif
      return
  600 continue
c     call UPRCAS(atname,8)
c     if(atname.eq.'$END') goto 120
      write(6,*) 'error reading ',gprnam,', check atom',natfmo
c     now kill slave nodes by broadcasting reset natfmo
      natfmo=1
      CALL DDI_BCAST(2400,'I',natfmo,1,MASTER)
      call abrt
      return
 9000 format(A80)
      end
C
C*MODULE fmoio   *DECK fmolmo
      SUBROUTINE fmolmo(gprnam1,gprnam2,taotyp,iaprjo,japrjo,CoreAO,
     *                  OccCor,nCBS,nCAO,shiftb,libish,orshft,modQbas)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*99 STR,STRF77
      dimension iargs(3),mCAO0(2)
      character*8 gprnam1,gprnam2,taotyp(*)
      logical GOPARR,DSKWRK,MASWRK,modQbas,dolend
      dimension iaprjo(MaxCAO,*),japrjo(MaxCAO,*),OccCor(MaxCAO,*),
     *          nCBS(*),nCAO(*),CoreAO(MaxCBS,MaxCAO,*),shiftb(MaxCAO,*)
     *         ,libish(*)
      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
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      data iargs/1,0,0/
c
c     this subroutine reads localised MOs for boundary orbitals.  
c     parstat: GroupMasterBcast
c     Two spellings of this group name are supported. 
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,gprnam1,IEOF)
      if(ieof.ne.0) then 
        CALL SEQREW(IR)
        CALL FNDGRP(IR,gprnam2,IEOF)
        if(ieof.ne.0) return 
      endif 
      if(modQbas) call viclr(nCAO,1,maxknd*3)
      ierr=0
      naotyp=0
      if(maswrk) then
  100   continue
          naotyp=naotyp+1
          if(naotyp.gt.MaxKnd) then
            write(iw,*) 'Two many LMO types',naotyp,MaxKnd
            call abrt
          endif
          ierr=0
          taotyp(naotyp)='        '
          READ(IR,9000,END=600,ERR=600) STR
          if(dolend(str)) goto 120
          call convsf77(80,str,strf77,3,iargs)
          READ(strf77,*,END=600,ERR=600) taotyp(naotyp),mCBS,mCAO
c         READ(IR,*,END=120,ERR=600) taotyp(naotyp),mCBS,mCAO
          call UPRCAS(taotyp(naotyp),8)
c         if(taotyp(naotyp).eq.'$END') goto 120
          if(mCBS.gt.MaxCBS.or.mCAO.gt.MaxCAO) then
             write(iw,*) 'Two many AOs/MOs in $FMOLMO',mCBS,MaxCBS,mCAO,
     *                   MaxCAO 
             call abrt
          endif
c         The else clause covers some dormant option to be activated later.
          if(mod(nguess,2).ne.333) then
          READ(IR,*,END=600,ERR=600) (iaprjo(i,naotyp),japrjo(i,naotyp),
     *                          (CoreAO(j,i,naotyp),j=1,mCBS),i=1,mCAO)
          else
          READ(IR,*,END=600,ERR=600) (iaprjo(i,naotyp),japrjo(i,naotyp),
     *                               (CoreAO(j,i,naotyp),j=1,mCBS),
     *                                OccCor(i,naotyp),i=1,mCAO)
          endif
c         READ(IR,*,END=600,ERR=600) (iaotyp(j,naotyp),j=1,mCBS) 
          nCBS(naotyp)=mCBS
          nCAO(naotyp)=mCAO
          call dacopy(mCAO,orshft,shiftb(1,naotyp),1)
          if(modQbas) then
            ishiftl=naotyp+maxknd
            ishiftr=naotyp+maxknd*2
            mCAO0(1)=ishiftl
            mCAO0(2)=ishiftr
            call convlmo(mCBS,mCAO,mCAO0,naotyp,iaprjo(1,naotyp),
     *                   japrjo(1,naotyp),CoreAO,libish)
            nCAO(ishiftl)=mCAO0(1)
            nCAO(ishiftr)=mCAO0(2)
          endif
          ierr=1
          if(mCBS.lt.mCAO) goto 600
        goto 100
      endif
  120 continue
      IF (GOPARR) CALL DDI_BCAST(2403,'I',naotyp,1,MASTER)
      naotyp=naotyp-1 
      if(naotyp.le.0) then
        if(maswrk) write(iw,*) 'No LMOs found'
        if(naotyp.lt.0) then
          write(6,*) 'Invalid naotyp',naotyp
          call abrt
        endif
      endif
      IF (GOPARR) then
        CALL DDI_BCAST(2404,'F',taotyp,naotyp,MASTER)
        CALL DDI_BCAST(2405,'I',iaprjo,MaxCAO*naotyp,MASTER)
        CALL DDI_BCAST(2406,'I',japrjo,MaxCAO*naotyp,MASTER)
        naotyp3=naotyp
        if(modQbas) naotyp3=maxknd*3 
        CALL DDI_BCAST(2407,'F',CoreAO,MaxCBS*MaxCAO*naotyp3,MASTER)
        if(mod(nguess,2).eq.1) 
     *  CALL DDI_BCAST(2408,'F',OccCor,MaxCAO*naotyp,MASTER)
        CALL DDI_BCAST(2409,'I',nCBS,naotyp,MASTER)
        CALL DDI_BCAST(2410,'I',nCAO,naotyp3,MASTER)
        CALL DDI_BCAST(2411,'F',shiftb,MaxCAO*naotyp,MASTER)
      endif
c     write(6,*) 'read ',naotyp,' LMO types'
      return
  600 continue
c     call UPRCAS(taotyp(naotyp),8)
c     if(taotyp(naotyp).eq.'$END') goto 120
      write(6,*) 'error reading LMO: check input',naotyp
      if(ierr.eq.0) write(6,*) 'format error (or insufficient data).' 
      if(ierr.eq.1) write(6,*) 'more MOs than AOs.'
c     if(ierr.eq.2) write(6,*) 'left-right assignment error.'
c     now kill slave nodes by broadcasting reset naotyp
      naotyp=0
      CALL DDI_BCAST(2403,'I',naotyp,1,MASTER)
      call abrt
      return
 9000 format(A80)
      end
C
C*MODULE fmoio   *DECK fmobon
      SUBROUTINE fmobon(gprnam,taotyp,taotypi,iabdfg,jabdfg,idxcao)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (MAXL=5)
      character*99 STR,STRF77
      dimension iargs(2+MAXL+1)
      character*8 gprnam,taotyp(*),taotypi(*)
      logical GOPARR,DSKWRK,MASWRK,dolend
      dimension iabdfg(*),jabdfg(*),idxCAO(MaxBnd,*)
      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
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
c
c     this subroutine reads boundary orbital indices.  
c     parstat: GroupMasterBcast
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,gprnam,IEOF)
      if(ieof.ne.0) return 
      iargs(1)=0
      iargs(2)=0
      do i=1,nlayer+1
        iargs(2+i)=1
      enddo
      nbdfg=0
      if(maswrk) then
  100   continue
          nbdfg=nbdfg+1
          if(nbdfg.gt.MaxBnd) then 
             write(iw,*) 'Two many cut bonds',nbdfg,MaxBnd
             call abrt
          endif
          READ(IR,9000,END=600,ERR=600) STR
          if(dolend(str)) goto 120
          call convsf77(80,str,strf77,2+nlayer+1,iargs)
          READ(strf77,*,END=600,ERR=600) iabdfg(nbdfg),jabdfg(nbdfg),
     *                                    (taotypi(i),i=1,nlayer+1)
c         READ(IR,*,END=120,ERR=120) iabdfg(nbdfg),jabdfg(nbdfg),
c    *                               (taotypi(i),i=1,nlayer+1)
          do i=1,nlayer+1
            call UPRCAS(taotypi(i),8)
          enddo
c         reorder indices in the canonical way: first negative
          if(jabdfg(nbdfg).lt.0) then
             itmp=jabdfg(nbdfg)
             jabdfg(nbdfg)=iabdfg(nbdfg)
             iabdfg(nbdfg)=itmp
             if(maswrk) write(6,*) 'Warning: indices reordered',nbdfg
c            This is nothing important (is it?) 
          endif
          do i=1,nlayer+1
            idxcao(nbdfg,i)=ifndtxt(naotyp,taotyp,taotypi(i))
          enddo
        goto 100
      endif
  120 continue
      IF (GOPARR) then
         CALL DDI_BCAST(2410,'I',nbdfg,1,MASTER)
         CALL DDI_BCAST(2411,'I',iabdfg,nbdfg,MASTER)
         CALL DDI_BCAST(2412,'I',jabdfg,nbdfg,MASTER)
         CALL DDI_BCAST(2413,'I',idxcao,MaxBnd*(nlayer+1),MASTER)
      endif
      nbdfg=nbdfg-1 
      if(nbdfg.le.0) then
        if(maswrk) write(iw,*) 'No bond cuts found.'
c       call abrt
      endif
      return
  600 continue
      write(6,*) 'error reading ',gprnam,' check entry ',nbdfg
      call abrt
      return
 9000 format(A80)
      end
c
C*MODULE fmoio   *DECK convsf77
      subroutine convsf77(lens,sinp,sout,nargs,iargs)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      character*99 sinp,sout
      character*1 quote,space
      dimension iargs(nargs) 
      data quote/''''/,space/' '/
c
c     Convert human string format into Fortran77 format by quoting all
c     string arguments listed in iargs. "99" above is obtained as 80+
c     at most 9 string arguments (since each adds 2 characters).
c     Only spaces are allowed to separate arguments (shame on tab users).
c     sinp: input string of length lens (<=80)
c     sout: output string
c     nargs: the expected number of space separated words in sinp
c     iargs(i) tells if word i is a string (iargs(i).ne.0).
c     nargs and iargs are arguments to convsf77.
c     Multiple space separation is allowed.
c
      if(nargs.gt.9) call abrt
      io=0
      iarg=0
      inside=0
      do i=1,lens
c
        if(sinp(i:i).eq.space.and.inside.ne.0) then
          inside=0
          if(iarg.le.nargs.and.iarg.gt.0) then
            if(iargs(iarg).ne.0) then
              io=io+1
              sout(io:io)=quote
c             end quote
            endif
          endif
        endif
c
        if(sinp(i:i).ne.space.and.inside.eq.0) then
          iarg=iarg+1
          inside=1
          if(iarg.le.nargs) then
            if(iargs(iarg).ne.0) then
              io=io+1
              sout(io:io)=quote
c             begin quote
            endif
          endif
        endif
c
        io=io+1
        sout(io:io)=sinp(i:i)
c       copy the rest, including possible NULL and garbage characters
      enddo
c     append end quote if needed at the very end
      if(inside.ne.0.and.iarg.le.nargs.and.iarg.gt.0) then
        if(iargs(iarg).ne.0) then
          io=io+1
          sout(io:io)=quote
c         end quote
        endif
      endif
      RETURN
      END
c
C*MODULE fmoio   *DECK ifndtxt
      function ifndtxt(n,text,texti) 
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 text(n),texti
c
      ifndtxt=0
      do i=1,n
         if(text(i).eq.texti) then
            ifndtxt=i
            return
         endif
      enddo
      write(6,*) texti,' not found in the stored area.'
      call abrt
      return
      end
C
C*MODULE fmoio   *DECK zsymnum
      subroutine zsymnum(atninp,atnout,imode)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (MXEL=137)
      character*8 atninp
      character*4 atnout,atmlab(MXEL)
C
      DATA ATMLAB/'H   ','HE  ','LI  ','BE  ','B   ','C   ',
     *            'N   ','O   ','F   ','NE  ','NA  ','MG  ',
     *            'AL  ','SI  ','P   ','S   ','CL  ','AR  ',
     *            'K   ','CA  ','SC  ','TI  ','V   ','CR  ',
     *            'MN  ','FE  ','CO  ','NI  ','CU  ','ZN  ',
     *            'GA  ','GE  ','AS  ','SE  ','BR  ','KR  ',
     *            'RB  ','SR  ','Y   ','ZR  ','NB  ','MO  ',
     *            'TC  ','RU  ','RH  ','PD  ','AG  ','CD  ',
     *            'IN  ','SN  ','SB  ','TE  ','I   ','XE  ',
     *            'CS  ','BA  ','LA  ','CE  ','PR  ','ND  ',
     *            'PM  ','SM  ','EU  ','GD  ','TB  ','DY  ',
     *            'HO  ','ER  ','TM  ','YB  ','LU  ','HF  ',
     *            'TA  ','W   ','RE  ','OS  ','IR  ','PT  ',
     *            'AU  ','HG  ','TL  ','PB  ','BI  ','PO  ',
     *            'AT  ','RN  ','FR  ','RA  ','AC  ','TH  ',
     *            'PA  ','U   ','NP  ','PU  ','AM  ','CM  ',
     *            'BK  ','CF  ','ES  ','FM  ','MD  ','NO  ',
     *            'LR  ','RF  ','HA  ','SG  ','NS  ','HS  ',
     *            'MT  ','110 ','111 ','112 ','113 ','114 ',
     *            '115 ','116 ','117 ','118 ','119 ','120 ',
     *            '121 ','122 ','123 ','124 ','125 ','126 ',
     *            '127 ','128 ','129 ','130 ','131 ','132 ',
     *            '133 ','134 ','135 ','136 ','137 '/
c
c     This subroutine either finds the atomic number given a text name (imode=0)
c     or returns the text name given an atomic number imode. 
C
      if(imode.eq.0) then
c       first try to see if atninp already contains an integer, possibly
c       followed by a decimal dot and some other numbers to be ignored.
        ipos=ifndchr(atninp,8,'.')
        if(ipos.eq.0) ipos=9
        read(UNIT=atninp(1:ipos-1),FMT='(I4)',ERR=100) iat
        imode=iat
        return
100     continue
        do i=1,MXEL
          if(atmlab(i).eq.atninp(1:4)) then
            imode=i
            return
          endif
        enddo
        write(6,*) 'Unknown element: ',atninp
        call abrt
      else
        atnout=atmlab(imode)
      endif
      return
      end
c
C*MODULE fmoio   *DECK fmoout
      subroutine fmoout(need,lbody,fmozan,mulfg,scffrg,fmoscf,fmoci,
     *                  fmodft,fmocc,mpnfmo,fmotd,modmol,molfrg,nstfmo,
     *                  ncvscf,mcnvfmo,ngrfmo,loadbf,loadgr,nstjob,
     *                  ichtot,nacut,savgrd,nxg,nyg,nzg,ngm,vdwrad,
     *                  grdpad,spargrid,nxyzg,ifgfmo0,orshft,mpcmit,
     *                  iexcit)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter (one=1.0D+00,maxpst=10)  
      logical savgrd,spargrid
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,nfmopcm
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      common /fmodim/ maxbnd,maxknd,maxcbs,maxcao,maxvec,maxl1,
     *                maxnat,maxabd,maxbas
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension lbody(*),fmozan(*),mulfg(*),scffrg(*),fmoscf(*),fmoci(*)
     *         ,fmodft(*),fmocc(*),mpnfmo(*),fmotd(*),molfrg(*),
     *          ncvscf(2),mcnvfmo(*),ngrfmo(maxpst,*),loadbf(maxpst,*),
     *          loadgr(maxpst,*),nstjob(*),ngm(3),vdwrad(*),iexcit(5)
c
      write(iw,8500)
      if(ifgfmo0.le.1) write(iw,8510)
      write(iw,8520)
      write(iw,9000) need,nfg,nlayer,nbody,natfmo,ichtot,maxbas,naotyp,
     *               nbdfg,convfg,maxl1
      if(nbdfg.ne.0) write(iw,8990) orshft
      if(ifgfmo0.gt.0) write(iw,9005) ifgfmo0
      nmols=0
      do ifg=1,nfg
        if(mulfg(ifg).ne.1) write(iw,8070) ifg,mulfg(ifg)
        if(scffrg(ifg).ne.scftyp) write(iw,8060) ifg,scffrg(ifg)
c       if(molfrg(ifg).ne.0) nmols=ifg
        if(molfrg(ifg).gt.0) nmols=ifg
      enddo
      nsca=0 
      do i=1,9
        if(espsca(i).ne.one) nsca=i
      enddo
      write(iw,8000)
      write(iw,8002) (lbody(i),i=1,nlayer)
      write(iw,8005) (fmoscf(i),i=1,nlayer)
      write(iw,8010) (mpnfmo(i),i=1,nlayer)
      write(iw,8020) (fmodft(i),i=1,nlayer)
      write(iw,8030) (fmoci(i),i=1,nlayer)
      write(iw,8040) (fmocc(i),i=1,nlayer)
      write(iw,8050) (fmotd(i),i=1,nlayer)
c     do ilay=1,nlayer
      if(nfmopcm.ne.0.and.mpcmit.ne.2) write(iw,8080) nfmopcm 
      if(nfmopcm.ne.0.and.mpcmit.eq.2) write(iw,8082) 1,nfmopcm
      if(iexcit(1).ne.0)  write(iw,8090) (iexcit(i),i=1,4) 
c
      write(iw,8100)
      if(respap(1).gt.0) write(iw,9010) respap(1)
      if(respap(1).lt.0) write(iw,9011) 
      if(respap(2).gt.0) write(iw,9015) respap(2) 
      if(respap(2).lt.0) write(iw,9016) 
      if(resppc(1).gt.0) write(iw,9020) resppc(1)
      if(resppc(1).lt.0) write(iw,9021) 
      if(resppc(2).gt.0) write(iw,9025) resppc(2)
      if(resppc(2).lt.0) write(iw,9026) 
      if(nbody.gt.1) then  
        if(resdim.gt.0) write(iw,9030) resdim
        if(resdim.lt.0) write(iw,9031) 
        if(rcorsd.gt.0) write(iw,9035) rcorsd
        if(rcorsd.lt.0) write(iw,9036) 
        if(nbody.gt.2) then  
         if(restri(1).gt.0) write(iw,9032) restri(1),restri(2),restri(3)
         if(restri(4).gt.0) write(iw,9033) restri(4)
        endif
      endif
      if(iand(IXESP,2).ne.0) write(iw,9091) 
      if(iand(IXESP,32).ne.0) write(iw,9095)
      if(iand(modESP,7).eq.0) write(iw,9096)
      if(iand(modESP,7).eq.1) write(iw,9097) 
      if(iand(modESP,7).eq.2) write(iw,9098) 
c     if(ifgdon.ne.0.and.ifgacc.ne.0) write(iw,9097) ifgdon,ifgacc
      call vdwrout(fmozan,vdwrad) 
      if(iand(modgrd,1).ne.0) write(iw,9040) 
      if(iand(modgrd,2).ne.0) write(iw,9042) 
c
      write(iw,8110)
      if(mod(nguess,  2).eq.1) write(iw,9070) 
      if(mod(nguess/2,2).eq.1) write(iw,9071) 
      if(mod(nguess/4,2).eq.1) write(iw,9072) 
      if(mod(nguess/8,2).eq.1) write(iw,9073) 
      if(mod(nguess/16,2).eq.1) write(iw,9074)
      if(mod(nguess/32,2).eq.1) write(iw,9075)
      if(mod(nguess/128,2).eq.1) write(iw,9077)
      if(iand(nguess,512).ne.0) write(iw,9079)
      if(iand(nguess,32768).ne.0) write(iw,9080)
      if(irststp.ne.0) write(iw,9120) irstlay,irststp,modrst
      if(savgrd) write(iw,9121)
c
      write(iw,8120)
      if(nbsse.eq.1) write(iw,9085) 
      if(nbsse.eq.2) write(iw,9086) 
      if(nbsse.eq.3) write(iw,9087) 
      if(iand(IXESP,1).ne.0) write(iw,9090) 
      if(iand(nfmopal,1).ne.0) write(iw,9150)
      if(nmols.ne.0) then
        if(iand(modmol,1).ne.0) then
          write(iw,9171) 
        else
          if(rcorsd.ne.0) write(iw,9170) 
        endif
        if(iand(modmol,2).ne.0) then
          write(iw,9172) (molfrg(i),i=1,nmols)
        else
          write(iw,9173) (molfrg(i),i=1,nmols)
        endif
        if(iand(modmol,4).ne.0.and.nbody.gt.2) write(iw,9174)
      endif
c
      write(iw,8130)
      if(nacut.ne.0) write(iw,9300) nacut 
c
      write(iw,8140)
      write(iw,9001) (mcnvfmo(i),i=1,nstfmo) 
      if(ncvscf(1).ne.9999) write(iw,9002) ncvscf(1) 
      if(ncvscf(2).ne.9999) write(iw,9003) ncvscf(2) 
      if(mod(modorb,  2).eq.0) write(iw,9100) 
      if(mod(modorb,  2).ne.0) write(iw,9101) 
      if(mod(modorb/2,2).ne.0) write(iw,9102) 
      if(nsca.gt.0) write(iw,9004) (espsca(i),i=1,nsca)
      if(coroff.ne.0) write(iw,9006) coroff
c
      write(iw,8150)
      if(mod(modpar,  2).ne.0) write(iw,9110) 
      if(mod(modpar/2,2).eq.0) write(iw,9111) 
      if(mod(modpar/2,2).ne.0) write(iw,9112) 
      if(iand(modpar,4).ne.0) write(iw,9113) 
      if(iand(modpar,8).ne.0) write(iw,9114) 
      if(iand(modpar,16).ne.0) write(iw,9115) 
      if(iand(modpar,64).ne.0) write(iw,9118) 
      if(iand(modpar,128).ne.0) write(iw,9119) 
      do i=1,nlayer  
        if(ngrfmo(1,i).ne.0) write(iw,9130) i,ngrfmo(1,i),ngrfmo(4,i),
     *                                        ngrfmo(6,i)
        if(ngrfmo(2,i).ne.0) write(iw,9132) i,ngrfmo(2,i),ngrfmo(5,i),
     *                                      ngrfmo(7,i)
        if(ngrfmo(3,i).ne.0) write(iw,9134) i,ngrfmo(3,i)
        if(nstjob(1).ne.0) write(iw,9140) i,nstjob(1),' monomer',
     *                                    loadbf(1,i),loadgr(1,i)
      enddo
c
      write(iw,8160)
      if(iand(nprfmo,3).eq.0) write(iw,9200) 
      if(iand(nprfmo,3).eq.1) write(iw,9202) 
      if(iand(nprfmo,4).ne.0) write(iw,9206) 
      if(iand(nprfmo,8).ne.0) write(iw,9208) 
      if(iand(modprp,1).ne.0) write(iw,9210)
      if(iand(modprp,4).ne.0) write(iw,9216) nxg,nyg,nzg 
      if(spargrid) write(iw,9220) (nxg*nyg*nzg*1.0D+02)/nxyzg,nxyzg
      if(ngm(1)*ngm(2)*ngm(3).ne.0) write(iw,9222) ngm(1),ngm(2),ngm(3)
      if(ngm(1)*ngm(2)*ngm(3).ne.0.or.iand(modprp,4).ne.0.or.spargrid) 
     *  write(iw,9230) grdpad 
      if(ipieda.eq.2) write(iw,9240) 
      if(ipieda.eq.1) write(iw,9250) 
c
      write(iw,9500) maxbnd,maxknd,maxcbs,maxcao
      write(iw,9900) 
c
 8000 format(1x,'Layer electron correlation information')
 8002 format(4x,' nbody=',8(1x,I2,5x))
 8005 format(4x,'SCFTYP=',8A8)
 8010 format(4x,'MPLEVL=',8(1x,I2,5x))
 8020 format(4x,'DFTTYP=',8A8)
 8030 format(4x,' CITYP=',8A8)
 8040 format(4x,' CCTYP=',8A8/)
 8050 format(4x,' TDDFT=',8A8/)
 8060 format(4x,'Fragment ',I5,' has SCFTYP=',A8)
 8070 format(4x,'Fragment ',I5,' has MUL=',I4)
 8080 format(1x,'FMO/PCM[',I1,'] model will be used.')
 8082 format(1x,'FMO/PCM[',I1,'(',I1,')] model will be used.')
 8090 format(1x,'Excited state options (TDDFT)',
     *       /5x,'Excited fragment:         ',I5,
     *       /5x,'Many-body excitation level:  ',I2,
     *       /5x,'Energy option:               ',I2, 
     *       /5x,'Dimer matching scheme:       ',I2) 
 8100 format(/1x,'FMO approximations.')
 8110 format(/1x,'FMO initial guess options.')
 8120 format(/1x,'FMO extensions.')
 8130 format(/1x,'FMO automation options.')
 8140 format(/1x,'FMO SCF converger options.')
 8150 format(/1x,'FMO parallelisation options.')
 8160 format(/1x,'FMO properties and print-out level.')
 8500 format(/1x,71(1H-),
     *      /18x,'The Fragment Molecular Orbital (FMO) method.')
 8510 format(23x,'The program is developed at the',
     *      /11x,'Research Institute for Computational Sciences (RICS)',
     *       /1x,'National Institute of Advanced Industrial Science ',
     *           'and Technology (AIST)',
     *       /5x,'AIST Tsukuba Central 2, Umezono 1-1-1, Tsukuba, ',
     *           '305-8568, Japan.',
     *      //30x,'Version 3.1',
     *       /20x,'coded by D. G. Fedorov and K. Kitaura,',
     *       /16x,'with contributions from M. Chiba and T. Nagata.',
     *       /7x,'All publications using the FMO method in GAMESS ',
     *           'should reference:'
     *       /11x,'D.G. Fedorov, K. Kitaura, J. Chem. Phys. ',
     *           '120, 6832 (2004).')
 8520 format(1x,71(1H-))
 9000 format(/1x,'Used memory (words): ',I8,
     *       /1x,'Number of fragments:',I6,/1x,'Number of layers:',I9,
     *       /1x,'N-body FMO method:  ',I6,/1x,'Number of atoms: ',I9,
     *       /1x,'Total charge: ',I12,/1x,'Number of basis sets:',I5,
     *       /1x,'Number of LMO types:',I6,/1x,'Number of boundaries:',
     *    I5,/1x,'MonomerSCF convergence:',E12.4,
     *       /1x,'Max AOs per frg:',I10)
 8990 format( 1x,'Projection operator constant:',E10.2)
 9001 format(4x,'SCF convergers:',/16I7/)
 9002 format(4x,'SCF convergers will be switched after',I3,
     *          ' monomer SCF iterations.')
 9003 format(4x,'SCF convergers will be reset to mconv(2) after',I3,
     *          ' monomer SCF iterations.')
 9004 format(4x,'SCF ESP scaling at initial iterations:',9F5.2/)
 9005 format(1x,'Free monomer:    ',I9)
 9006 format(4x,'DFT will be turned off until monomer SCF converges to',
     *          E11.4)
 9010 format(4x,'Using AO pop approximation for the 2e ESP with the VdW'
     *         ,' factor',F7.2)
 9011 format(4x,'Using AO pop approximation for the 2e ESP of all ',
     *           'fragments.')
 9015 format(4x,'Using AO pop approximation for separated ',
     *          'dimers with the VdW factor',F7.2)
 9016 format(4x,'Using AO pop approximation for the energy of all ',
     *          'dimers .')
 9020 format(4x,'Using point charge approximation for the 2e ESP ',
     *           'with the VdW factor',F7.2)
 9021 format(4x,'Using point charge approximation for the 2e ESP ',
     *           'of all fragments.')
 9025 format(4x,'Using point charge approximation for the energy of ',
     *           'of separated dimers with the VdW factor',F7.2)
 9026 format(4x,'Using point charge approximation for the energy of ',
     *           'all dimers.')
 9030 format(4x,'Using ES-dimer approximation with the VdW factor',F7.2)
 9031 format(4x,'Using ES-dimer approximation for all fragments.')
 9032 format(4x,'Ignore trimers using thresholds T1, T2 and T3 if ',
     *      /4x,'RD>T1 and RM>T2 or RM>T3',3F7.2,/4x,'where RM is the ',
     *       'monomer-dimer distance and the monomer is chosen so that',
     *      /4x,'the remaining dimer has the shortest interfragment',
     *      /4x,'distance RD among all 3 dimers in the trimer.') 
 9033 format(4x,'Ignore correlation for trimers separated by the VdW ',
     *          'factor',F7.2)
 9035 format(4x,'No dynamic correlation for dimers with the VdW',
     *          ' factor ',F7.2)  
 9036 format(4x,'No dynamic correlation for all dimers.')
 9040 format(4x,'FMO gradient: retain ESP in the Lagrangian.')
 9042 format(4x,'FMO gradient: add ESP derivative terms.')
 9070 format(4x,'Using free molecule guess for monomers.')
 9071 format(4x,'Using monomer density to form dimer guess.') 
 9072 format(4x,'Adding HOP projector into Huckel guess.')
 9073 format(4x,'Applying HOP projector to the dimer initial orbitals.')
 9074 format(4x,'Run RHF to get DFT or MCSCF initial dimer orbitals.')
 9075 format(4x,'Compute MCSCF canonical orbitals based on the additive' 
     *      /4x,'monomer densities and use such orbitals for dimers.') 
 9077 format(4x,'Orbitals will not be reused during geometry ',
     *          'optimisation.')
 9079 format(4x,'MCSCF orbitals will be manually reordered.')
 9080 format(4x,'Use orbitals from lower levels as initial guess.')
 9085 format(4x,'BSSE correction: many body counter poise, ',
     *          'not self-consistent.')
 9086 format(4x,'BSSE correction: many body ES counter poise, ',
     *          'not self-consistent.')
 9087 format(4x,'BSSE correction: counter poise in vaccuum, ',
     *          'not self-consistent.')
 9090 format(4x,'Adding exchange terms to ESPs.')
 9091 format(4x,'No 2nd order ESP energy terms.')
 9095 format(4x,'No gradient contribution from separated dimers.') 
c9097 format(4x,'Charge transfer calculation from ifg=',I5,
c    *                                       ' to jfg=',I5)
 9096 format(4x,'The uniform distances will be used in ESPs.')
 9097 format(4x,'The consistent distances will be used in ESPs for ',
     *          'unconnected n-mers.')
 9098 format(4x,'The consistent distances will be used in ESPs for ',
     *          'all n-mers.')
 9100 format(4x,'Density will be exchanged between fragments.')
 9101 format(4x,'Orbitals will be exchanged between fragments.')
 9102 format(4x,'Energy will be exchanged between fragments.')
 9110 format(4x,'Smart parallel load balancing: large jobs first.')
 9111 format(4x,'ESPs will be shell parallelised.')
 9112 format(4x,'ESPs will be fragment parallelised.')
 9113 format(4x,'Broadcast all fragments at once.')
 9114 format(4x,'Do the fragment initialisation in the reverse order.')
 9115 format(4x,'Do not parallelise projector operator contributions.')
 9118 format(4x,'Initial density file will be broadcast from grand ',
     *          'master.')
 9119 format(4x,'Rearrange separated dimers.')
 9120 format(4x,'Restart layer: ',I2,' step: ',I2,' record ',I9)
 9121 format(4x,'Gradient restart information will be dumped to F38.') 
 9130 format(4x,'Layer',I2,', monomer(SCF,MP2,MCSCF) groups are:',3I5)
 9132 format(4x,'Layer',I2,',   dimer(SCF,ESD,MCSCF) groups are:',3I5)
 9134 format(4x,'Layer',I2,',  trimer(SCF) groups are:          ',I5)
 9140 format(/1x,'Layer',I2,' static LB for',I6,A8,'(s) larger than ',
     *           I5,' on',I5,' group(s).',/)
 9150 format(4x,'G94-like interfragment distances: rounded off to .1A.')
 9170 format(4x,'Correlation will be applied only to selected n-mers.')
 9171 format(4x,'Only selected n-mers will be computed (n>1).') 
 9172 format(4x,'Intramolecular model (within selected).',
     *      /4x,'The following monomers are selected: ',99I4)
 9173 format(4x,'Intermolecular model (between selected and the rest).',
     *      /4x,'Total binding energies will be printed for monomers:',
     *          99I4) 
 9174 format(4x,'Economy mode: no trimer arrays.') 
 9200 format(4x,'Voluminous print-out level.')
 9202 format(4x,'Compact print-out level.')
 9206 format(4x,'Fragment connexion information will be printed.')
 9208 format(4x,'Atomic charges will be printed.')
 9210 format(4x,'The total electron density in AO basis will be dumped.'
     *      )
 9216 format(4x,'The total electron density will be punched on the ',
     *          'grid:',3I6)
 9220 format(4x,'Sparse grid takes',F6.1,'% of the full grid: (',I10,
     *          ' words).')
 9222 format(4x,'Monomer grid will have dimensions of',3I6)
 9230 format(4x,'Grid padding parameter in vdW radii is',F9.3)
 9240 format(4x,'Pair interaction energy decomposition analysis (PIEDA)'
     *         ,' will be done (full).')
 9250 format(4x,'Pair interaction energy decomposition analysis (PIEDA)'
     *         ,' will be done (PL only).')
 9300 format(4x,'Number of atoms per fragment',I3)
 9500 format(/1x,'Array dimensions maxbnd,maxknd,maxcbs,maxcao are:',
     *       /17x,5I7)
 9900 format(/1x,50(1H-)//)
      return
      END
c
C*MODULE fmoio   *DECK fmoprop
      SUBROUTINE fmoprop(nder,ichfmo,nefmo,l0fmo,l1fmo,nfg2d,nfg2,nfg3,
     *                   nenm,nend,nent,needdm,ichfg,frgnam,layfrg,indat
     *                  ,iabdfg,jabdfg,itrlay,fmozan,scffrg,fmoscf,
     *                   lcorrel,modmol,exclmol,molfrg,emolfrg,emon,edim
     *                  ,etrim,edimq,esolv,emocdr,fmode,fmoq,primul,
     *                   enucfmo,atmulq,nfmoelm,fmoelm,nelm,modef,correl
     *                  ,nedimes,nedimex,ifgfmo0,doddcor,isgddi0,extracc
     *                  ,nextracc,didcc,etotdft,nedft,ctdim,n0bda,i0bda,
     *                   e0bda,efmo0,epl0ds,eint0,dopleda,rappri,nappri,
     *                   erapp,nrapp,ext2lay,DEFTF,TORQF,dolat,nunint,
     *                   nsymeq,dotd,iexcit,eexcit,texcit,nstmono,osmd,
     *                   isumd,eexfg)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      logical primul,masout,needdm,exclmol,usedij,correl,doddcor,GOPARR,
     *        DSKWRK,MASWRK,isgddi,parout,INITGDDI,isgddi0,didcc,skipc,
     *        dopleda,ext2lay,fmoq,dotd,dolat,SG1
      character*1 corri,corrs(4)
      character*3 corstr2
      character*6 corstr(3)
      character*8 charger
      parameter(MXATM=2000,MXFRG=1050,MXDFG=5,MXDPPT=MXFRG*MXDFG*12,
     *          mxrt=100,MXPTPT=100,MXGRID=10,sawdust=1.0D-25,
     *          TOKCAL=627.51D+00,one=1.0D+00,zero=0.0D+00)
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /ELPROP/ ELDLOC,ELMLOC,ELPLOC,ELFLOC,
     *                IEDEN,IEMOM,IEPOT,IEFLD,MODENS,
     *                IEDOUT,IEMOUT,IEPOUT,IEFOUT,
     *                IEDINT,IEMINT,IEPINT,IEFINT
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT,INITGDDI
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMOPT/ RABI,RASC,REFPOL,THRSLS,DENSLS,WB,WA,ETA2,GD,EVAC,
     *                RHOW,PM,AREATL,AREAKP,BONDRY,OMEGA,RET,FRO,EPSINF,
     *                EPS,DR,RSOLV,VMOL,TCE,STEN,DSTEN,CMF,TABS,IDIRCT,
     *                IPCDER,IDP,ICOMP,IFIELD,ICAV,IDISP,IPRINT,IRETCAV,
     *                ICENT,IFAST,NEVAL,IEFPOL,KEEPSM,IMGABI,IMGASC,NADD
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,nfmopcm
      COMMON /PCMPRT/ GCAVP,GCAVS,GDISP,GRP,EHFGAS
      COMMON /POINTS/ NPOINT,IPUNIT,XPOINT(MXPTPT),YPOINT(MXPTPT),
     *                ZPOINT(MXPTPT)
      COMMON /STNBUF/ STNPNT(4,2*MXATM),BIGEXP,NPTSTN,NBUFFM
      COMMON /XYZPRP/ XP,YP,ZP,
     *                DMX,DMY,DMZ,
     *                QXX,QYY,QZZ,QXY,QXZ,QYZ,
     *                QMXX,QMYY,QMZZ,QMXY,QMXZ,QMYZ,
     *                OXXX,OXXY,OXXZ,OXYY,OYYY,OYYZ,
     *                OXZZ,OYZZ,OZZZ,OXYZ,
     *                OMXXX,OMXXY,OMXXZ,OMXYY,OMYYY,
     *                OMYYZ,OMXZZ,OMYZZ,OMZZZ,OMXYZ
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      common /fmorun/ espscf,e0scf(2),IDAFMO,icurlay,icurfg,jcurfg,
     *                kcurfg,ncursh,ngau,icurpop,ifmostp,moncor,nat1e,
     *                needr,modrst,norbproj,nunesp,iskipesp,icurunt
      dimension ichfg(*),frgnam(*),layfrg(*),indat(*),iabdfg(*),
     *          jabdfg(*),itrlay(*),fmozan(*),scffrg(*),fmoscf(*),
     *          lcorrel(*),molfrg(*),emolfrg(nfg,3),emon(nfg,4,*),
     *          edim(nfg2d,*),etrim(nfg3,*),edimq(*),esolv(*),
     *          emocdr(3,*),fmode(3,natfmo,*),atmulq(natfmo,*),
     *         fmoelm(nfmoelm,*),nelm(0:*),dumbuf(3),extracc(nextracc,*)
     *         ,etotdft(*),ctdim(*),ctmax(4),ictmax(4),e0bda(4,nbdfg,*),
     *          efmo0(nfg,2,*),epl0ds(nfg,3),eint0(4),rappri(3),
     *          erapp(nappri,2),nrapp(nappri,2),DEFTF(3,NFRG,*),
     *          TORQF(3,NFRG,*),nsymeq(*),iexcit(5),eexcit(mxrt,2),
     *          texcit(3,mxrt,2),osmd(mxrt,2),isumd(mxrt+1),
     *          eexfg(2,nfg,mxrt)
      data corrs/'C','N','S','M'/
      DATA RMC/8HMCSCF   /,check/8HCHECK   /
c
c     Compute properties for FMO.
c     parstat: GlobalNone 
c     Group masters do redundant (each group does the same) calculation
c     of properties and broadcast the results onto slaves.
c
c     ifmostp=6
      if(maswrk) then 
      call timit(1)
      write(iw,9100)
      if(isgddi0) then 
        masout=iand(nprfmo,3).le.1.and.meglob.eq.0
      else
        masout=iand(nprfmo,3).le.1
      endif
      usedij=iand(IXESP,2).eq.0
      corstr(1)='uncorr' 
      if(correl) then
        if(dotd) then
          corstr(1)='exc   '
          corstr(2)='ground'
          corstr(3)='excit '
        else
          corstr(1)='corr  '
          corstr(2)='uncorr'
          corstr(3)='delta'
        endif
      endif 
      ndiv2=0 
      ndiv3=0
c
      do ilay=1,nlayer
        iter=itrlay(ilay)
        if(iter.le.mxitfg) then
          if(iter.eq.0) then
            write(iw,9106) ilay
          else
            write(iw,9109) ilay,iter,mxitfg,convfg
          endif
        else
          write(iw,9108) ilay,mxitfg,convfg
        endif
      enddo
      corstr2=' un'
      if(doddcor) corstr2='   '
      write(iw,9220) enucfmo,natfmo,nefmo,ichfmo,l1fmo,l0fmo
      write(iw,9110) 
      if(nfmoelm.gt.0) write(iw,9111) 
     *                   XPOINT(npoint),YPOINT(npoint),ZPOINT(npoint)
      if(correl) write(iw,9112) corstr2
      if(didcc) write(iw,9113)
      ENUCR=enucfmo
c
      if(nbody.lt.1.or.dotd.and.iexcit(3).eq.0) goto 1000 
c
      if(masout) then
      write(iw,9114)
c       if(ipieda.ne.0) then
        if(ipieda.gt.1) then
          if(dopleda) then
            write(iw,8117)
          else
            if(nenm.gt.1) write(iw,8115)
            if(nenm.le.1.and.iemom.ge.1) write(iw,8116)
          endif
        else
          if(nenm.gt.1) write(iw,9115)
          if(nenm.le.1.and.iemom.ge.1) write(iw,9116)
        endif
      endif
      emontot=zero
      emontotu=zero
      e0tot=zero
      e0totu=zero
      call vclr(fmoelm(1,nfg+3),1,nfmoelm) 
      epl0s=zero
      epl0d=zero
      epld=zero
      epl0di=zero
      epld_di=zero
      epld_pld=zero
      do ifg=1,nfg
        if(masout) then
        do ilay=1,layfrg(ifg)
          nelm1=nelm(1)
          if(iemom.lt.1.or.ilay.lt.layfrg(ifg)) nelm1=0
c         if(ipieda.ne.0) then
          if(ipieda.gt.1) then
            eipl0d=epl0ds(ifg,1)
            eipl0s=epl0ds(ifg,2)
            eipl0=eipl0d+eipl0s
            if(ilay.eq.layfrg(ifg)) e0tot=e0tot+efmo0(ifg,1,ilay)
            if(correl) then
              eipld=emon(ifg,2,ilay)-efmo0(ifg,2,ilay)
              eipl0di=epl0ds(ifg,3)
              eipld_di=emon(ifg,1,ilay)-efmo0(ifg,1,ilay)-eipld-eipl0di
              if(ilay.eq.layfrg(ifg)) e0totu=e0totu+efmo0(ifg,2,ilay)
            else
              eipld=emon(ifg,1,ilay)-efmo0(ifg,1,ilay)
              eipl0di=zero
              eipld_di=zero
            endif
            eipld_pld=eipld-eipl0d
            epl0s=epl0s+eipl0s
            epl0d=epl0d+eipl0d
            epl0di=epl0di+eipl0di
            epld_di=epld_di+eipld_di
            epld_pld=epld_pld+eipld_pld 
            epld=epld+eipld
            if(dopleda) then
              write(iw,8118) ifg,frgnam(ifg),ilay,eipl0d*tokcal,
     *                       eipl0s*tokcal,eipl0*tokcal,eipl0di*tokcal,
     *                       (fmoelm(ielm,ifg),ielm=nelm(0)+1,nelm1)
            else
            if(nenm.gt.1) then
              write(iw,9117) ifg,frgnam(ifg),ilay,eipl0d*tokcal,
     *                       eipl0s*tokcal,eipl0*tokcal,eipl0di*tokcal,
     *                       eipld_pld*tokcal,eipld_di*tokcal,
     *                       (fmoelm(ielm,ifg),ielm=nelm(0)+1,nelm1)
            else
              write(iw,9119) ifg,frgnam(ifg),ilay,eipl0d*tokcal,
     *                      eipl0s*tokcal,eipl0*tokcal,eipld_pld*tokcal,
     *                       (fmoelm(ielm,ifg),ielm=nelm(0)+1,nelm1)
            endif
            endif
          else 
            if(nenm.gt.1) then
              write(iw,9118) ifg,frgnam(ifg),ilay,(emon(ifg,k,ilay),k=1,
     *                     nenm),(fmoelm(ielm,ifg),ielm=nelm(0)+1,nelm1)
            else
              write(iw,9120) ifg,frgnam(ifg),ilay,(emon(ifg,k,ilay),k=1,
     *                     nenm),(fmoelm(ielm,ifg),ielm=nelm(0)+1,nelm1)
            endif
          endif 
        enddo
        if(iemom.ge.2) 
     *    write(iw,8122) (fmoelm(ielm,ifg),ielm=nelm(1)+1,nelm(2))
        if(iemom.ge.3) 
     *    write(iw,8123) (fmoelm(ielm,ifg),ielm=nelm(2)+1,nelm(3))
        endif
        ilay=layfrg(ifg)
        emontot=emontot+emon(ifg,1,ilay)
        if(nenm.ge.2) emontotu=emontotu+emon(ifg,2,ilay)
        call daxpy(nfmoelm,one,fmoelm(1,ifg),1,fmoelm(1,nfg+3),1)
      enddo
c
      IF (IEFPFMO.EQ.1) THEN
        emontot=emontot+REPNUCEFP
        emontotu=emontotu+REPNUCEFP
      END IF
      etot=emontot
      etotu=emontotu
c
      if(ifgfmo0.ne.0) then
c       Save some data for FMO0, based on just fragment 1.
        ETOT=emon(1,1,layfrg(1))
        ESCF=emon(1,2,layfrg(1))
c       This is a despicable way to store ES moments. 
        if(nfmoelm.ne.0) call dcopy(nfmoelm,fmoelm(1,1),1,DMX,1)
      endif
      if(ipieda.eq.2) then
        write(iw,9200) corstr(1),0,e0tot
        if(correl) write(iw,9200) 'uncorr',0,e0totu
        if(correl) write(iw,9200) 'delta ',0,e0tot-e0totu
      endif
      write(iw,9200) corstr(1),1,emontot
c     for rcorsd.eq.0 the data are not accumulated for uncorr. 
      if(nenm.ge.2.and.rcorsd.ne.0) then
        write(iw,9200) corstr(2),1,emontotu
        write(iw,9200) corstr(3),1,emontot-emontotu
        if(didcc) call fmoccp(1,emontotu,extracc(1,1))
      endif
c     if(ifgdon.ne.0.and.ifgacc.ne.0) then
c        ilayd=layfrg(ifgdon)
c        ilayc=layfrg(ifgacc)
c        edonacc=emon(ifgdon,1,ilayd)+emon(ifgacc,1,ilayc)
c        write(iw,9420) corstr2,edonacc
c        if(nenm.gt.1) then
c          edonacc=emon(ifgdon,2,ilayd)+emon(ifgacc,2,ilayc)
c          write(iw,9420) 'uncorr',edonacc
c        endif
c     endif
      if(iemom.ge.1)
     *  write(iw,9230) 1,(fmoelm(ielm,nfg+3),ielm=nelm(0)+1,nelm(1)),
     *                  DNRM2(3,fmoelm(1,nfg+3),1)
      if(iemom.ge.2) 
     *  write(iw,9231) (fmoelm(ielm,nfg+3),ielm=nelm(1)+1,nelm(2))
      if(iemom.ge.3) 
     *  write(iw,9232) (fmoelm(ielm,nfg+3),ielm=nelm(2)+1,nelm(3))
c
      if(nder.gt.0.and.modef.ne.1) then
        call fmogout(1,masout,indat,fmozan,fmode,DEFTF,TORQF)
      endif
C
      if(nbody.lt.2) goto 1000 
c
      if(masout) then
      write(iw,9121)
      if(nbsse.eq.0) then
        if(ipieda.ne.0) then
          if(rcorsd.eq.0) then
            if(dopleda) then 
              write(iw,8124)
            else 
              write(iw,9124)
            endif 
          else 
            if(dopleda) then
              write(iw,8125)
            else
              write(iw,9125)
            endif
          endif
        else 
          if(dolat) then
            if(rcorsd.eq.0) then
              write(iw,8120)
            else
              write(iw,8121)
            endif
          else
            if(rcorsd.eq.0) then
              write(iw,9122)
            else
              write(iw,9123)
            endif
          endif
        endif 
      else if(nbsse.eq.1.or.nbsse.eq.2) then
        write(iw,9126)
      else if(nbsse.eq.3) then
        if(nend.gt.6) write(iw,9128)
        if(nend.le.6) write(iw,9129)
      endif
      endif
      ijfg=0
      e0mon=0
      e0monu=0
      edimtot=zero
      edimtotu=zero
      edimcon=zero
      edimconu=zero
      edimijq=zero
      edimijqu=zero
      eddim=zero
      eddimu=zero
      debsse=zero
      debsseu=zero
      debsseuij=zero
      ees=zero
      eex=zero
      ect=zero
      edi=zero
      eesbda=zero
      eexbda=zero
      ectbda=zero
      edibda=zero
      ebdaes=zero
      ebdaex=zero
      ebdact=zero
      ebdadi=zero
      rij=-one
      sepmax=zero 
c     ires=0
c     -1 is supposed to make it look unknown which it is if so printed.
      if(molfrg(1).ne.-1) call vclr(emolfrg,1,nfg*3)
      do ifg=1,nfg
        ilay=layfrg(ifg)
        njfg=ifg-1
        do 500 iu=0,nunint
        icurunt=iu
        if(iu.gt.0) then
          if(nsymeq(iu).le.0) goto 500
        endif
        do jfg=1,njfg
c       do jfg=1,ifg-1
          ijfg=ijfg+1
          ijlay=min(ilay,layfrg(jfg))
c         if(molfrg(1).ge.0) ires=molfrg(ifg)+molfrg(jfg)
          ires=molfrg(ifg)+molfrg(jfg)
          if(exclmol) then
            skipc=ires.ge.0.and.ires.ne.2
          else
            skipc=ires.eq.0.or.ires.eq.2
          endif
          if(iand(modmol,1).ne.0.and.skipc) goto 100 
          if(needr.ne.0) rij=fmodist(ifg,0,0,jfg) 
c         write(6,*) 'www',needr,rij,ifg,jfg  
          sepmax=max(sepmax,rij)
          irij=int(rij*1.0D+02)
          ijcharge=ichfg(ifg)*ichfg(jfg)
          memon=1
c         memon1=0
          memon1=2
c         if(correl) memon1=1
          if(iand(lcorrel(ijlay),1).ne.0) memon1=1
          if(resdim.ne.0) then
            if(rij.gt.resdim) memon1=3
            if(rcorsd.ne.0) then
c             ires=molfrg(ifg)+molfrg(jfg)
              if(rij.gt.rcorsd.or.skipc) then
                memon=2
                memon1=max(memon1,2)
              endif
            endif 
          endif 
          if((scffrg(ifg).eq.rmc.or.scffrg(jfg).eq.rmc).and.
     *       fmoscf(ijlay).eq.rmc.and.(resdim.eq.0.or.rij.le.resdim)) 
     *      memon1=4 
          memonb=memon+2
          if(dopleda) memon=memonb
c         corri=corrs(memon)
          corri=corrs(memon1) 
        edimij1=edim(ijfg,1)-emon(ifg,memon,ijlay)-emon(jfg,memon,ijlay)
          edimij2u=edim(ijfg,2)
          edimij2=edim(ijfg,4)
          if(edim(ijfg,1).eq.zero) ndiv2=ndiv2+1
c
c         This is the latest fix to replace delta-D (MP2) by RHF.
c         to put back, uncontract the if below and remove the line next to it.
c
c         if(.not.doddcor.or.memon.eq.2) edimij2=edimij2u
          if(fmoq) edimijq=edimq(ijfg)
          edimijqu=edimijq
          edimij2=edimij2u
          eddim=eddim+edimij2
          eddimu=eddimu+edimij2u
c         if(nend.lt.3) then
          if(rcorsd.eq.0) then
            edimij1u=edimij1
          else
            memon2=2 
            if(dopleda) memon2=memon2+2
            edimij1u=edim(ijfg,3)-emon(ifg,memon2,ijlay)
     *                           -emon(jfg,memon2,ijlay)
          endif
          if(abs(rij).lt.1.0D-08) then
            edimcon=edimcon+edimij1+edimij2+edimijq
            edimconu=edimconu+edimij1u+edimij2u+edimijqu
          endif
c         if(nend.gt.4) then 
          if(nbsse.ne.0) then
            edimij3=edim(ijfg,5)
            edimij4=edim(ijfg,6)
            if(edimij3.eq.sawdust) edimij3=zero
            if(edimij4.eq.sawdust) edimij4=zero
            if(nend.le.6) then
              edimijcb=edim(ijfg,6)
            else
              edimijcb=edim(ijfg,7)
            endif
          endif
          edimes=edim(ijfg,nedimes)
          edimex=edim(ijfg,nedimex)
          if(masout) then
          if(nbsse.eq.0) then 
            detot=edimij1+edimij2+edimijq
            detotu=edimij1u+edimij2u+edimijqu
c           detots=detot
c           Remove artificial tens of hartrees for connected monomers.
c           If distances are not available, use plain energy to judge.
            if(ipieda.eq.0.and.(needr.ne.0.and.rij.eq.zero.or.
     *         needr.eq.0.and.edimij1.lt.-1.0D+01)) then
              detot=edimij2+edimijq
              detotu=edimij2u+edimijqu
            endif
            ctij=zero
            if(primul) ctij=ctdim(ijfg)
c           rij means two fragments are connected by a bond and thus
c           edimij1 contains some meaningless (well,?) value.
            if(ipieda.ne.0) then
c
              edimdi=edimij1-edimij1u
              if(rcorsd.eq.0) edimdi=zero
              edimct=detot-edimes-edimex-edimdi
              if(dopleda.and.memon1.ne.3) then 
                epl0ij=epl0ds(ifg,1)+epl0ds(ifg,2)+
     *                 epl0ds(jfg,1)+epl0ds(jfg,2)
                edimct=edimct-epl0ij
c               In PL0 runs polarisation is double counted.
              else
                epl0ij=zero
              endif
              edimij1b=edimij1
              if(n0bda.ne.0) then
                e0bda(1,i0bda,1)=edimes
                e0bda(2,i0bda,1)=edimex
                e0bda(3,i0bda,1)=edimct
                e0bda(4,i0bda,1)=edimdi
              else if(nbdfg.ne.0) then
                ebdaes=ebdaes+edimes
                ebdaex=ebdaex+edimex
                ebdact=ebdact+edimct
                ebdadi=ebdadi+edimdi
                call bdasub(ifg,jfg,correl,edimij1b,detot,detotu,edimes,
     *                      edimex,edimct,edimdi,e0mon,e0monu,indat,
     *                      iabdfg,jabdfg,e0bda(1,1,ijlay),ires)
                ebdaes=ebdaes-edimes
                ebdaex=ebdaex-edimex
                ebdact=ebdact-edimct
                ebdadi=ebdadi-edimdi
                if(ires.eq.0) then
                  eesbda=eesbda+edimes
                  eexbda=eexbda+edimex
                  ectbda=ectbda+edimct
                  edibda=edibda+edimdi
                endif
              endif
              ees=ees+edimes
              eex=eex+edimex
              ect=ect+edimct
              edi=edi+edimdi
              if(rcorsd.eq.0) then
                if(dopleda) then
                write(iw,9132) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edimij1b*tokcal,edimij2*tokcal,
     *                         detot*tokcal,epl0ij*tokcal,
     *                         edimes*tokcal,edimex*tokcal,edimct*tokcal
                else
                write(iw,9134) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edimij1b*tokcal,edimij2*tokcal,
     *                         detot*tokcal,
     *                         edimes*tokcal,edimex*tokcal,edimct*tokcal
                endif
              else
                if(dopleda) then
                write(iw,9132) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edimij1b*tokcal,edimij2*tokcal,
     *                         detot*tokcal,epl0ij*tokcal,
     *                         edimes*tokcal,edimex*tokcal,edimct*tokcal
     *                        ,edimdi*tokcal
                else
                write(iw,9134) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edimij1b*tokcal,edimij2*tokcal,
     *                         detot*tokcal,
     *                         edimes*tokcal,edimex*tokcal,edimct*tokcal
     *                        ,edimdi*tokcal
                endif
              endif 
            else  
              if(dolat) then
              if(rcorsd.eq.0) then
                write(iw,9161) ifg,jfg,iu,nsymeq(iu),corri,ijlay,
     *                         ijcharge,rij,ctij,edim(ijfg,1),edimij1,
     *                         edimij2,detot*tokcal
              else
                write(iw,9163) ifg,jfg,iu,nsymeq(iu),corri,ijlay,
     *                         ijcharge,rij,ctij,edim(ijfg,1),
     *                         edim(ijfg,3),edimij1,edimij1u,
     *                         edimij2,detot*tokcal
              endif
              else
              if(rcorsd.eq.0) then
                write(iw,9131) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edim(ijfg,1),edimij1,edimij2,
     *                         detot*tokcal
c    *                         edim(ijfg,1),edimij1,edimij2,edimijq,
              else
                write(iw,9133) ifg,jfg,corri,ijlay,ijcharge,rij,ctij,
     *                         edim(ijfg,1),edim(ijfg,3),edimij1,
     *                         edimij1u,edimij2,detot*tokcal
c    *                         edimij1u,edimij2,edimijq,detot*tokcal
              endif 
              endif 
            endif 
            debsseij=zero
          else if(nbsse.eq.1.or.nbsse.eq.2) then 
            debsseij=edimij3+edimij4-edimij1-edimij2
            if(edimij3+edimij4.eq.zero) debsseij=zero
            detot=edimij1+edimij2+debsseij
            write(iw,9140) ifg,jfg,corri,ijlay,edim(ijfg,1),edimij1,
     *            edimij1u,edimij2,edimij3,edimij4,debsseij,detot*tokcal
          else if(nbsse.eq.3) then 
            debsseij=emon(ifg,memonb,ijlay)+emon(jfg,memonb,ijlay)-
     *               edimij3
            debsseuij=emon(ifg,4,ijlay)+emon(jfg,4,ijlay)-edimijcb
c           take care of separated dimers that ignore BSSE
            if(edimij3.eq.zero) debsseij=zero
            if(edimij3.eq.zero.or.nend.lt.7) debsseuij=zero
            detot=edimij1+edimij2+debsseij
            if(nend.ge.7) then 
              write(iw,9150) ifg,jfg,corri,ijlay,irij,edim(ijfg,1),
     *          edimij1,edimij1u,edimij2,debsseij,debsseuij,detot*tokcal
            else
              write(iw,9152) ifg,jfg,corri,ijlay,irij,edim(ijfg,1),
     *           edimij1,edimij2,debsseij,detot*tokcal
            endif
          endif
          endif
          if(dopleda.and.memon1.ne.3) then
c           Replace energies by PL0 values to get the PIEDA energy.
c           memon1=3 means separated dimers, in which case
c           edimij1 is measured against FMO0, and there is no need to
c           adjust it. 
            edimij1=edim(ijfg,1)-emon(ifg,memon-2,ijlay)
     *                          -emon(jfg,memon-2,ijlay)
            edimij2=edim(ijfg,4)
            edimij1u=edimij1
            edimij2u=edimij2
            if(rcorsd.ne.0)
     *        edimij1u=edim(ijfg,3)-emon(ifg,2,ijlay)-emon(jfg,2,ijlay)
          endif
          edimtot=edimtot+edimij1+edimij2+edimijq
          edimtotu=edimtotu+edimij1u+edimij2u+edimijqu
c         write(6,*) 'wwwEE',ifg,jfg,edimij1,edimij1u,edimij2
          if(nappri.ne.0) then 
            if(rij.le.rappri(1)) then
              iappri=1
            else
              iappri=int((rij-rappri(1))/rappri(3))+2
            endif
            if(iappri.le.nappri) then
              erapp(iappri,1)=erapp(iappri,1)+edimij1+edimij2+edimijq
              erapp(iappri,2)=erapp(iappri,2)+edimij1u+edimij2u+edimijqu
              nrapp(iappri,1)=nrapp(iappri,1)+1 
              if(memon1.eq.1) nrapp(iappri,2)=nrapp(iappri,2)+1
            endif
          endif
          debsse=debsse+debsseij
          debsseu=debsseu+debsseuij
c         if(ires.eq.0.or.ires.eq.2) then
          if(.not.exclmol.and.ires.eq.1) then
            if(molfrg(ifg).gt.0) then
              emolfrg(ifg,1)=emolfrg(ifg,1)+detot
              emolfrg(ifg,2)=emolfrg(ifg,2)+detotu
              if(ipieda.ne.0) emolfrg(ifg,3)=emolfrg(ifg,3)+edimes
            endif
            if(molfrg(jfg).gt.0) then
              emolfrg(jfg,1)=emolfrg(jfg,1)+detot
              emolfrg(jfg,2)=emolfrg(jfg,2)+detotu
              if(ipieda.ne.0) emolfrg(jfg,3)=emolfrg(jfg,3)+edimes
            endif
          endif
 100    continue
        enddo
        njfg=nfg
  500 continue
      enddo
      if(molfrg(1).ge.0.and..not.exclmol) then
        write(iw,9160) 
        do ifg=1,nfg
          if(molfrg(ifg).gt.0) then
            if(ipieda.ne.0) then
              write(iw,9175) ifg,frgnam(ifg),(emolfrg(ifg,i)*tokcal,
     *                       i=3,1,-1)
            else
              write(iw,9170) ifg,frgnam(ifg),(emolfrg(ifg,i)*tokcal,
     *                       i=2,1,-1)
            endif
          endif
        enddo
c       No total 2-body properties are available.
        goto 1000 
      endif
c     if(nbsse.eq.0) then
c       etot=edimtot-(nfg-2)*emontot
c     else
      etot=emontot+edimtot
      etotu=emontotu+edimtotu
c     At this point fmoelm(1,nfg+3) contains the 1st order moments
c                   fmoelm(1,nfg+2) contains the 3rd order correction 
c                   fmoelm(1,nfg+1) contains the 2nd order correction 
      call daxpy(nfmoelm,one,fmoelm(1,nfg+3),1,fmoelm(1,nfg+1),1)
c     endif
c     enucr=ENUC(natfmo,fmoZan,fmoC)
c
c     for BSSE runs, the total energy (as will be used during optimisation)
c     includes BSSE correction, but the gradient doesn't!
c 
      if(ndiv2.ne.0.and.runtyp.ne.check) write(iw,9500) ndiv2,2
      if(ipieda.ne.0.and.nbdfg.ne.0.and.n0bda.eq.0.and.e0mon.ne.0) then
        write(iw,9201) corstr(1),1,emontot+e0mon
        if(nenm.ge.2.and.rcorsd.ne.0) then
          write(iw,9201) 'uncorr',1,emontotu+e0monu
          write(iw,9201) 'delta ',1,emontot-emontotu+(e0mon-e0monu)
        endif
      endif
      if(nbsse.eq.0) then
        write(iw,9200) corstr(1),2,etot
c       if(doddcor) write(iw,9200) 'corr-D',2,etot-eddim+eddimu
c       corr-D has no value because to get it one needs E'I with RHF delta-D. 
        if(nenm.gt.1) then
          write(iw,9200) corstr(2),2,etotu
          write(iw,9200) corstr(3),2,etot-etotu
          if(didcc) call fmoccp(2,etotu,extracc(1,2))
        endif
      else
        write(iw,9210) corstr(1),etot,etot+debsse
        etot=etot+debsse
        if(nend.ge.7) then
          write(iw,9210) 'uncorr',etotu,etotu+debsseu 
          escf=escf+debsseu
        endif
      endif
c     EBB may be defined without distances known
c     but then the code should get the connecting
c     information from elsewhere.
      if(nbdfg.ne.0.and.needr.ne.0) then
        write(iw,9202) corstr(1),2,emontot+edimcon
        if(correl) write(iw,9202) 'uncorr',2,emontotu+edimconu
      endif
      if(nappri.ne.0) then
        do i=2,nappri
          erapp(i,1)=erapp(i,1)+erapp(i-1,1)
          erapp(i,2)=erapp(i,2)+erapp(i-1,2)
          nrapp(i,1)=nrapp(i,1)+nrapp(i-1,1)
          nrapp(i,2)=nrapp(i,2)+nrapp(i-1,2)
        enddo
        write(iw,*) ' ' 
        do i=1,nappri
          riapp=rappri(1)+rappri(3)*(i-1)
c         sepmax=min(sepmax,resdim)
          if(riapp.gt.sepmax.and.riapp-rappri(3).le.sepmax) riapp=sepmax
          if(riapp.le.sepmax.and.sepmax.ne.0) then 
          if(correl) then
           write(iw,9205) riapp,nrapp(i,2),nrapp(i,1),2,
     *       erapp(i,1)+emontot-erapp(i,2)-emontotu,erapp(i,2)+emontotu
          else
           write(iw,9206) riapp,nrapp(i,1),2,erapp(i,1)+emontot
          endif
          endif
        enddo
        write(iw,*) ' ' 
      endif
c
      if(ipieda.ne.0) then
        if(dopleda) then
          eint0(1)=ees
          eint0(2)=eex
          eint0(3)=ect
          eint0(4)=edi
          epl0=epl0d+epl0s
          eib0=epl0+epl0di+ees+eex+ect+edi
          e1int=eesbda+eexbda+ectbda+edibda
          write(iw,9203) corstr(1),2,e0tot+e0mon+eib0
          if(correl) write(iw,9203) 'uncorr',2,e0totu+e0monu+
     *                                         eib0-epl0di-edi
          if(correl) write(iw,9203) 'delta ',2,e0tot+e0mon-e0totu-e0monu
     *                                        +epl0di+edi
          if(ipieda.eq.2) then
            write(iw,9600) epl0d*tokcal,
     *                     epl0s*tokcal,
     *                         epl0*tokcal,
     *                         epl0di*tokcal,
     *                         ees*tokcal,
     *                         eex*tokcal,
     *                         ect*tokcal,
     *                         edi*tokcal,
     *                             eib0*tokcal,(eib0-epl0di-edi)*tokcal
            if(nbdfg.ne.0) then
             write(iw,9620) eesbda*tokcal,eexbda*tokcal,
     *                     ectbda*tokcal,edibda*tokcal,e1int*tokcal,
     *                     ebdaes*tokcal,ebdaex*tokcal,ebdact*tokcal,
     *                     ebdadi*tokcal,e0mon*tokcal
            endif
          endif
        else
          ees0=eint0(1)
          eex0=eint0(2)
          ect0=eint0(3)
          edi0=eint0(4)
          epls_pls=ees-ees0-epl0s
          epl_pl=epld_pld+epls_pls
          epls_ex=eex-eex0
          epls_ct=ect+epl_pl-ect0
          epls_di=edi-edi0
          epl0=epl0d+epl0s
          epl_di=epld_di+epls_di
          eib0=epl0+epl0di+ees0+eex0+ect0+edi0
          eib=eib0+epls_ex+epls_ct+epl_di
          eint=ees+eex+ect+edi
          e1int=eesbda+eexbda+ectbda+edibda
          epl=epl0+epl_pl
          epldcorr=epl0d+epld_pld+epl0di+epld_di
          if(ipieda.eq.2) then
          write(iw,9203) corstr(1),2,e0tot+e0mon+eib0
          if(correl) write(iw,9203) 'uncorr',2,e0totu+e0monu+
     *                                         eib0-epl0di-edi0
          if(correl) write(iw,9203) 'delta ',2,e0tot+e0mon-e0totu-e0monu
     *                                        +epl0di+edi0
            write(iw,9600) epl0d*tokcal,
     *                     epl0s*tokcal,
     *                         epl0*tokcal,
     *                         epl0di*tokcal,
     *                         ees0*tokcal,
     *                         eex0*tokcal,
     *                         ect0*tokcal,
     *                         edi0*tokcal,
     *                             eib0*tokcal,(eib0-epl0di-edi0)*tokcal
     *                            ,epls_ex*tokcal,
     *                             epls_ct*tokcal,
     *                         epld_di*tokcal,
     *                         epls_di*tokcal,
     *                             epl_di*tokcal,
     *                                 eib*tokcal,
     *                                 (eib-epl0di-edi0-epl_di)*tokcal
c
            write(iw,9605) epld_pld*tokcal,epls_pls*tokcal,
     *                     epl_pl*tokcal,epl*tokcal,(ect0-epl_pl)*tokcal
     *                    ,epldcorr*tokcal
          endif
c
          write(iw,9610) ees*tokcal,eex*tokcal,ect*tokcal,edi*tokcal,
     *                   eint*tokcal
          if(nbdfg.ne.0) then
             write(iw,9620) eesbda*tokcal,eexbda*tokcal,
     *                     ectbda*tokcal,edibda*tokcal,e1int*tokcal,
     *                     ebdaes*tokcal,ebdaex*tokcal,ebdact*tokcal,
     *                     ebdadi*tokcal,e0mon*tokcal
             endif
        endif
      endif
c
c     write(iw,9229+i) (fmoelm(ielm,nfg+1),ielm=nelm(i-1)+1,nelm(i))
      if(iemom.ge.1) 
     *  write(iw,9230) 2,(fmoelm(ielm,nfg+1),ielm=nelm(0)+1,nelm(1)),
     *                  DNRM2(3,fmoelm(1,nfg+1),1)
      if(iemom.ge.2) 
     *  write(iw,9231) (fmoelm(ielm,nfg+1),ielm=nelm(1)+1,nelm(2))
      if(iemom.ge.3) 
     *  write(iw,9232) (fmoelm(ielm,nfg+1),ielm=nelm(2)+1,nelm(3))
c
      if(primul) then
        write(iw,9420)
        nctsiz=min(nfg-1,4)
        ctmall=zero
        ctdall=zero
        do ifg=1,nfg
          cttot=zero
          call vclr(ctmax,1,nctsiz)
          call viclr(ictmax,1,nctsiz)
          do jfg=1,nfg
            if(ifg.ne.jfg) then
              if(ifg.gt.jfg) then
                ctij=-ctdim(((ifg-1)*(ifg-2))/2+jfg)
                ctdall=ctdall+abs(ctij)
              else 
                ctij= ctdim(((jfg-1)*(jfg-2))/2+ifg)
              endif 
              cttot=cttot+ctij
              imin=1
              do ii=2,nctsiz
                if(abs(ctmax(ii)).lt.abs(ctmax(imin))) imin=ii
              enddo
c             Replace the smallest element in ctmax by ctij
              if(abs(ctmax(imin)).lt.abs(ctij)) then
                ctmax(imin)=ctij
                ictmax(imin)=jfg
              endif
            endif
          enddo
          call dasort(nctsiz,ctmax,ictmax)
          write(iw,9425) ifg,ichfg(ifg),cttot,
     *                   (ictmax(ii),ctmax(ii),ii=1,nctsiz)
          ctmall=ctmall+abs(cttot)
        enddo
        write(iw,9430) ctmall,ctdall
      endif
c
      if(nder.gt.0.and.modef.ne.1) then
        call fmogout(2,masout,indat,fmozan,fmode,DEFTF,TORQF)
      endif
c
      if(nfg3.eq.0) goto 1000 
c     At present FMO3 does not work with BSSE
      if(masout) then
      write(iw,9250)
      if(needdm) then
        write(iw,9260)
      else
        write(iw,9265)
      endif
      endif 
      etrimtot1=zero
      etrimtot2=zero
      etrimtot2a=zero
      etrimtot1u=zero
      ijkfg=0
      rmin=-one
      rmax=-one
      rrij=zero
      rrik=zero
      rrjk=zero
      if(nappri.ne.0) call vclr(erapp,1,nappri*2)
      if(nappri.ne.0) call viclr(nrapp,1,nappri*2)
      sepmax=zero
      do ifg=1,nfg
        ilay=layfrg(ifg)
        do jfg=1,ifg-1
          ijlay=min(ilay,layfrg(jfg))
          do kfg=1,jfg-1
            if(needr.ne.0) call fmodist3(ifg,jfg,kfg,rmin,rmax)
            sepmax=max(sepmax,rmax)
            ijkfg=ijkfg+1
            ijklay=min(ijlay,layfrg(kfg))
            memon=1
            metrim=1
            if(correl) then
              ires=molfrg(ifg)+molfrg(jfg)+molfrg(kfg)
              if(exclmol) then
                skipc=ires.ge.0.and.ires.ne.3
              else
                skipc=ires.ge.0.and.ires.ne.1
              endif
              if(rmax.gt.restri(4).and.restri(4).ne.0 .or. skipc) then
                 memon=2
                 metrim=nent
              endif
            endif
            if(ext2lay) then
              memon=2
              metrim=1
            endif
            memon1=memon
            if(iand(lcorrel(ijklay),1).eq.0) memon1=2
c           if(.not.correl) memon1=2
            corri=corrs(memon1)
            emonijk=emon(ifg,memon,ijklay)+emon(jfg,memon,ijklay)+
     *              emon(kfg,memon,ijklay)
            indij=(ifg*ifg-3*ifg)/2+jfg+1
            indik=(ifg*ifg-3*ifg)/2+kfg+1
            indjk=(jfg*jfg-3*jfg)/2+kfg+1
            eijk=etrim(ijkfg,metrim)
c
c           For trimers it is more difficult to identify divergence. 
            if(etrim(ijkfg,1).eq.one) then
              eijk=one
              ndiv3=ndiv3+1
            endif
c
            ddijk=etrim(ijkfg,2)
c           if(ddijk.ne.0.or..not.usedij) then
            if(eijk.ne.0.or..not.usedij) then
              if(eijk.eq.one) eijk=zero
c             deijk=eijk+emonijk-edim(indij,medimij)-edim(indik,medimik)-
c    *                           edim(indjk,medimjk)
c             Actually one only has to set memij,memik,memjk properly.
c             It makes no difference with dimers (medij,medik,medjk) because 
c             for them uncorrelated values are copied to the correlated location
c             so if one accesses the wrong place by mistake there is no 
c             problem. But monomers always (well, almost) have different 
c             corr/uncorr values so proper indices memij,memik,memjk must be 
c             computed.
              medij=1
              medik=1
              medjk=1
              memij=1
              memik=1
              memjk=1
c             The code below has not yet been changed to accomodate exclmol.
              if(correl) then
                if(resdim.ne.0) then
                  rrij=fmodist(ifg,0,0,jfg)
                  rrik=fmodist(ifg,0,0,kfg)
                  rrjk=fmodist(jfg,0,0,kfg)
                endif
                ires=molfrg(ifg)+molfrg(jfg)
                if(memon.eq.2.or.(rrij.gt.rcorsd.and.rcorsd.ne.0).or.
     *             ires.eq.0.or.ires.eq.2) then
                  medij=3
                  memij=2
                endif
                ires=molfrg(ifg)+molfrg(kfg)
                if(memon.eq.2.or.(rrik.gt.rcorsd.and.rcorsd.ne.0).or.
     *             ires.eq.0.or.ires.eq.2) then
                  medik=3
                  memik=2
                endif
                ires=molfrg(jfg)+molfrg(kfg)
                if(memon.eq.2.or.(rrjk.gt.rcorsd.and.rcorsd.ne.0).or.
     *             ires.eq.0.or.ires.eq.2) then
                  medjk=3
                  memjk=2
                endif
              endif
              deij=edim(indij,medij)-
     *               emon(ifg,memij,ijklay)-emon(jfg,memij,ijklay)
              deik=edim(indik,medik)-
     *               emon(ifg,memik,ijklay)-emon(kfg,memik,ijklay)
              dejk=edim(indjk,medjk)-
     *               emon(jfg,memjk,ijklay)-emon(kfg,memjk,ijklay)
              deijk=eijk-emonijk-deij-deik-dejk
c
              ddijk=ddijk-edim(indij,2)-edim(indik,2)-edim(indjk,2)
              etrimtot1=etrimtot1+deijk
              if(rcorsd.ne.0) then
                deijku=etrim(ijkfg,nent)+
     *         emon(ifg,2,ijklay)+emon(jfg,2,ijklay)+emon(kfg,2,ijklay)-
     *          edim(indij,3)-edim(indik,3)-edim(indjk,3)
                etrimtot1u=etrimtot1u+deijku
              else
                deijku=deijk
              endif
              if(masout) then
              if(needdm) then
c               Use V-consistent ddDVijk Tr(delta-delta-Density-ijk * V-ijk).
c               However, the 3-body effect that is printed as etot is computed
c               from the pure 3-body quantity ddVijk. This might be confusing,
c               but the other possibility dddijk might be equally so.
c               In the ddDVijk formula one should ignore dDij and use only dddijk;
c               but then the notation of purely three-body density effect is lost.
                dddijk=etrim(ijkfg,3) 
                detot=deijk+ddijk
                etrimtot2=etrimtot2+dddijk
                etrimtot2a=etrimtot2a+ddijk
                write(iw,9275) ifg,jfg,kfg,corri,ijklay,rmin,rmax,eijk
     *                        ,deijk,deijku,ddijk,dddijk,detot*tokcal
              else
                detot=deijk+ddijk
                etrimtot2=etrimtot2+ddijk
                write(iw,9280) ifg,jfg,kfg,corri,ijklay,rmin,rmax,eijk
     *                        ,deijk,deijku,ddijk,detot*tokcal
              endif
              endif
              if(nappri.ne.0) then
                if(rmax.le.rappri(1)) then
                  iappri=1
                else
                  iappri=int((rmax-rappri(1))/rappri(3))+2
                endif
                if(iappri.le.nappri) then
                  erapp(iappri,1)=erapp(iappri,1)+deijk+ddijk
                  erapp(iappri,2)=erapp(iappri,2)+deijku+ddijk
                  nrapp(iappri,1)=nrapp(iappri,1)+1
                  if(memon1.eq.1) nrapp(iappri,2)=nrapp(iappri,2)+1
                endif
              endif
            endif
          enddo
        enddo
      enddo
c
      etota=etot+etrimtot1+etrimtot2a
      etotau=etotu+etrimtot1u+etrimtot2a
c     if(needr2.ne.0) then
c     As mentioned above, the ddDVijk equation should have two-body dDij
c     subtracted as it was added in the two-body sum.
      if(needdm) then
        etot=etot-eddim
        etotu=etotu-eddimu
      endif
      etotsav=etot
      etotusav=etotu
      etrimtot=etrimtot1+etrimtot2
      etrimtotu=etrimtot1u+etrimtot2
      etot1=etot+etrimtot1
      etot1u=etotu+etrimtot1u
      etot=etot+etrimtot
      etotu=etotu+etrimtotu
      if(ndiv3.ne.0.and.runtyp.ne.check) write(iw,9500) ndiv3,3
      write(iw,9200) corstr(1),3,etot
      if(needdm) then
c       needdm prints three energies:
c       Euncorr is the EI'+dE'IJ+dE'IJK+ddDijkVijk
c       E no D  is the EI'+dE'IJ+dE'IJK
c       E D23   is the EI'+dE'IJ+dE'IJK+dDijVij
c                         +(dDijkVijk-dDijVij-dDikVik-dDjkVjk)
c               where dDijk=Dijk-(Di+Dj+Dk).
c       "E no D" is useless, Euncorr is proper only for no ritrim approximation.
c       "Euncorr" is equal to "E no D" without approximations.
c       "E no D" is the most sound one at present. It is used as Euncorr for 
c       needdm=.false. runs (production).
c       Its drawback is a fairly high sensitivity to ES approximations,
c       as it is based on Vij and Vijk.
c
        write(iw,9200) ' no D ',3,etot1
        write(iw,9200) ' D23  ',3,etota
      endif
      if(rcorsd.ne.0) then
        write(iw,9200) 'uncorr',3,etotu
        write(iw,9200) 'delta ',3,etot-etotu
        if(needdm) then
          write(iw,9200) 'uncnoD',3,etot1u
          write(iw,9200) 'uncD23',3,etotau
        endif
        if(didcc) call fmoccp(3,etotu,extracc(1,3))
      endif
c
      if(nappri.ne.0) then
        do i=2,nappri 
          erapp(i,1)=erapp(i,1)+erapp(i-1,1)
          erapp(i,2)=erapp(i,2)+erapp(i-1,2)
          nrapp(i,1)=nrapp(i,1)+nrapp(i-1,1)
          nrapp(i,2)=nrapp(i,2)+nrapp(i-1,2)
        enddo
        write(iw,*) ' '
        do i=1,nappri 
          riapp=rappri(1)+rappri(3)*(i-1)
          if(restri(3).ne.0) sepmax=min(sepmax,restri(3))
          if(riapp.gt.sepmax.and.riapp-rappri(3).le.sepmax) riapp=sepmax
          if(riapp.le.sepmax.and.sepmax.ne.0) then 
          if(correl) then
           write(iw,9205) riapp,nrapp(i,2),nrapp(i,1),3,
     *       erapp(i,1)+etotsav-erapp(i,2)-etotusav,erapp(i,2)+etotusav
          else
           write(iw,9206) riapp,nrapp(i,1),3,erapp(i,1)+etotusav
          endif
          endif
        enddo
        write(iw,*) ' '
      endif
c
c     write(iw,9200) ' + D2 ',3,etot1+eddim
      call daxpy(nfmoelm,one,fmoelm(1,nfg+1),1,fmoelm(1,nfg+2),1)
      if(iemom.ge.1) 
     *  write(iw,9230) 3,(fmoelm(ielm,nfg+2),ielm=nelm(0)+1,nelm(1)),
     *                  DNRM2(3,fmoelm(1,nfg+2),1)
      if(iemom.ge.2) 
     *  write(iw,9231) (fmoelm(ielm,nfg+2),ielm=nelm(1)+1,nelm(2))
      if(iemom.ge.3) 
     *  write(iw,9232) (fmoelm(ielm,nfg+2),ielm=nelm(2)+1,nelm(3))
c     Do not output the gradient for optimisation (it will be done elsewhere)
c
      if(nder.gt.0.and.modef.ne.1) then
        call fmogout(3,masout,indat,fmozan,fmode,DEFTF,TORQF)
      endif
c
 1000 continue
c
      IF(nfmopcm.ne.0) THEN
c
c     At present multilayer properties are not available (one layer only).
c
        write(iw,9350)
        ges=zero
        if(nbody.ge.1) then
          write(iw,9360)
          do ifg=1,nfg
            write(iw,9365) ifg,frgnam(ifg),esolv(ifg)*tokcal,
     *                     (emocdr(i,ifg),i=1,3)
            ges=ges+esolv(ifg)
          enddo
          write(iw,9355) 1,ges*tokcal
c         write(6,*) 'wwwpcm-esi',ges 
        endif
        if(nbody.ge.2) then
          write(iw,9370)
          ijfg=0
          rij=-one
          do ifg=1,nfg
            do jfg=1,ifg-1
              ijfg=ijfg+1
              if(needr.ne.0) then
c               separated dimers have no esolv set (since monomer 
c               values were not global summed before edimer).
                rij=fmodist(ifg,0,0,jfg)
                if(resdim.ne.0.and.rij.gt.resdim)
     *            esolv(nfg+ijfg)=esolv(ifg)+esolv(jfg)
c                 note that these values will be used in trimers below (FMO3)
              endif
              eesiij=esolv(nfg+ijfg)-esolv(ifg)-esolv(jfg)
              ges=ges+eesiij
              write(iw,9375) ifg,jfg,esolv(nfg+ijfg)*tokcal,
     *                               eesiij*tokcal
            enddo
          enddo
          write(iw,9355) 2,ges*tokcal
        endif
        if(nbody.ge.3) then
          write(iw,9380)
          ijkfg=0
          rrij=zero
          rrik=zero
          rrjk=zero
          do ifg=1,nfg
            do jfg=1,ifg-1
              do kfg=1,jfg-1
                ijkfg=ijkfg+1
                eijk=esolv(nfg+nfg2+ijkfg)
                if(eijk.ne.zero) then
                  indij=(ifg*ifg-3*ifg)/2+jfg+1
                  indik=(ifg*ifg-3*ifg)/2+kfg+1
                  indjk=(jfg*jfg-3*jfg)/2+kfg+1
                  deijk=eijk+esolv(ifg)+esolv(jfg)+esolv(kfg)
     *               -esolv(nfg+indij)-esolv(nfg+indik)-esolv(nfg+indjk)
                  ges=ges+deijk
                  write(iw,9385) ifg,jfg,kfg,eijk*tokcal,deijk*tokcal
                endif
              enddo
            enddo
          enddo
          write(iw,9355) 3,ges*tokcal
        endif
        GCAV=zero
        IF(ICAV.EQ.1) gcav=gcavp 
        IF(IDISP.EQ.0) THEN
          GDISP=ZERO
          GREP=ZERO
        ELSE
          GREP=GRP
        END IF
        gint=ges+(gcavp+gdisp+grep)/tokcal
        etot=etot+gint-ges
        if(correl) etotu=etotu+gint-ges
        write(iw,9390) GES*tokcal,GCAV,GDISP,GREP,GINT*tokcal
        corstr2=' un'
        if(correl) corstr2='   '
        write(iw,9395) corstr2,etot,corstr2,etot-gint
        if(correl) write(iw,9395) ' un',etotu,' un',etotu-gint
      END IF
c
      if(ifgfmo0.eq.0) then
        escf=etot
        if(nend.ge.7) escf=etotu
      endif 
c
      if(ndftfg.eq.1) then
        write(iw,9450) (etotdft(i),i=1,nedft)
      endif
      if(primul.and.nbody.gt.0.and.(.not.dotd.or.iexcit(3).gt.0)) then
c       charges are not summed properly in PL0
        nbodyc=nbody
c       if(dopleda) nbodyc=1 
c       Choose the charge meister.
        charger='Mulliken'
        if(NPTSTN.ne.0) charger='Stone' 
        write(iw,9400) charger
        do i=1,natfmo
          atmulq(i,2)=atmulq(i,2)+atmulq(i,1)
          atmulq(i,3)=atmulq(i,3)+atmulq(i,2)
          write(iw,9410) i,indat(i),fmozan(i),
     *                   (-atmulq(i,j)+fmozan(i),j=1,nbodyc)
        enddo
      endif
c     The end of if(masout).
      endif
c
c     Print FMO-TDDFT RESULTS
c
      if(dotd) call fmotdprnt(frgnam,eexcit,texcit,osmd,iexcit,nstmono,
     *                        isumd,eexfg,needr.ne.0)
c
c     Zero out the energy for diverged runs
c
      if(ndiv2+ndiv3.ne.0) then
        etot=zero 
        escf=zero
      endif
c
c     Now we broadcast the results to slaves: ESCF, ETOT, ENUCR.
c
      IF(GOPARR) then
        dumbuf(1)=etot
        dumbuf(2)=escf
        dumbuf(3)=enucr
        CALL DDI_BCAST(2421,'F',dumbuf,3,MASTER)
        etot=dumbuf(1)
        escf=dumbuf(2)
        enucr=dumbuf(3)
      endif
      if(maswrk) write(iw,9900)
      call timit(1)
      RETURN
 9100 format(//9x,'FMO properties',/9x,15(1H#),//)
 9106 format(1x,'No monomer SCF done in layer',I2)
 9108 format(1x,'Layer',I2,' DIVERGED (max iterations',I3,', conv ',
     *          E7.2,').')
 9109 format(1x,'Layer',I2,' converged in',I3,' iterations (max',I3,
     *          ', conv ',E7.2,').',/)
 9110 format(/1x,'Energy values (E) are given with ',
     *           'external ESP subtracted.')
 9111 format(/1x,'D(dipole), Q(quadrupole) and O(octopole) moments are',
     *           ' printed in the order:',
     *       /1x,'X,Y,Z; XX,YY,ZZ,XY,XZ,YZ; XXX,XXY,XXZ,XYY,YYY,YYZ,',
     *           'XZZ,YZZ,ZZZ,XYZ,',
     *       /1x,'computed at the point: ',3F13.8)
 9112 format(/1x,
     *         'Both uncorrelated and correlated energies are printed.',
     *      /1x,'Gradient is printed only as correlated values. ',
     *      /1x,'Other properties correspond to',A3,'correlated ',
     *          'wavefunction.')
 9113 format(1x,'For FMO-CC all values labelled "correlated" are for ',
     *          'the highest CC level.')
 9114 format(/1x,'One-body FMO properties.',/1x,23(1H=)/)
 9115 format(26x,'Ecorr',11x,'Euncorr',8x,'DX',8x,'DY',8x,'DZ')
 9116 format(28x,'E',11x,'DX',7x,'DY',7x,'DZ')
 8115 format(22x,'EPL0d',4x,'EPL0s',4x,'EPL0',4x,'EPL0DI',2x,'EPLd*PLd',
     *       2x,'EPLd*DI',4x,'DX',8x,'DY',8x,'DZ')
 8116 format(22x,'EPL0d',5x,'EPL0s',5x,'EPL0',5x,'EPLd*PLd',
     *        8x,'DX',7x,'DY',7x,'DZ')
 8117 format(22x,'EPL0d',4x,'EPL0s',4x,'EPL0',4x,'EPL0DI',5x,
     *           'DX',8x,'DY',8x,'DZ')
 8118 format(1x,I4,'(',A8,',L',I1,')',4F9.3,3F10.5)
 9117 format(1x,I4,'(',A8,',L',I1,')',6F9.3,3F10.5)
 9118 format(1x,I4,'(',A8,',L',I1,')',2F17.9,3F10.5)
 9119 format(1x,I4,'(',A8,',L',I1,')',4F9.3,3F10.5)
 9120 format(1x,I4,'(',A8,',L',I1,')', F17.9,3F10.5)
 8122 format(8x,'Q=',6F10.5) 
 8123 format(8x,'O=',5F10.5,/8x,2x,5F10.5)
 9121 format(/1x,'Two-body FMO properties.',/1x,24(1H=)//,
     *        1x,'DL: D=C dynamically correlated (MP2,CI), D=N not ',
     *          'dynamically correlated',/1x,'(RHF,DFT). D=S separated',
     *          ' dimer: semiclassical ','interaction (ES), D=M MCSCF.',
     *       /1x,'L stands for layer, Z is the monomer charge product,',
     *           ' R is the interfragment',/1x,'distance relative to ',
     *         'van-der-Waals radii (-1.00 is printed if distances are',
     *       /1x,'not computed). dDIJ*VIJ is the density polarisation ',
     *           'contribution.',/1x,'Q(I->J) is the charge transfer ',
     *           'amount, printed as zero if not available.',
     *       /1x,'Positive values correspond to I in IJ having extra ',
     *           'negative charge.')
 9122 format(/4x,'I    J DL  Z    R   Q(I->J)',8x,'E',11x,
     *       'EIJ-EI-EJ   dDIJ*VIJ',6x,'tot',/1x,81(1H-))
c9122 format(/4x,'I    J DL  Z    R   Q(I->J)',8x,'E',11x,
c    *       'EIJ-EI-EJ   dDIJ*VIJ    dDIJ*dVIJ',6x,'tot',/1x,90(1H-))
c9123 format(/4x,'I    J DL  Z    R   Q(I->J)',8x,'Ecorr',12x,'Euncorr',
c    *       6x,'EIJ-EI-EJ,corr/uncorr  dDIJ*VIJ    dDIJ*dVIJ tot,corr',
c    *       /1x,118(1H-))
 9123 format(/4x,'I    J DL  Z    R   Q(I->J)',8x,'Ecorr',12x,'Euncorr',
     *        6x,'EIJ-EI-EJ,corr/uncorr  dDIJ*VIJ,unc tot,corr',
     *       /1x,112(1H-))
 8120 format(/3x,'I   J  IU  SF DL  Z    R   Q(I->J)',8x,'E',11x,
     *       'EIJ-EI-EJ   dDIJ*VIJ',6x,'tot',/1x,87(1H-))
 8121 format(/3x,'I   J  IU  SF DL  Z    R   Q(I->J)',8x,'Ecorr',12x,
     *      'Euncorr',6x,'EIJ-EI-EJ,corr/uncorr  dDIJ*VIJ,unc tot,corr',
     *       /1x,118(1H-))
 9124 format(/4x,'I    J DL  Z    R   Q(I->J)  EIJ-EI-EJ dDIJ*VIJ',4x,
     *           'total',5x,'Ees ',5x,'Eex ',3x,'Ect+mix ',/1x,
     *           84(1H-))
 8124 format(/4x,'I    J DL  Z    R   Q(I->J)  EIJ-EI-EJ dDIJ*VIJ',4x,
     *           'total',3x,'EPL0I+J',4x,'Ees0',5x,'Eex0',4x,'Ectmix0',
     *       /1x,94(1H-))
 9125 format(/4x,'I    J DL  Z    R   Q(I->J)  EIJ-EI-EJ dDIJ*VIJ',4x,
     *        'total',5x,'Ees ',5x,'Eex ',3x,'Ect+mix ',3x,'Edisp ',
     *       /1x,93(1H-))
 8125 format(/4x,'I    J DL  Z    R   Q(I->J)  EIJ-EI-EJ dDIJ*VIJ',4x,
     *        'total',3x,'EPL0I+J',4x,'Ees0',5x,'Eex0',4x,'Ectmix0',3x,
     *        'Edisp0',/1x,103(1H-))
 9126 format(/1x,'Notation: BSSE corr. is a sum of int columns 3+4-1-2',
     *        8x,'/--- mon basis ---\\ /+++ dim basis +++\\ ',4x,
     *        'I   J DL',10x,'E',9x,'EIJ-EI-EJ   dDIJ*VIJ    EIJ-EI-EJ',
     *           ' dDIJ*VIJ BSSE corr.',4x,'tot',/1x,93(1H-))
 9128 format(/4x,'I   J DL    R',11x,'E',9x,'EIJ-EI-EJ corr/uncorr ',
     *       'dDIJ*VIJ',4x,'BSSE corr/uncorr',4x,'tot',/1x,88(1H-))
 9129 format(/4x,'I   J DL    R',11x,'E',9x,'EIJ-EI-EJ   dDIJ*VIJ',4x,
     *       'BSSE corr.',4x,'tot',/1x,79(1H-))
 9131 format(2I5,1x,A1,I1,I3,F7.2,F8.4,F17.9,F13.8,F12.8,F9.3)
 9132 format(2I5,1x,A1,I1,I3,F7.2,F8.4,F10.3,F9.3,F10.3,F9.3,F10.3,
     *       4F9.3)
 9133 format(2I5,1x,A1,I1,I3,F7.2,F8.4,2F17.9,2F13.8,F12.8,F9.3)
 9134 format(2I5,1x,A1,I1,I3,F7.2,F8.4,F10.3,F9.3,2F10.3,5F9.3)
 9140 format(2I5,1x,A1,I1,1x,F17.9,5F12.8,F9.3)
 9150 format(2I5,1x,A1,I1,I7,F17.9,5F12.8,F9.3)
 9152 format(2I5,1x,A1,I1,I7,F17.9,3F12.8,F9.3)
 9160 format(/1x,'Binding (BE) to the selected fragments in kcal/mol',
     *           ' (final is correlated or SCF).',/)
 9161 format(4I4,1x,A1,I1,I3,F7.2,F8.4,F17.9,F13.8,F12.8,F9.3)
 9163 format(4I4,1x,A1,I1,I3,F7.2,F8.4,2F17.9,2F13.8,F12.8,F9.3)
 9170 format(1x,I5,'(',A8,'): SCF',F11.4,', final BE',F11.4)
 9175 format(1x,I5,'(',A8,'): ES',F11.4,', SCF',F11.4,', final BE',
     *          F11.4)
 9200 format(/1x,'Total energy of the molecule: E',A6,'(',I1,')=',F19.9)
 9201 format(/1x,'Total energy, BDA corrected E',A6,'(',I1,')=',F19.9)
 9202 format(/1x,'The backbone energy EBB',A6,'(',I1,')=',F19.9)
 9203 format(/1x,'The total PIEDA energy E',A6,'(',I1,')=',F19.9)
 9205 format(1x,'R=',F6.2,' N=',2I6,' (C,N) ',I1,'-body E=',F14.9,F19.9)
 9206 format(1x,'R=',F6.2,' N=',I6,' (N) ',I1,'-body E=',F19.9)
 9210 format(//1x,'Total energy E',A6,'(2) without BSSE correction: ',
     *        F19.9/1x,24x,'with    BSSE correction: ',F19.9)
 9220 format(/1x,'Nuclear repulsion energy:   ',F22.9/,
     *        1x,'Total number of atoms:             ',I9/,
     *        1x,'Total number of electrons:         ',I9/,
     *        1x,'Total charge:                      ',I9/,
     *        1x,'Total number of basis functions:   ',I9/,
     *        1x,'Total number of molecular orbitals:',I9)
 9230 format(1x,'Dipole moment D(xyz),DA(',I1,')=',4F13.7)
 9231 format(1x,'Quadrupole moment Q(XX,YY,ZZ)=',3F12.7/,
     *       1x,'                   (XY,XZ,YZ)=',3F12.7)
 9232 format(1x,'Octopole moment O(XXX,XXY,XXZ,XYY)=',4F11.6/,
     *       1x,'                 (YYY,YYZ,XZZ,YZZ)=',4F11.6/,
     *       1x,'                         (ZZZ,XYZ)=',2F11.6)
 9250 format(/1x,'Three-body FMO properties.',/1x,24(1H=)/)
c9260 format(/4x,'I   J   K DL    RMIN   RMAX   Ecorr',10x,'deltaEIJK',
c    *        3x,'dDIJK*VIJK',2x,'ddDIJK*VIJK',3x,'tot',/1x,91(1H-))
 9260 format(/4x,'I   J   K DL   RMIN   RMAX       Ecorr',7x,'deltaEIJK'
     *        ,',corr/uncorr',3x,'dDIJK*VIJK',2x,'ddDIJK*VIJK',3x,'tot',
     *       /1x,102(1H-))
 9265 format(/4x,'I   J   K DL   RMIN   RMAX       Ecorr',7x,'deltaEIJK'
     *          ,',corr/uncorr',3x,'dDIJK*VIJK',4x,'tot',/1x,91(1H-))
 9275 format(1x,3I4,1x,A1,I1,2F7.2,F17.9,4F12.8,F9.3)
 9280 format(1x,3I4,1x,A1,I1,2F7.2,F17.9,3F12.8,F9.3)
 9350 format(/8x,'FMO-PCM properties.',/1x,38(1H=),
     *       /1x,'Ges is the energy of electrostatic interaction,',/1x,
     *       'included in the individual n-mer and total FMO energies.',
     *       /1x,'Internal energies are for solute only excluding all ',
     *           'interactions with solvent.')  
 9355 format(/1x,'Total Ges(',I1,')=',F10.3,' kcal/mol.')
 9360 format(/4x,'I',17x,'Ges',5x,'Gcav',5x,'Gdisp',4x,'Grep',
     *           ' (all in kcal/mol)')
 9365 format(I5,' (',A8,') ',4F9.3)
 9370 format(/4x,'I    J',5x,'Ges,    dGes (kcal/mol)',/1x,38(1H-))
 9375 format(2I5,2F9.3)
 9380 format(/4x,'I    J    K',5x,'Ges,    dGes (kcal/mol)',/1x,38(1H-))
 9385 format(3I5,2F9.3)
 9390 FORMAT(/1X,'ELECTROSTATIC INTERACTION  =',F12.3,' KCAL/MOL'/
     *        1X,'PIEROTTI CAVITATION ENERGY =',F12.3,' KCAL/MOL'/
     *        1X,'DISPERSION FREE ENERGY     =',F12.3,' KCAL/MOL'/
     *        1X,'REPULSION FREE ENERGY      =',F12.3,' KCAL/MOL'/
     *        1X,'TOTAL INTERACTION          =',F12.3,' KCAL/MOL')
 9395 format(/1x,'    Free ',A3,'correlated energy in solvent=',F19.9,
     *       /1x,'Internal ',A3,'correlated energy in solvent=',F19.9)
 9400 format(/1x,'n-body ',A8,' atomic charges Q(n)',
     *       //4x,'IAT  IFG   Z',7x,'Q(1)',8x,'Q(2)',8x,'Q(3)')
 9410 format(1x,2I5,F6.1,3F12.6)
 9420 format(/1x,'Charge transfer for each fragment:',
     *      //2x,'IFG QFG  DeltaQ     and its contributions from JFG, ',
     *           'Q(JFG->IFG).',/1x,78(1H-))
 9425 format(I5,I3,F9.4,' =',4(I5,'->',F8.4))
 9430 format(/1x,'Total absolute monomer transf. charge   =',F12.6,
     *       /1x,'Total amount of absolute transf. charge =',F12.6)
c9420 format(1x,'Donor+acceptor ',A6,' energy sum is ',F19.9)
 9450 format(/1x,'DFT exchange+correlation energy=',F19.9,/,
     *       1x,'Total electron number          =',F19.9)
 9500 format(/1x,I6,' of ',I1,'-body energies diverged!!!',
     *       ' The total properties are not meaningful.')
 9600 format(/1x,'PIEDA summary',
     *       /5x,'Main relations:'
     *       /5x,'EIB0 = EPL0 + EPL0DI + EES0 + EEX0 + ECTmix0 + EDI0',
     *       /5x,'EIB  = EIB0 + EPL*EX + EPL*CTmix + EPL*DI',
     *       /5x,'EIB  = Eint + EPL0d + EPL0DI + EPLd*PLd + EPLd*DI',
     *       /5x,'EIB  = E(FMO2) - E(FMO0) - E(BDA)', 
     *       /5x,'EIB0 - interaction relative to free state using free',
     *           ' state densities,',
     *       /5x,'EIB  - interaction relative to free state using PL ',
     *           'state densities,',/,
     *       /1x,'Polarisation (destabilisation)  EPL0d',F12.3,
     *       /1x,'Polarisation (stabilisation)    EPL0s',F12.3,
     *       /5x,'Polarisation (total, free state)    EPL0   ',F12.3,
     *       /5x,'Pol.-induced dispersion(free state) EPL0DI ',F12.3,
     *       /5x,'Electrostatic (free state)          EES0   ',F12.3,
     *       /5x,'Exchange (free state)               EEX0   ',F12.3,
     *       /5x,'Charge transfer+all mixing(free s.) ECTmix0',F12.3,
     *       /5x,'Dispersion (free state)             EDI0   ',F12.3,
     *       /9x,'Internal binding (free state)        EIB0     ',F12.3,
     *       /9x,'(Uncorr. int. bind.,free s.) EIB0-EPL0DI-EDI0 ',F12.3,
     *           ')',/,
     *       /9x,'Polarisation-exchange                EPL*EX   ',F12.3,
     *       /9x,'Polarisation-charge transfer         EPL*CTmix',F12.3,
     *       /5x,'Polarisation-dispersion (destab)    EPLd*DI',F12.3,
     *       /5x,'Polarisation-dispersion (stab)      EPLs*DI',F12.3,
     *       /9x,'Polarisation-dispersion (total)      EPL*DI   ',F12.3,
     *      /13x,'Internal binding                     EIB      ',F12.3,
     *      /13x,'(Uncorr. int. bind.) EIB-EPL0DI-EDI0-EPL*DI   ',F12.3,
     *           ')',/)
 9605 format(/1x,'Polarisation-polarisation (destab)  EPLd*PLd',F12.3,
     *       /1x,'Polarisation-polarisation (stab)    EPLs*PLs',F12.3,
     *       /1x,'Polarisation-polarisation (total)    EPL*PL ',F12.3,
     *       /5x,'Polarisation (total)        EPL=EPL0+EPL*PL ',F12.3,
     *       /1x,'(EPL*PL up to two-body terms is included in ECTmix0,',
     *       /1x,'the rest of EPL*PL is in EPL*CTmix. E(CT+mix) does',
     *           ' not include any of EPL*PL.)',
     *       /1x,'Estimate of free state CT without PL*PL, E(CT0+mix0)',
     *           F12.3,
     *       /1x,'Total monomer destabilisation,   E(PLd+corr)',F12.3)
 9610 format(/1x,'Interaction energy relative to PL state:'
     *       /1x,'Eint = EES + EEX + E(CT+mix) + EDI',
     *       /1x,'Eint = E(FMO2) - E(FMO1) - E(BDA)',
     *       /1x,'Eint - interaction relative to PL state using PL ',
     *           'state densities.',/,
     *       /5x,'Electrostatic (PL state, incl. EPLs)  EES  ',F12.3,
     *       /5x,'Exchange (PL state)                   EEX  ',F12.3,
     *       /5x,'Charge transfer (PL state)        E(CT+mix)',F12.3,
     *       /5x,'Dispersion (PL state)                 EDI  ',F12.3,
     *       /9x,'Total interaction (PL state)         Eint  ',F12.3,/)
 9620 format(/1x,'Unconnected interaction energy relative to PL state:',
     *       /1x,'E''int = E''ES + E''EX + E''(CT+mix) + E''DI',
     *       /1x,'E''int = E''(FMO2) - E''(FMO1)',
     *       /1x,'E''int - interaction relative to PL state using PL ',
     *           'state densities, ',
     *       /1x,'excluding contributions from dimers between which ',
     *       /1x,'a covalent bond is fractioned.',/,
     *       /5x,'Electrostatic (PL state, incl. EPLs)  E''ES ',F12.3,
     *       /5x,'Exchange (PL state)                   E''EX ',F12.3,
     *       /5x,'Charge transfer (PL state)      E''(CT+mix) ',F12.3,
     *       /5x,'Dispersion (PL state)                 E''DI ',F12.3,
     *       /9x,'Total interaction (PL state)         E''int ',F12.3,/,
     *       /1x,'BDA energies',
     *       /5x,'Electrostatic (BDA)              EESBDA',F16.3,
     *       /5x,'Exchange (BDA)                   EEXBDA',F16.3,
     *       /5x,'Charge transfer (BDA)      E(CT+mix)BDA',F16.3,
     *       /5x,'Dispersion (BDA)                 EDIBDA',F16.3,
     *       /9x,'Total interaction (BDA)         EintBDA',F16.3,/)
 9900 format(/1x,'Done with FMO properties.')
      END
c
C*MODULE fmoio   *DECK vdwrout
      SUBROUTINE vdwrout(fmozan,vdwrad)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 c8dum
      character*4 ATOMNM
      PARAMETER (MAXNZ=137,units=0.52917724924D+00)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension fmozan(*),vdwrad(MAXNZ),nout(MAXNZ)
      data ATOMNM/'    '/,c8dum/' '/
c
c     Output van der Waals radii.
c
      call viclr(nout,1,MAXNZ)
c
      write(iw,9000)
      do iat=1,natfmo
        iz=int(fmozan(iat)+0.1D+00)
        if(nout(iz).eq.0) then
          call zsymnum(c8dum,ATOMNM,iz)
          write(iw,9100) ATOMNM,vdwrad(iz)*units
          nout(iz)=1
        endif
      enddo
      return
 9000 format(4x,'Van der Waals radii in Angstrom are:')
 9100 format(7x,'R(',A3,')=',F8.4) 
      END
c
C*MODULE fmoio   *DECK eminout
      subroutine eminout(imode,ilay,ifgfmo0,layfrg,emon,some)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 gprnam,end
      LOGICAL GOPARR,DSKWRK,MASWRK,some
c     parameter (sawdust=1.0D-25)
      dimension layfrg(*),emon(nfg,4,nlayer)
      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
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
c
      dimension ewrk(4)
c
      data gprnam/' $FMOENM'/,end/' $END'/
c
c     read (imode=0) or write (imode=1) monomer energies
c     ilay is ignored for imode=0.
c
c     nen=1
c     if(nbsse.eq.3) nen=nen+1
c     if(rcorsd.ne.0) nen=nen+nen
      nen=4
      if(imode.eq.0) then
        CALL SEQREW(IR)
        CALL FNDGRP(IR,gprnam,IEOF)
        if(ieof.ne.0) then
          if(maswrk) write(iw,*) gprnam,' not found'
c         abort if not doing monomer SCF
          if(irststp.gt.2.or.irststp.eq.2.and.nbsse.eq.3) call abrt
c         for BSSE one needs free monomer energies
          return
        endif
        if(maswrk) then
  100   continue
          READ(IR,9010,END=120,ERR=120) ifg,ilay,(ewrk(ien),ien=1,nen)
          if(ifg.gt.nfg.or.ilay.gt.nlayer) then
            if(maswrk) write(6,*) 'bad indices read in',gprnam,ifg,ilay
            call abrt
          endif
c         this is to avoid double counting; assign all energies to the
c         grand master. Avoid parallel confusion by assigning some sawdust
c         to other masters.
c         At present no need for sawdusting monomers.
          if(some) write(iw,9000) ifg,ilay,(ewrk(ien),ien=1,nen)
c         if(mygroup.ne.0) call dacopy(nen,sawdust,ewrk,1)
          call dcopy(nen,ewrk,1,emon(ifg,1,ilay),nfg)
c         write(iw,9000) ifg,ilay,(emon(ifg,ien,ilay),ien=1,nen)
        goto 100
  120   continue 
        endif
        IF (GOPARR) CALL DDI_BCAST(2421,'F',emon,nfg*4*nlayer,MASTER)
      else
        if(maswrk) then
          write(ip,9100) gprnam 
          if(ifgfmo0.ne.0) then
            if(layfrg(1).ge.ilay) write(ip,9010) ifgfmo0,ilay,
     *                                   (emon(1,ien,ilay),ien=1,nen)
          else
            do ifg=1,nfg
              if(layfrg(ifg).ge.ilay) write(ip,9010) ifg,ilay,
     *                                   (emon(ifg,ien,ilay),ien=1,nen)
            enddo
          endif
          write(ip,9100) END
        endif
      endif
      return
 9000 format(1x,'E ifg=',I5,'ilay=',I2,' E=',4F14.10)
 9010 format(1x,I5,I2,4F18.10)
 9100 format(a8)
      end
c
C*MODULE fmoio   *DECK edin
      subroutine edin(nfg2,nen,layfrg,edim,some)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 gprnam
      parameter (zero=0.0D+00,sawdust=1.0D-25)
      dimension layfrg(*),edim(nfg2,*)
      LOGICAL GOPARR,DSKWRK,MASWRK,isgddi,parout,INITGDDI,some
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT,INITGDDI
      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
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      dimension ewrk(9)
      data gprnam/' $FMOEND'/
c
c     read dimer energies
c     this does not work with dolat (change the leading dim in edim).
c
c     nen=nen0
c     if(nbsse.ne.0.and.IRSTSTP.eq.5) nen=nen-2
c     if(imode.eq.0) then
      CALL SEQREW(IR)
      CALL FNDGRP(IR,gprnam,IEOF)
      if(ieof.ne.0) then
        write(iw,*) gprnam,' not found: all dimers will be run.'
c       abort if not doing dimer SCF
        if(irststp.gt.4) then
          write(6,*) 'bad irststp',irststp
          call abrt
        endif
        return
      endif
      if(maswrk) then
      nread=0
  100 continue
c       if(nen.le.4) then
        READ(IR,9210,END=120,ERR=120) ifg,jfg,ilay,(ewrk(ien),ien=1,nen)
c       else
c       READ(IR,9215,END=120,ERR=120) ifg,jfg,ilay,(ewrk(ien),ien=1,nen)
c       endif
        if(ifg.le.1.or.ifg.gt.nfg.or.jfg.gt.nfg.or.jfg.gt.ifg.or.
     *     ilay.gt.nlayer) then
          write(6,*) 'bad indices read in',gprnam,ifg,jfg,ilay
          call abrt
        endif
        nread=nread+1
c       ignore energies for lower levels
        if(ilay.lt.layfrg(ifg).and.ilay.lt.layfrg(jfg)) goto 100
        if(some) then
          if(nen.le.3) then 
            write(iw,9000) ifg,jfg,ilay,ewrk(1)
          else
            write(iw,9000) ifg,jfg,ilay,ewrk(1),(ewrk(ii),ii=3,4)
          endif
        endif
c       this is to avoid double counting and prevent recompututing
c       if(meglob.ne.0) call dacopy(nen,sawdust,ewrk,1)
        if(mygroup.ne.0) call dacopy(nen,sawdust,ewrk,1)
c       ilay is ignored in saving data. This is because we expect the 
c       layers be given in ascending order, so that higher layers if exist
c       should replace the data for the lower ones. No check is done.
c       In fact ilay is punched for the user to check the data.
        if(ifg.lt.jfg) then
          write(6,*) 'bad index pair',ifg,jfg
          call abrt
        endif 
        ijfg=((ifg-1)*(ifg-2))/2+jfg
        call dcopy(nen,ewrk,1,edim(ijfg,1),nfg2)
        if(nbsse.ne.0.and.IRSTSTP.eq.5) then
          edim(ijfg,4)=zero
          edim(ijfg,5)=zero
cnb       correlation (6)? 
        endif
      goto 100
  120 continue 
      write(iw,9100) nread,nfg2
      endif
      IF (GOPARR) then
        CALL DDI_BCAST(2422,'F',edim,nfg2*nen,MASTER)
c       Slaves in group zero should have energies replaced by sawdust to avoid
c       double counting.
c       This code is not correct for BSSE restarts.
        if(nbsse.ne.0.and.IRSTSTP.eq.5) call abrt 
c
        if(mygroup.eq.0.and..not.maswrk) then
          do ijfg=1,nfg2
            if(edim(ijfg,1).ne.zero) 
     *        call dcopy(nen,sawdust,1,edim(ijfg,1),nfg2)
          enddo
        endif
      endif
      return
 9000 format(1x,'E ifg=',I5,' jfg=',I5,' ilay=',I2,' E=',3F14.7)
 9100 format(1x,i10,' out of',I10,' dimer energies read in.')
 9210 format(1x,2I5,I3,3F22.10,2(/14x,3F22.10))
c9215 format(1x,2I5,I3,2F22.10,/14x,2F22.10)
      end
c
C*MODULE fmoio   *DECK etin
      subroutine etin(nfg3,nent,layfrg,etrim,some)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*8 gprnam
      parameter (zero=0.0D+00,sawdust=1.0D-25)
      dimension layfrg(*),etrim(nfg3,*)
      LOGICAL GOPARR,DSKWRK,MASWRK,isgddi,parout,INITGDDI,some
      COMMON /GDDI/   ISCOPE,NGROUPS,MYGROUP,MEGLOB,npglob,nnglob,JBTYP,
     *                ISGDDI,PAROUT,INITGDDI
      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
      common /fmoopt/ espsca(9),RESPAP(2),rESPPC(2),rESDIM,restri(4),
     *                rcorsd,respct,convfg,cnvdmp,coroff,IXESP,mxitfg,
     *                nguess,NBSSE,modorb,modpar,irststp,irstlay,nprfmo,
     *                nfmopal,modprp,maxl1c,ipieda,modgrd,modesp,ivmul
      dimension ewrk(9)
      data gprnam/' $FMOENT'/
c
c     read trimer energies
c
      CALL SEQREW(IR)
      CALL FNDGRP(IR,gprnam,IEOF)
      if(ieof.ne.0) then
        write(iw,*) gprnam,' not found: all trimers will be run.'
c       abort if not doing dimer SCF
        if(irststp.gt.9) then
          write(6,*) 'bad irststp',irststp
          call abrt
        endif
        return
      endif
      if(maswrk) then
      nread=0
  100 continue
        READ(IR,9210,END=120,ERR=120) ifg,jfg,kfg,ilay,
     *                                (ewrk(ien),ien=1,nent)
        if(ifg.le.1.or.ifg.gt.nfg.or.jfg.gt.ifg.or.kfg.gt.jfg.or.
     *     ilay.gt.nlayer) then
          write(6,*) 'bad indices read in',gprnam,ifg,jfg,kfg,ilay
          call abrt
        endif
        nread=nread+1
c       ignore energies for lower levels
        if(ilay.lt.layfrg(ifg).and.ilay.lt.layfrg(jfg)) goto 100
        if(some) then
          if(nent.le.3) then 
            write(iw,9000) ifg,jfg,kfg,ilay,ewrk(1)
          else
            write(iw,9000) ifg,jfg,kfg,ilay,ewrk(1),(ewrk(ii),ii=3,4)
          endif
        endif
c       this is to avoid double counting and prevent recompututing
        if(mygroup.ne.0) call dacopy(nent,sawdust,ewrk,1)
c       ilay is ignored in saving data. This is because we expect the 
c       layers be given in ascending order, so that higher layers if exist
c       should replace the data for the lower ones. No check is done.
c       In fact ilay is punched for the user to check the data.
        ijkfg=((ifg-1)*(ifg-2)*(ifg-3))/6+((jfg-1)*(jfg-2))/2+kfg
        call dcopy(nent,ewrk,1,etrim(ijkfg,1),nfg3)
      goto 100
  120 continue 
      write(iw,9100) nread,nfg3
      endif
      IF (GOPARR) then
        CALL DDI_BCAST(2422,'F',etrim,nfg3*nent,MASTER)
c       Slaves in group zero should have energies replaced by sawdust to avoid
c       double counting.
c
        if(mygroup.eq.0.and..not.maswrk) then
          do ijkfg=1,nfg3
            if(etrim(ijkfg,1).ne.zero) 
     *        call dcopy(nent,sawdust,1,etrim(ijkfg,1),nfg3)
          enddo
        endif
      endif
      return
 9000 format(1x,'E ifg=',I5,' jfg=',I5,' kfg=',I5,' ilay=',I2,' E=',
     *       3F14.7)
 9100 format(1x,i10,' out of',I10,' trimer energies read in.')
 9210 format(1x,3I5,I3,2F22.10,2F16.10)
      end
c
C*MODULE fmoio   *DECK PRTRILe
      SUBROUTINE PRTRILe(D,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      DIMENSION D(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT SYMMETRIC MATRIX -D- OF DIMENSION -N- -----
C     THE ROWS WILL BE LABELED WITH BASIS FUNCTION TAGS.
C
      IF (MASWRK) THEN
      MAX = 5
      IF (NPRINT .EQ. 6) MAX = 10
      MM1 = MAX - 1
      DO 120 I0=1,N,MAX
         IL = MIN(N,I0+MM1)
         WRITE(IW,9008)
         WRITE(IW,9028) (I,I=I0,IL)
         WRITE(IW,9008)
         IL = -1
         DO 100 I=I0,N
            IL=IL+1
            J0=I0+(I*I-I)/2
            JL=J0+MIN(IL,MM1)
            WRITE(IW,9048) I,BFLAB(I),(D(J),J=J0,JL)
  100    CONTINUE
  120 CONTINUE
      END IF
      RETURN
 9008 FORMAT(1X)
 9028 FORMAT(15X,10(4X,I4,3X))
 9048 FORMAT(I5,2X,A8,10E15.6)
      END
C*MODULE fmoio   *DECK fmoccp
      SUBROUTINE fmoccp(ibody,escf,extracc)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*13 CMET(9)
      LOGICAL CNVR,QDR,QDCR
      COMMON /CCENGY/ ENRG,EREF,EMP2,eccn,ETOT(6),ECORR(6),
     *                DIAGS(3),AMPMX(5,2),IAMPMX(5,4,2),XO1,XO2,
     *                DIFMAX,DIFFENG,ITER,CNVR
      COMMON /CCINFO/ TSH,NH,NP,MET,MEM,ICONV,MAXIT,IREST,IDISC
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON/T4TOT/ETT2,ETS2,EQ1A,EQ1B,EQ2A,EQ2B,EQ3A,EQ3B,EQ4A,
     *EQ4B,O4A,O4B,ECORQ(8),ETOTQ(8)
      COMMON/QUADRUPLE/QDR,QDCR
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      DATA CMET/'LCCD,','CCD,','CCSD,','CCSD[T],','CCSD(T),',
     *          'R-CCSD[T],','R-CCSD(T),','CR-CCSD[T],','CR-CCSD(T),'/
      data CRCCL/8HCR-CCL  /
      dimension extracc(*)
C
c     This is a mutilated clone of DPRINT.
c     Note that it is assumed that there is no 2nd order correlated density
c     contribution of the form Tr(deltaDcorr*V) to n-mers, which is normally
c     so unless some new method such as FMO1-RCC (relaxed CC) is developed.
c
      IMET=MET+1
      IF(MET.GT.2) IMET=3
c
      EREF=escf
      EMP2=extracc(1)
      eccn=extracc(2)
      call dcopy(6,extracc(3),1,ECORR,1)
      if(qdr) call dcopy(7,extracc(9),1,ECORQ(2),1)
c
      write(iw,9000)
      WRITE(iw,1000) ibody,'MP2,         ',EMP2,eref+EMP2
      WRITE(iw,1000) ibody,CMET(IMET),eccn,eref+ECCN
      IF(MET.GT.2) THEN
        WRITE(iw,1000) ibody,'CCSD[T],     ',ECORR(1),eref+ECORR(1)
        WRITE(iw,1000) ibody,'CCSD(T),     ',ECORR(2),eref+ECORR(2)
        IF(MET.GT.3) THEN
          WRITE(iw,1000) ibody,'R-CCSD[T],   ',ECORR(3),eref+ECORR(3)
          WRITE(iw,1000) ibody,'R-CCSD(T),   ',ECORR(4),eref+ECORR(4)
          IF(MET.GT.4.AND..NOT.(QDR.AND..NOT.QDCR)) THEN
            WRITE(iw,1000) ibody,'CR-CCSD[T],  ',ECORR(5),eref+ECORR(5)
            WRITE(iw,1000) ibody,'CR-CCSD(T),  ',ECORR(6),eref+ECORR(6)
          END IF
        END IF
      END IF
c
      IMET=MET+1
      IF(MET.GT.2) IMET=MET+2
      IF(QDR) THEN
        WRITE(iw,1000) ibody,'CCSD(TQ),B   ',ECORQ(2),eref+ECORQ(2)
        WRITE(iw,1000) ibody,'R1-CCSD(TQ),A',ECORQ(3),eref+ECORQ(3)
        WRITE(iw,1000) ibody,'R1-CCSD(TQ),B',ECORQ(4),eref+ECORQ(4)
        WRITE(iw,1000) ibody,'R2-CCSD(TQ),A',ECORQ(5),eref+ECORQ(5)
        WRITE(iw,1000) ibody,'R2-CCSD(TQ),B',ECORQ(6),eref+ECORQ(6)
        IF(QDCR) THEN
          WRITE(iw,1000) ibody,'CR-CCSD(TQ),A',ECORQ(7),eref+ECORQ(7)
          WRITE(iw,1000) ibody,'CR-CCSD(TQ),B',ECORQ(8),eref+ECORQ(8)
        END IF
      END IF
      if(CCTYP.EQ.CRCCL) then
        e1=eccn+extracc(16)
        e2=eccn+extracc(17)
        WRITE(iw,1000) ibody,'CCSD(2)_T,   ',e1,e1+eref
        WRITE(iw,1000) ibody,'CR-CCSD(T)_L,',e2,e2+eref
      END IF
c
      write(iw,9000)
      return
C
 1000 FORMAT(1X,'FMO',I1,'-',A13,' corr. energy=',F15.9,
     *       ' total E=',F19.9)
 9000 FORMAT(1X)
      END
c
C*MODULE fmoio   *DECK fmogout
      SUBROUTINE fmogout(ibody,masout,indat,fmozan,fmode,DEFTF,TORQF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      logical masout
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      parameter (one=1.0D+00,MXFRG=1050,MXFGPT=12000,MXDFG=5,
     *           MXDPPT=MXFRG*MXDFG*12)
      COMMON /EFPFMO/ IEFPFMO,ISWENUC,REPNUCEFP,EFPEFP,RUNEFP
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      Common /fmoinf/ nfg,nlayer,natfmo,nbdfg,naotyp,nbody
      dimension indat(*),fmozan(*),fmode(3,natfmo,*),DEFTF(3,nfrg,*),
     *          TORQF(3,nfrg,*),DEFTA(3),dnam(3)
      DATA DNAM /4HE'X ,4HE'Y ,4HE'Z /
c
      write(iw,9300) ibody
c     CALL EGOUT(fmode,NATfmo)
c     EGOUT cannot be used due to ANAM and BNAM being undefined.
      if(masout) then
        WRITE(IW,9000) (DNAM(J),J=1,3)
        DO I=1,NATfmo
          WRITE(IW,9010) I,indat(i),fmozan(i),(fmode(J,I,ibody),J=1,3)
        enddo
        if(IEFPFMO.eq.1) then
          do i=1,3
            DEFTA(i)=ddot(nfrg,DEFTF(i,1,ibody),3,one,0)
          enddo
          WRITE (IW,8000) ibody
          CALL EFOUT(DEFTF(1,1,ibody),DEFTA,TORQF(1,1,ibody))
        endif
      endif
      if(iefpfmo.eq.1) then
        call EGMAX(fmode(1,1,ibody),natfmo*3,DEFTF(1,1,ibody),
     *             TORQF(1,1,ibody),NFRG,FMAX,FRMS)
      else
        call EGMAX(fmode(1,1,ibody),natfmo*3,DEFT,TORQ,NFRG,FMAX,FRMS)
      endif
      WRITE (IW,9030) ibody,FMAX,FRMS
c
      RETURN
 8000 FORMAT(/21X,29('-')/
     *        21X,'FRAGMENT GRADIENT INFORMATION, FMO',I1,/
     *        21X,29('-')/)
 9000 FORMAT(/2X,'ATOM# FRG#   Z',12X,A4,13X,A4,13X,A4)
 9010 FORMAT(1X,2I5,F6.1,3F18.9)
 9030 FORMAT(/6X,'(',I1,') ',19HMAXIMUM GRADIENT = ,F10.7,4X,
     *           'RMS GRADIENT =',F10.7)
 9300 format(/1x,'Energy gradient (hartree/bohr), no BSSE: G(',I1,')')
      END
