PROGRAM FCF
c     compile by FORTRAN 90
c     Frank-Condon factors
c     maxb ... maximum excitation
c     M3 ... number of oscillators
c     nk0 ... numer of k-points
c     dimension of the m-vector batch
      IMPLICIT none
      integer*4 maxb,nk0,M3,maxe,time,tim,NB0
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 maxn,ii,jj,j,k,kk,nk,n,iq,nv(M3),mv(M3),LEXCL,
     1N7,m1,jq,ic,i,imaxu,INDX(M3),GEXCL,EEXCL,IWR,nlim,I2e,J1e,J2e,
     1state(M3),staten(M3),IEX,it,NEXC,NNEXC,sm(maxe),sn(maxe),
     1ie0,mm,units,NP,ifc,kt,kmax,NOH,ig,NG1,NG2,NE1,NE2,I1g,I2g,I1e,
     1kkr(maxb,2*nk0+1),IK,IG1,I0,IE,IG2,NDIM,NST
      real*8 pi,p,ank,ak,akfac,tp,stp,air,aii,ai,BOXF,OH,Cwork(M3),
     1xkr(maxb,2*nk0+1),xki(maxb,2*nk0+1),q,xq,fik,box,Determ,
     1w(M3,M3),ww(M3,M3),U(M3),UU(M3),e,es(M3),fcf1,fcf2,
     1maxu,fcf3,wt(M3,M3),RSCALE,
     1TEMP,Boltzmann,ktcm,EDIF,STRENGTH,EN,
     1EE,wcm,wev,wnm,zpeg,zpee,Wi(M3,M3),fcf4,fcf5,fcf6,
     1om,XMIN,XMAX,DX,y,yt,
     1ykr(maxb,2*nk0+1),yki(maxb,2*nk0+1),ifct
      real*8,allocatable::akb(:),eeb(:),S(:),BIGC(:)
      integer*4,allocatable::mvb(:)
      common/ho/p(maxb,maxb),ank(maxb),e(M3)
      logical lex,ltest,lspec,ltab,LG,lanharm
      character*4 so,unitc(3)
      data unitc/'  eV','  nm','cm-1'/
      integer*4 ix,iy,iz
      common /randc/ ix, iy, iz
      integer*4 N7g,Ng,LEXCLg,NHBIGg,MVg,KGg,ISBg
      integer*4,allocatable::ICCg(:,:),INg(:,:),ISg(:),ISIg(:,:)
      real*8 WDUMM
      real*8,allocatable::CFIg(:),EAHg(:),Wg(:)
      integer*4 N7e,Ne,LEXCLe,NHBIGe,MVe,KGe,ISBe,IDUM
      integer*4,allocatable::ICCe(:,:),INe(:,:),ISe(:),ISIe(:,:)
      real*8,allocatable::CFIe(:),EAHe(:),We(:)
c     Qs excited state normal mode coordinates, qs dimensionless
c     Ql ground state normal mode coordinates, qs dimensionless
c     Qs,i=wi,j Ql,j + Ui
c     qs,i=sqrt(esi/elj)wi,j Ql,j +sqrt(esi) Ui =wwik qli + uui
c
      NG1=1
      NG2=1
      NE1=0
      NE2=0
      ix=1
      iy=1000
      iz=99
      N7=0
      N=0
      tim=time()
      NB0=200
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      stp=dsqrt(tp)
      ai=0
      do 108 i=1,M3
      do 108 j=1,M3
108   w(i,j)=0.0d0
      do 109 i=1,M3
      nv(i)=1
      mv(i)=1
      u(i)=0.0d0
      w(i,i)=1.0d0
      e(i)=1.0d0
109   es(i)=1.0d0
      ie0=300
      ltest=.false.
      N7=1
      ifc=1
      LEXCL=20
      nlim=6
      IWR=1
      NK=6
      GEXCL=1
      EEXCL=2
      M1=10
      BOXF=1.0d0
      TEMP=298.0d0
      NOH=3
      EDIF=7.0d0
      STRENGTH=1.0d0
      lspec=.true.
      LG=.false.
      ltab=.true.
      XMIN=0.3d0
      NP=1000
      XMAX=2.0d0
      DX=0.006d0
      kmax=4
c     UNITS: 1(eV) 2(nm) 3(cm-1):
      UNITS=1
      RSCALE=1.0d0
      maxn=0


      inquire(file='g.scr',exist=lanharm)
      if(lanharm)then
c      g.scr is EIGEN.SCR from s4 for ground state
       OPEN(4,FILE='g.scr',FORM='UNFORMATTED')
       READ(4)N7g,Ng,LEXCLg,NHBIGg,MVg,KGg,ISBg
       write(6,6401)N7g,Ng,LEXCLg,NHBIGg,MVg,KGg,ISBg
6401   format('N7:',i4,'N:',i4,'LEXCL:',i4,'NHBIG:',i8,'MV:',
     1 i4,'KG:',i4,'ISB:',i4)
       allocate(ICCg(NHBIGg,LEXCLg),INg(NHBIGg,LEXCLg),ISg(NHBIGg),
     1 CFIg(NHBIGg),ISIg(KGg,4),EAHg(NHBIGg))
       READ(4)(IDUM,I=1,NHBIGg) ,(ISg(I),I=1,NHBIGg),
     1 ((ICCg(I,J),I=1,NHBIGg),J=1,LEXCLg) ,
     1 ((INg(I,J),I=1,NHBIGg),J=1,LEXCLg),
     2 (EAHg(I),I=1,MVg),(Wg(I),I=N7g,Ng),(WDUMM,I=N7g,Ng),
     3 ((IDUM,K=1,2),J=1,ISBg),((ISIg(J,K),K=1,4),J=1,KGg)
       CLOSE(4)
       write(6,*)'Anharmonic parameters read from g.scr'

       OPEN(4,FILE='e.scr',FORM='UNFORMATTED')
       READ(4)N7e,Ne,LEXCLe,NHBIGe,MVe,KGe,ISBe
       write(6,6401)N7e,Ne,LEXCLe,NHBIGe,MVe,KGe,ISBe
       allocate(ICCe(NHBIGe,LEXCLe),INe(NHBIGe,LEXCLe),ISe(NHBIGe),
     1 CFIe(NHBIGe),ISIe(KGe,4),EAHe(NHBIGe))
       READ(4)(IDUM,I=1,NHBIGe) ,(ISe(I),I=1,NHBIGe),
     1 ((ICCe(I,J),I=1,NHBIGe),J=1,LEXCLe) ,
     1 ((INe(I,J),I=1,NHBIGe),J=1,LEXCLe),
     2 (EAHe(I),I=1,MVe),(We(I),I=N7e,Ne),(WDUMM,I=N7e,Ne),
     3 ((IDUM,K=1,2),J=1,ISBe),((ISIe(J,K),K=1,4),J=1,KGe)
       CLOSE(4)
       write(6,*)'Anharmonic parameters read from e.scr'
       NE1=1
       NE2=NHBIGe
      endif

      inquire(file='FCF.OPT',exist=lex)
      if(lex)then
       write(6,*)'FCF.OPT found'
       open(9,file='FCF.OPT')
92     read(9,919,end=9191,err=9191)so
919    format(a4)
       if(so.eq.'TEMP')read(9,*)TEMP
       if(so.eq.'RSCA')read(9,*)RSCALE
       if(so.eq.'SPEC')read(9,*)lspec
       if(so.eq.'TABL')read(9,*)ltab
       if(so.eq.'UNIT')read(9,*)UNITS
       if(so.eq.'IE0 ')read(9,*)ie0
       if(so.eq.'IWR ')read(9,*)IWR
       if(so.eq.'IFC ')read(9,*)ifc
       if(so.eq.'GAUS')read(9,*)LG
       if(so.eq.'TEST')read(9,*)ltest
       if(so.eq.'EDIF')read(9,*)EDIF
       if(so.eq.'GEXC')read(9,*)GEXCL
       if(so.eq.'EEXC')read(9,*)EEXCL
       if(so.eq.'STRE')read(9,*)STRENGTH
       if(so.eq.'LEXC')read(9,*)LEXCL
       if(so.eq.'NLIM')read(9,*)nlim
       if(so.eq.'NK  ')read(9,*)NK
       if(so.eq.'NP  ')read(9,*)NP
       if(so.eq.'XMIN')read(9,*)XMIN
       if(so.eq.'NB0 ')read(9,*)NB0
       if(so.eq.'XMAX')read(9,*)XMAX
       if(so.eq.'KMAX')read(9,*)kmax
       if(so.eq.'DX  ')read(9,*)DX
       if(so.eq.'BOXF')read(9,*)BOXF
       if(so.eq.'M1  ')read(9,*)M1
       if(so.eq.'N7  ')read(9,*)N7
       if(so.eq.'N   ')read(9,*)N
       if(so.eq.'NOH ')read(9,*)NOH
       goto 92
9191   close(9)
      endif

      allocate(S(NP))

      inquire(file='fg.inp',exist=lex)
      if(lex)then
       inquire(file='fe.inp',exist=lex)
       if(lex)then
        write(6,*)'fg.inp and fe.inp found'
        call ALLIGN(N7,N,WW,UU,e,es,WI,IWR,WT)
       endif
      endif


      inquire(file='NM.TXT',exist=lex)
      if(lex)then
       write(6,*)'NM.TXT found'
       open(9,file='NM.TXT')
       read(9,*)(nv(i),i=N7,N)
       read(9,*)(mv(i),i=N7,N)
       close(9)
       maxn=1
       do 112 i=N7,N
       if(mv(i).gt.maxn)maxn=mv(i)
112    if(nv(i).gt.maxn)maxn=nv(i)
      endif


      inquire(file='E.TXT',exist=lex)
      if(lex)then
       write(6,*)'E.TXT found'
       open(9,file='E.TXT')
       do 110 iq=N7,N
110    read(9,*)e(iq),es(iq)
       close(9)
      endif

      inquire(file='U.TXT',exist=lex)
      if(lex)then
       write(6,*)'U.TXT found'
       open(9,file='U.TXT')
       read(9,*)(uu(iq),iq=N7,N)
       close(9)
      endif
      
      inquire(file='W.TXT',exist=lex)
      if(lex)then
       write(6,*)'W.TXT found'
       open(9,file='W.TXT')
       do 114 iq=N7,N
       read(9,*)(ww(iq,jq),jq=N7,N)
       do 114 jq=N7,N
       wi(iq,jq)=ww(iq,jq)
114    wt(iq,jq)=ww(iq,jq)
       close(9)
      endif

      maxu=0.0d0
      imaxu=1
      do 1131 iq=N7,N
      if(dabs(uu(iq)).gt.maxu)then
       maxu=dabs(uu(iq))
       imaxu=iq
      endif
1131  continue

      write(6,6500)TEMP,unitc(UNITS),IFC,lspec,ltab,ltest
6500  format('    TEMP: ',f8.2,' UNITS:',4x,a4,'  IFC:',i8,/,
     1       '    SPEC: ',l8  ,'  TABL:',l8   ,' TEST:',l8)
      write(6,6501)LG,XMIN,XMAX,NP,DX,STRENGTH
6501  format('   GAUSS: ',l8  ,' XMIN:',g8.2  ,'  XMAX:',g8.2,/,
     1       '      NP: ',i8  ,'   DX:',g8.2  ,'  STRE:',g8.2)
      write(6,6502)GEXCL,EEXCL,EDIF,LEXCL,NK,kmax
6502  format('    GEXC: ',i8  ,' EEXC:',i8    ,'  EDIF:',g8.2,/,
     1       '    LEXC: ',i8  ,'   NK:',g8.2  ,'  KMAX:',g8.2)
      OH=10.0d0**(-NOH*(N-N7+1)*2)
      write(6,6503)ie0,IWR,NLIM,NB0,BOXF,M1,N7,N,NOH,OH
6503  format('     IE0: ',i8  ,'  IWR:',i8    ,'  NLIM:',i8,/,
     1       '     NB0: ',i8  ,' BOXF:',f8.2  ,'    M1:',i8,/,
     1       '      N7: ',i8  ,'    N:',i8    ,'   NOH:',i8,/,
     1       '      OH: ',g8.2)
      ktcm=0.69507d0*TEMP
c     maximu exc number/Hermite polynomial degree:
      if(nk.gt.nk0)call report('too many k-points')
      if(LEXCL.gt.maxe)call report('Too many excitations required')
      if(maxn.eq.0)maxn=max(EEXCL,GEXCL)+1
      if(kmax.gt.2*NK+1)call report('KMAX > 2NK+1')

      write(6,*)'Ground and excited state oscillator energies:'
      do 111 iq=N7,N
111   write(6,6508)iq,e(iq),es(iq)
6508  format(i4,2f10.1)
      write(6,6045)maxn-1
6045  format(' Max n    :',i5)
      if(maxn.gt.maxb)call report('maxn too high')
      write(6,6009)imaxu,maxu
6009  format(' Max shift:',i5,f15.4)
      box=BOXF*(dsqrt(2.0d0*dble(maxn))*3.0d0*sqrt(2.0d0)+maxu)

c     iniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniini
      do 105 k=1,maxn
      akfac=1.0d0
      do 101 kk=1,k-1
101   akfac=akfac*dble(kk)
105   ank(k)=dsqrt(dsqrt(1.0d0/pi)/2.0d0**(k-1)/akfac)
      
      do 103 ii=1,maxn
      do 103 jj=1,maxn
103   p(jj,ii)=0.0d0
      p(1,1)=1.0d0
      k=1
102   k=k+1
      p(k,1  )=-p(k-1,2)
      do 104 j=2,k-2
104   p(k,j  )=(2.0d0*p(k-1,j-1)-dble(j)*p(k-1,j+1))
      if(k-2.gt.0)p(k,k-1)=2.0d0*p(k-1,k-2)
      p(k,k  )=2.0d0*p(k-1,k-1)
      if(k.lt.maxn)goto 102
c     iniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniini

c
c     transformation of all basis:x(q) .. harmonic oscillator fcion of q
c     xk=int(x(q) exp(-ik.q),q=-inf...inf)
      do 1 k=-nk,nk
      ak=tp*dble(k)/box
      do 1 j=1,maxn
      xkr(j,k+nk+1)=0.0d0
      xki(j,k+nk+1)=0.0d0
      do 1 ii=1,j
c     int(q^n exp(-q^2/2-ik.q)/sqrt(2*pi),q=-inf...inf)
c     air,aii .. real an imag. parts
      call intq(ii-1,ak,air,aii)
      xkr(j,k+nk+1)=xkr(j,k+nk+1)+p(j,ii)*ank(j)*air*stp/box
1     xki(j,k+nk+1)=xki(j,k+nk+1)+p(j,ii)*ank(j)*aii*stp/box
c
c     normalize coeffs
c     do 5 j=1,maxn
c     akfac=0.0d0
c     do 6 k=-nk,nk
c     akfac=akfac+xkr(j,k+nk+1)**2+xki(j,k+nk+1)**2
c     if(akfac.ne.0.0d0)then
c      do 7 k=-nk,nk
c      xki(j,k+nk+1)=xki(j,k+nk+1)/dsqrt(akfac)
c      xkr(j,k+nk+1)=xkr(j,k+nk+1)/dsqrt(akfac)
c     endif
c     continue
c
      write(6,*)'trafo done'
      open(30,file='0.prn')
      open(39,file='n.prn')
      open(40,file='0f.prn')
      open(49,file='nf.prn')
      do 2 iq=-500,500
      q=dble(iq)/1000.d0*box
      xq=0.0d0
      do 3 kk=-nk,nk
      ak=tp*dble(kk)/box
3     xq=xq+xkr(1,kk+nk+1)*cos(ak*q)-xki(1,kk+nk+1)*sin(ak*q)
      write(30,3001)q,xq
3001  format(2f20.3)
      write(40,3001)q,fik(1,q)
      xq=0.0d0
      do 39 kk=-nk,nk
      ak=tp*dble(kk)/box
39    xq=xq+xkr(maxn,kk+nk+1)*cos(ak*(q-maxu))
     1     -xki(maxn,kk+nk+1)*sin(ak*(q-maxu))
      write(39,3001)q,xq
2     write(49,3001)q,fik(maxn,q-maxu)
      close(30)
      close(40)
      close(39)
      close(49)
      write(6,*)' HO transformations done'


c     order the coefficients according to their size:
      if(ifc.eq.9)then
       do 401 j=1,maxn
       do 4012 k=1,2*nk+1
       ykr(j,k)=xkr(j,k)
       yki(j,k)=xki(j,k)
4012   kkr(j,k)=k-nk-1
       do 401 k=1,2*nk+1
       do 401 kk=k+1,2*nk+1
       if(ykr(j,kk)**2+yki(j,kk)**2.gt.ykr(j,k)**2+yki(j,k)**2)then
        yt=ykr(j,kk)
        ykr(j,kk)=ykr(j,k)
        ykr(j,k)=yt
        yt=yki(j,kk)
        yki(j,kk)=yki(j,k)
        yki(j,k)=yt
        kt=kkr(j,kk)
        kkr(j,kk)=kkr(j,k)
        kkr(j,k)=kt
       endif
401    continue
       write(6,*)'Coefficients ordered according to their size'
      endif

      call DTRM(WT,M3,N-N7+1,Determ,INDX,Cwork)
      write(6,6111)Determ,N7,N
6111  format(' Transformation determinant:',f10.3,' N7:',i2,' N:',i3)
      Determ=dsqrt(dabs(Determ))
      if(ltest)then
       write(6,*)
       write(6,6001)
       write(6,6002)(nv(iq)-1,iq=N7,N)
       write(6,6005)
       write(6,6002)(mv(iq)-1,iq=N7,N)
       write(6,6004)
       if(ifc.eq.1)then
        ak=fcf1(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
       else
        if(ifc.eq.2)then
         ak=fcf2(nv,mv,nk,ww,uu,LEXCL,nlim,xkr,xki,N7,N,box,
     1   m1,ai,OH)
        else
         if(ifc.eq.3)then
          ak=fcf3(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
         else
          if(ifc.eq.4)then
           ak=fcf4(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
          else
           if(ifc.eq.5)then
            ak=fcf5(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
           else
            if(ifc.eq.6)ak=fcf6(nv,mv,
     1      nk,wi,uu,xkr,xki,N7,N,box,ai,OH,ie0)
           endif
          endif
         endif
        endif
       endif
       ak=ak*Determ
       ai=ai*Determ
       write(6,6003)ifc,ak,ai
       stop 
      endif
      it=0
      if(ltab)then
       open(9,file='FCF.TAB')
       write(9,*)
       write(9,*)
       write(9,*)
      endif
      if(lspec)then
       do 400 iq=1,NP
400    S(iq)=0.0d0
      endif

      zpeg=0.0d0
      zpee=0.0d0
      do 52 iq=N7,N
      zpeg=zpeg+e(iq)/2.0d0
52    zpee=zpee+es(iq)/2.0d0

      ifct=0.0d0
      if(lanharm)then
c      produce ground states:
       do 503 ig=NG1,NG2
       EN=EAHg(ig)*219470.0d0
       Boltzmann=exp(-EN/ktcm)
       if(Boltzmann.lt.OH)goto 503
       if(iwr.ne.0)write(6,6112)Boltzmann
c      find symmetry:
       KK=0
       DO 1001 IK=1,KGg
1001   IF(IG.GE.ISIg(IK,3).AND.IG.LE.ISIg(IK,4))KK=IK
       I1g=ISIg(KK,1)
       I2g=ISIg(KK,2)
       WRITE(6,302)IG,KK,EN
302    format(' Ground state ',i5,', symmetry ',i2,/,
     1        '       energy ',f10.2,' cm-1')
       CALL LOADC('cfig.scr',IG,CFIG,IG2-IG1+1)

C      Loop over all symmetries
       DO 501 I0=1,KGe
       I1e=ISIe(I0,1)
       I2e=ISIe(I0,2)
       J1e=ISIe(I0,3)
       J2e=ISIe(I0,4)
       NDIM=I2e-I1e+1
       NST=J2e-J1e+1
       IF(I1e.NE.0)THEN
        write(6,*)'Symmetry ',I0,NST,NDIM
        allocate(BIGC(NDIM*NST),akb(NDIM),eeb(NDIM))

        DO 816 IE=J1e,J2e
        ee=EAHe(IE)*219470.0d0
        wcm=8065.5d0*EDIF+ee-en+zpee-zpeg
        wnm=1.0d7/wcm
        wev=wcm/8065.5d0
        if(units.eq.1)then
         om=wev
        else
         if(units.eq.2)then
          om=wnm
         else
          if(units.eq.3)then
           om=wcm
          else
           call report('unknown units')
          endif
         endif
        endif
        ifct=ifct+1.0d0
        eeb(IE)=om
        CALL LOADC('cfie.scr',IE,CFIe,NDIM)
        DO 816 I=1,NDIM
816     BIGC(I+(IE-1)*NDIM)=CFIe(I)

        call fcf8ah(nk,w,u,xkr,xki,N7g,Ng,N7e,Ne,box,
     1  akb,OH,ie0,NST,RSCALE,INg,ISg,ISe,INe,CFIG,NHBIGg,LEXCLg,
     1  NHBIGe,LEXCLe,I1g,I2g,I1e,I2e,CFIE,maxn,ICCg,ICCe)

        do 3041 ii=1,NDIM
        if(dabs(akb(ii)).gt.OH)then
         ak=akb(ii)
         ak=ak*Determ
         om=eeb(ii)
         y=STRENGTH*ak**2*Boltzmann
         if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0)
         it=it+1
         if(ltab)write(9,6011)it,om,y,IG,II
6011     format(i4,f18.9,g17.6,' ',2I8)
        endif
3041    continue
        deallocate(BIGC,akb,eeb)
       ENDIF
501    continue
503    continue
       if(ltab)write(9,90011)
       if(ltab)close(9)
       if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,0.0d0,0.0d0,1)
       write(6,6505)ifct
       write(6,6504)(time()-tim)
       deallocate(S)
       stop
c      end of anharmonic
      endif

      if(ifc.ge.7.and.ifc.le.9)then
       allocate(mvb(NB0*M3),akb(NB0),eeb(NB0))
      endif

c     produce ground states
c     loop start 11111111111111111111111111111111111111111111111111
      do 31 Nnexc=0,GEXCL
      do 411 iex=1,Nnexc
411   sn(iex)=N7
51    DO  19661 I=N7,N
19661 STATEn(I)=0
      do 44 iex=1,Nnexc
44    staten(sn(iex))=staten(sn(iex))+1
      EN=0.0d0
      DO 45 I=N7,N
45    EN=EN+dble(staten(I))*E(I)
      Boltzmann=exp(-EN/ktcm)
      if(Boltzmann.lt.OH)goto 444
      if(iwr.ne.0)write(6,6112)Boltzmann
6112  format(' exp(-E/kt):',f10.6)
c     do 45 iex=1,Nnexc
c5    if(staten(sn(iex)).gt.m1)goto 444
      mm=0
c     produce excited states
c     22222222222222222222222222222222222222222222222222222222222222
      do 311 Nexc=0,EEXCL
      do 4111 iex=1,Nexc
4111  sm(iex)=N7
511   DO  661 I=N7,N
661   STATE(I)=0
      do 441 iex=1,Nexc
441   state(sm(iex))=state(sm(iex))+1
      ee=0.0d0
      DO 46 I=N7,N
46    ee=ee+dble(state(I))*ES(I)
c     do 45 iex=1,Nnexc
c5    if(staten(sn(iex)).gt.m1)goto 344
      do 47 iq=N7,N
      nv(iq)=staten(iq)+1
47    mv(iq)=state(iq) +1

      wcm=8065.5d0*EDIF+ee-en+zpee-zpeg
      wnm=1.0d7/wcm
      wev=wcm/8065.5d0
      if(units.eq.1)then
       om=wev
      else
       if(units.eq.2)then
        om=wnm
       else
        if(units.eq.3)then
         om=wcm
        else
         call report('unknown units')
        endif
       endif
      endif
      ifct=ifct+1.0d0
      if(ifc.eq.1)then
       ak=fcf1(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
      else
       if(ifc.eq.2)then
        ak=fcf2(nv,mv,nk,ww,uu,LEXCL,nlim,xkr,xki,N7,N,box,
     1  m1,ai,OH)
       else
        if(ifc.eq.3)then
         ak=fcf3(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
        else
         if(ifc.eq.4)then
          ak=fcf4(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
         else
          if(ifc.eq.5)then
           ak=fcf5(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH)
          else
      if(ifc.eq.6)then
       ak=fcf6(nv,mv,nk,wi,uu,xkr,xki,N7,N,box,ai,OH,ie0)
      else
       if(ifc.eq.7.or.ifc.eq.8.or.ifc.eq.9)then
        if(GEXCL.ne.0)call report('ifc 7, 8 and 9 only for ground 0')
        mm=mm+1
        eeb(mm)=om
        do 301 iq=N7,N
301     mvb(mm+(iq-1)*NB0)=mv(iq)
        if(mm.eq.NB0)then
         if(ifc.eq.7)then
          call fcf7(nv,mvb,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,
     1    m1,akb,OH,mm,NB0)
         else
          if(ifc.eq.8)then
           call fcf8(nv,mvb,nk,wi,uu,xkr,xki,N7,N,box,
     1     akb,OH,ie0,mm,NB0,RSCALE)
          else
           call fcf9(nv,mvb,nk,kmax,wi,uu,LEXCL,xkr,xki,N7,N,box,
     1     akb,OH,mm,ykr,yki,kkr,NB0)
          endif
         endif
         do 303 ii=1,mm
         if(dabs(akb(ii)).gt.OH)then
          ak=akb(ii)
          ak=ak*Determ
          om=eeb(ii)
          y=STRENGTH*ak**2*Boltzmann
          if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0)
          if(ltab)then
           it=it+1
           write(9,900)it,om,y
           write(9,6001)
6001       format('<',$) 
           write(9,6002)(nv(iq)-1,iq=N7,N)
6002       format(i1,$) 
           write(9,6005)
6005       format('|',$) 
           write(9,6002)(mvb(ii+(iq-1)*NB0)-1,iq=N7,N)
           write(9,6004)
6004       format('>') 
          endif
         endif
303      continue
         mm=0
        endif
       else
        call report('unknown IFC')
       endif
      endif
          endif
         endif
        endif
       endif
      endif
      if(ifc.ne.7.and.ifc.ne.8.and.ifc.ne.9.and.dabs(ak).gt.OH)then
       ak=ak*Determ
       ai=ai*Determ
       if(iwr.ne.0)write(6,6003)ifc,ak,ai
6003   format(' FCF',i1,' = ',2f10.6)
       if(ltest)stop
       y=STRENGTH*ak**2*Boltzmann
       if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0)
       it=it+1
       if(ltab)then
        write(9,900)it,om,y
900     format(i4,f18.9,g17.6,' ',$)
        write(9,6001)
        write(9,6002)(nv(iq)-1,iq=N7,N)
        write(9,6005)
        write(9,6002)(mv(iq)-1,iq=N7,N)
        write(9,6004)
       endif
      endif

c344   continue
      if(Nexc.eq.0)goto 311
      do 811 ic=Nexc,1,-1
811   if(sm(ic).lt.N)goto 911
      goto 311
911   do 1011 i=ic+1,Nexc
1011  sm(i)=sm(ic)+1
      sm(ic)=sm(ic)+1
      goto 511
311   continue
c     loop end 222222222222222222222222222222222222222222222222222222 
c    
      if(ifc.eq.7.or.ifc.eq.8.or.ifc.eq.9)then
       if(ifc.eq.7)then
        call fcf7(nv,mvb,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,
     1  m1,akb,OH,mm,NB0)
       else
        if(ifc.eq.8)then
         call fcf8(nv,mvb,nk,wi,uu,xkr,xki,N7,N,box,
     1   akb,OH,ie0,mm,NB0,RSCALE)
        else
         call fcf9(nv,mvb,nk,kmax,wi,uu,LEXCL,xkr,xki,N7,N,box,
     1   akb,OH,mm,ykr,yki,kkr,NB0)
        endif
       endif
       do 304 ii=1,mm
       if(dabs(akb(ii)).gt.OH)then
        ak=akb(ii)
        ak=ak*Determ
        om=eeb(ii)
        y=STRENGTH*ak**2*Boltzmann
        if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0)
        it=it+1
        if(ltab)then
         write(9,900)it,om,y
         write(9,6001)
         write(9,6002)(nv(iq)-1,iq=N7,N)
         write(9,6005)
         write(9,6002)(mvb(ii+(iq-1)*NB0)-1,iq=N7,N)
         write(9,6004)
        endif
       endif
304    continue
      endif

444   if(Nnexc.eq.0)goto 31
      do 81 ic=Nnexc,1,-1
81    if(sn(ic).lt.N)goto 91
      goto 31
91    do 402 i=ic+1,Nnexc
402   sn(i)=sn(ic)+1
      sn(ic)=sn(ic)+1
      goto 51
31    continue
c     loop end 11111111111111111111111111111111111111111111111111

      if(ifc.ge.7.and.ifc.le.9)then
       deallocate(mvb,akb,eeb)
      endif

      if(ltab)write(9,90011)
      if(ltab)close(9)
90011 format(60(1h-))

      if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,0.0d0,0.0d0,1)

      write(6,6505)ifct
6505  format(f20.0,' factors')
      write(6,6504)(time()-tim)
6504  format(' Time: ',i8,' sec. ')
      deallocate(S)
      stop
      end
c     =======================================================
      subroutine tabprnfs(S,wsmin,wsmax,LG,npoint,hw,wo,x,ic)
      implicit none
      integer*4 ic,npoint,i
      real*8 wsmax,wsmin,s(*),x,dw,dd,hw,xx,wo,sf
      logical LG

      dw=(wsmax-wsmin)/dble(npoint-1)
      if(ic.eq.0)then
c      disperse value over spectrum 
       if(LG)then
        dd=hw/2.0d0/dsqrt(log(2.0d0))
       else
        dd=hw/2.0d0
       endif
       do 4 i=1,npoint
       xx=wsmin+dw*dble(i-1)
4      s(i)=s(i)+sf(dd,LG,xx,wo)*x
      else
c      write spectrum
       call wrspec2('S.PRN',npoint,wsmin,dw,s)
      endif
      return
      end
c     ==================================================================
      subroutine wrspec2(sn,nw,wmin,dw,spec)
      implicit none
      character*(*) sn
      integer*4 nw,iw
      real*8 wmin,dw,spec(*),w
      open(8,file=sn)
      do 1 iw=1,nw
      w=wmin+dble(iw-1)*dw
1     write(8,805)w,spec(iw)
805   format(f18.7,g15.6)
      close(8)
      write(6,*)sn,' written'
      return
      end
c     ==================================================================
      function sf(d,g,x,x0)
      implicit none
      logical g
      real*8 pi,spi,sf,d,x,x0,dd
      pi=4.0d0*atan(1.0d0)
      spi=dsqrt(pi)
      dd=((x-x0)/d)**2
c     write(6,*)d,x,x0,g
      if(g)then
       if(dd.lt.20.0d0)then
        sf=exp(-dd)/d/spi
       else
        sf=0.0d0
       endif
      else
       if(dd.lt.1.0d20)then
        sf=1.0d0/d/(dd+1.0d0)/pi
       else
        sf=0.0d0
       endif
      endif
      return
      end
c     ------------------------------------------------------------
      function fcf1(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box,
     1m1,fcsumi,OH)
c     via delta-f, standard
      implicit none
      integer*4 maxb,nk0,M3,maxe
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL,ii
      real*8 w(M3,M3),u(*),fcf1,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,prl,pil,box,fim,pi,tp,aw(M3),
     1ptr,pti,OH,t,xkrn,xkin,ak,air,aii,ank,p,stp,e
      common/ho/p(maxb,maxb),ank(maxb),e(M3)
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      stp=dsqrt(tp)
      fcsumr=0.0d0
      fcsumi=0.0d0
c     distribute k-"excitations" on electronic excited coeff
C     "Ground state" zero coeffs. - always:
c     a.b = (ReaReb-ImaImb) + i (ReaImb+ImaReb):
c     ps=xs1(0)*xs2(0)...xsN(0)
      prs= xkr(mv(N7),nk+1)
      pis= xki(mv(N7),nk+1)
      do 1 i=N7+1,N
      t  = prs*xkr(mv(i),nk+1)-pis*xki(mv(i),nk+1)
      pis= prs*xki(mv(i),nk+1)+pis*xkr(mv(i),nk+1)
      prs=t
1     if(abs(prs)+abs(pis).lt.OH)goto 1001
c     pl=xl1(0)*xl2(0)...xlN(0)
      prl= xkr(nv(N7),nk+1)
      pil= xki(nv(N7),nk+1)
      do 2 i=N7+1,N
      t  = prl*xkr(nv(i),nk+1)-pil*xki(nv(i),nk+1)
      pil= prl*xki(nv(i),nk+1)+pil*xkr(nv(i),nk+1)
      prl=t
2     if(abs(prl)+abs(pil).lt.OH)goto 1001
c     ps* x pl:
c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*prl+pis*pil)*box
      fcsumi=fcsumi+(prs*pil-pis*prl)*box
c
c     distribute LEXCL onto the excited k-expansions N7...N
c     loop start 22222222222222222222222222222222222222222222222222
1001  do 3 Nexc=1,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the phase and product exp(if)xs...
      fim=0.0d0
      do 10 I=N7,N
10    fim=fim+u(I)*kstate(I)*tp/box
      prs=cos(fim)
      pis=sin(fim)
      do 11 i=N7,N
      t  = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i))
      pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c     write(6,*)'ks',kstate(1),prs,pis
c
c     calculate the aw products (skip 2 pi /L multiplier)
      do 12 i=N7,N
      aw(i)=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    aw(i)=aw(i)+dble(kstate(j))*w(j,i)*tp/box
c
c     calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)...
      ptr=1.0d0
      pti=0.0d0
      do 111 i=N7,N
      ak=aw(i)
      xkrn=0.0d0
      xkin=0.0d0
      do 6 ii=1,nv(i)
      call intq(ii-1,ak,air,aii)
      xkrn=xkrn+p(nv(i),ii)*ank(nv(i))*air*stp
6     xkin=xkin+p(nv(i),ii)*ank(nv(i))*aii*stp
c     write(6,*)xkrn,xkin
c     write(6,*)xkr(nv(i),nk+1+kstate(i)),xki(nv(i),nk+1+kstate(i))
      t  = ptr*xkrn-pti*xkin
      pti= ptr*xkin+pti*xkrn
      ptr=t
111   if(abs(ptr)+abs(pti).lt.OH)goto 333
c     write(6,*)'kl',kstate(1),ptr,pti

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*ptr+pis*pti)
      fcsumi=fcsumi+(prs*pti-pis*ptr)

c     find index to be changed
333   do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      fcf1=fcsumr
      return
      end
c     ------------------------------------------------------------
      function fcf2(nv,mv,nk,w,u,LEXCL,nlim,xkr,xki,N7,N,box,
     1m1,fcsumi,OH)
c     left sin ... very slow, results virtually sae as fcf1
c     not recommended for normal computation
      implicit none
      integer*4 maxb,nk0,M3,maxe
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),nlim,N7,i,nk,N,Nexc,Nnexc,iex,m1,j,
     1state(M3),staten(M3),si(maxe),sn(maxe),kstate(M3),kstaten(M3),
     1naw(M3),ic,LEXCL
      real*8 w(M3,M3),u(*),fcf2,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,prl,pil,box,fim,pi,tp,aw(M3),
     1dk,sf,ptr,pti,OH,t
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     nlim .. take ground ks only within this limit around excited ks
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      fcsumr=0.0d0
      fcsumi=0.0d0
c     distribute k-"excitations" on electronic excited coeff
C     "Ground state" zero coeffs. - always:
c     a.b = (ReaReb-ImaImb) + i (ReaImb+ImaReb):
c     ps=xs1(0)*xs2(0)...xsN(0)
      prs= xkr(mv(N7),nk+1)
      pis= xki(mv(N7),nk+1)
      do 1 i=N7+1,N
      t  = prs*xkr(mv(i),nk+1)-pis*xki(mv(i),nk+1)
      pis= prs*xki(mv(i),nk+1)+pis*xkr(mv(i),nk+1)
      prs=t
1     if(abs(prs)+abs(pis).lt.OH)goto 1001
c     pl=xl1(0)*xl2(0)...xlN(0)
      prl= xkr(nv(N7),nk+1)
      pil= xki(nv(N7),nk+1)
      do 2 i=N7+1,N
      t  = prl*xkr(nv(i),nk+1)-pil*xki(nv(i),nk+1)
      pil= prl*xki(nv(i),nk+1)+pil*xkr(nv(i),nk+1)
      prl=t
2     if(abs(prl)+abs(pil).lt.OH)goto 1001
c     ps* x pl:
c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*prl+pis*pil)*box
      fcsumi=fcsumi+(prs*pil-pis*prl)*box
c
c     distribute LEXCL onto the excited k-expansions N7...N
c     loop start 22222222222222222222222222222222222222222222222222
1001  do 3 Nexc=1,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the phase and product exp(if)xs...
      fim=0.0d0
      do 10 I=N7,N
10    fim=fim+u(I)*kstate(I)*tp/box
      prs=cos(fim)
      pis=sin(fim)
      do 11 i=N7,N
      t  = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i))
      pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c     write(6,*)'ks',kstate(1),prs,pis
c
c     calculate the aw products (skip 2 pi /L multiplier)
      do 12 i=N7,N
      aw(i)=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    aw(i)=aw(i)+dble(kstate(j))*w(j,i)
c     closest integer value:
      do 121 i=N7,N
121   naw(i)=NINT(aw(i))
c
c     take n-k vectors in the vicinity of aw
c     generate excitations on the n(ground state) site
c     loop start 11111111111111111111111111111111111111111111111111
      do 31 Nnexc=1,LEXCL
      do 411 iex=1,Nnexc
411   sn(iex)=N7
51    DO  19661 I=N7,N
19661 STATEn(I)=0
      do 44 iex=1,Nnexc
44    staten(sn(iex))=staten(sn(iex))+1
      do 45 iex=1,Nnexc
45    if(staten(sn(iex)).gt.m1)goto 444

      do 421 I=N7,N
      if(mod(staten(I),2).eq.0)then
       kstaten(I)=-staten(I)/2
      else
       kstaten(I)=(staten(I)+1)/2
      endif
c     take only states that are close to aw
      if(abs(kstaten(i)).gt.nk)goto 444
421   if(iabs(kstaten(I)-naw(I)).gt.nlim)goto 444
c
c     calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)...
      ptr=1.0d0
      pti=0.0d0
      do 111 i=N7,N
      t  = ptr*xkr(nv(i),nk+1+kstaten(i))-pti*xki(nv(i),nk+1+kstaten(i))
      pti= ptr*xki(nv(i),nk+1+kstaten(i))+pti*xkr(nv(i),nk+1+kstaten(i))
      ptr=t
      if(abs(ptr)+abs(pti).lt.OH)goto 444
      dk=dble(kstaten(i))-aw(i)
      if(dabs(dk).gt.OH)then
       sf=2.0d0*sin(pi*dk)/(tp*dk/box)
      else
       sf=box
      endif
      ptr=ptr*sf
111   pti=pti*sf
c     write(6,*)'kl',kstaten(1),ptr,pti

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*ptr+pis*pti)
      fcsumi=fcsumi+(prs*pti-pis*ptr)
c
444   do 81 ic=Nnexc,1,-1
81    if(sn(ic).lt.N)goto 91
      goto 31
91    do 101 i=ic+1,Nnexc
101   sn(i)=sn(ic)+1
      sn(ic)=sn(ic)+1
      goto 51
31    continue
c     loop end 11111111111111111111111111111111111111111111111111

c     find index to be changed
333   do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      fcf2=fcsumr
      return
      end
c     ------------------------------------------------------------
c     ------------------------------------------------------------
      function fcf3(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box,
     1m1,fcsumi,OH)
c     ground state first, otherwise same as fcf1
      implicit none
      integer*4 maxb,nk0,M3,maxe
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL,ii
      real*8 w(M3,M3),u(*),fcf3,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1ptr,pti,OH,t,xkrn,xkin,ak,air,aii,ank,p,stp,e
      common/ho/p(maxb,maxb),ank(maxb),e(M3)
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      stp=dsqrt(tp)
      fcsumr=0.0d0
      fcsumi=0.0d0

c     distribute k-"excitations" on electronic excited coeff
c     start with ground as it may need less k-points
      do 3 Nexc=0,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the product xl...
      prs=1.0d0
      pis=0.0d0
      do 11 i=N7,N
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      ptr=1.0d0
      pti=0.0d0
      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
      ak=ak*tp/box
c     calculate the phase:
      fim=fim+u(i)*ak
c
c     calculate Fourier expansion of excited at ak:
      xkrn=0.0d0
      xkin=0.0d0
      do 6 ii=1,mv(i)
      call intq(ii-1,ak,air,aii)
      xkrn=xkrn+p(mv(i),ii)*ank(mv(i))*air*stp
6     xkin=xkin+p(mv(i),ii)*ank(mv(i))*aii*stp
      t  = ptr*xkrn-pti*xkin
      pti= ptr*xkin+pti*xkrn
      ptr=t
111   if(abs(ptr)+abs(pti).lt.OH)goto 333

      t  = ptr*cos(fim)-pti*sin(fim)
      pti= ptr*sin(fim)+pti*cos(fim)
      ptr=t

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*ptr+pis*pti)
      fcsumi=fcsumi+(prs*pti-pis*ptr)

c     find index to be changed
333   if(Nexc.eq.0)goto 3
      do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      fcf3=fcsumr
      return
      end
c     ------------------------------------------------------------
      function fcf4(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box,
     1m1,fcsumi,OH)
c     fcf1, with an interpolation for k-coefficients
      implicit none
      integer*4 maxb,nk0,M3,maxe
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL,km
      real*8 w(M3,M3),u(*),fcf4,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      fcsumr=0.0d0
      fcsumi=0.0d0
c
c     distribute LEXCL onto the excited k-expansions N7...N
c     loop start 22222222222222222222222222222222222222222222222222
      do 3 Nexc=0,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the phase and product exp(if)xs...
      fim=0.0d0
      do 10 I=N7,N
10    fim=fim+u(I)*kstate(I)
      fim=fim*tp/box
      prs=cos(fim)
      pis=sin(fim)
      do 11 i=N7,N
      t  = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i))
      pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)...
      ptr=1.0d0
      pti=0.0d0
      do 111 i=N7,N
      ak=0.0d0
      do 12 j=N7,N
c     w-first index (j) excited
12    ak=ak+dble(kstate(j))*w(j,i)
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 333
      if(abs(km+1).gt.nk)goto 333
      dk=ak-dble(km)
      dko=1.0d0-dk
c
c     Fourier coefficients by interpolation:
      xkrn=xkr(nv(i),nk+1+km)*dko+dk*xkr(nv(i),nk+1+km+1)
      xkin=xki(nv(i),nk+1+km)*dko+dk*xki(nv(i),nk+1+km+1)

      t  = ptr*xkrn-pti*xkin
      pti= ptr*xkin+pti*xkrn
      ptr=t
111   if(abs(ptr)+abs(pti).lt.OH)goto 333

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*ptr+pis*pti)
      fcsumi=fcsumi+(prs*pti-pis*ptr)

c     find index to be changed
333   if(Nexc.eq.0)goto 3
      do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      fcf4=fcsumr*box**(N-N7+1)
      return
      end
c     ------------------------------------------------------------
      function fcf5(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box,
     1m1,fcsumi,OH)
c     fcf3 with interpolation for k-coefficients
      implicit none
      integer*4 maxb,nk0,M3,maxe
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,km,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL
      real*8 w(M3,M3),u(*),fcf5,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      fcsumr=0.0d0
      fcsumi=0.0d0

c     distribute k-"excitations" on electronic excited coeff
c     start with ground as it may need less k-points
      do 3 Nexc=0,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the product xl...
      prs=1.0d0
      pis=0.0d0
      do 11 i=N7,N
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      ptr=1.0d0
      pti=0.0d0
      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 333
      if(abs(km+1).gt.nk)goto 333
      dk=ak-dble(km)
      dko=1.0d0-dk

c     calculate the phase:
      fim=fim+u(i)*ak*tp/box
c     write(6,*)fim,ak,u(i)
c
c     Fourier coefficients by interpolation:
      xkrn=xkr(mv(i),nk+1+km)*dko+dk*xkr(mv(i),nk+1+km+1)
      xkin=xki(mv(i),nk+1+km)*dko+dk*xki(mv(i),nk+1+km+1)
c
      t  = ptr*xkrn-pti*xkin
      pti= ptr*xkin+pti*xkrn
      ptr=t
111   if(abs(ptr)+abs(pti).lt.OH)goto 333
c     read(5,*)ic


      t  = ptr*cos(fim)-pti*sin(fim)
      pti= ptr*sin(fim)+pti*cos(fim)
      ptr=t

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=fcsumr+(prs*ptr+pis*pti)
      fcsumi=fcsumi+(prs*pti-pis*ptr)

c     find index to be changed
333   if(Nexc.eq.0)goto 3
      do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      fcf5=fcsumr*box**(N-N7+1)
      return
      end
c     ------------------------------------------------------------
      function fcf6(nv,mv,nk,w,u,xkr,xki,N7,N,box,
     1fcsumi,OH,ie0)
c     integration by random number generation
      implicit none
      integer*4 maxb,nk0,M3
      parameter (maxb=100,M3=100,nk0=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,j,km,
     1kstate(M3),k,
     1ie0
      real*8 w(M3,M3),u(*),fcf6,fcsumr,fcsumi,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko,
     1it,ip,im,amax,amaxi,frandom,y
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      integer*4 ix,iy,iz,itc
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      fcsumr=0.0d0
      fcsumi=0.0d0
      amax=1.0d0
      do 200 i=N7,N
      amaxi=0.0d0
      do 201 k=-nk,nk
      if(dabs(xkr(nv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(nv(i),nk+1+k))
201   if(dabs(xkr(mv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(mv(i),nk+1+k))
200   amax=amax*amaxi**2

      it=0
      ip=0
      im=0

      do 4000 itc=1,ie0

      do 43 I=N7,N
43    kstate(I)=int(-nk+int(dble(nk+nk)*frandom()))
      y=-amax+(amax+amax)*frandom()
c
      prs=1.0d0
      pis=0.0d0
      ptr=1.0d0
      pti=0.0d0

c     calculate the product xl...
      do 11 i=N7,N
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 4000
      if(abs(km+1).gt.nk)goto 4000
      dk=ak-dble(km)
      dko=1.0d0-dk

c     calculate the phase:
      fim=fim+u(i)*ak*tp/box
c
c     Fourier coefficients by interpolation:
      xkrn=xkr(mv(i),nk+1+km)*dko+dk*xkr(mv(i),nk+1+km+1)
      xkin=xki(mv(i),nk+1+km)*dko+dk*xki(mv(i),nk+1+km+1)
c
      t  = ptr*xkrn-pti*xkin
      pti= ptr*xkin+pti*xkrn
      ptr=t
111   if(abs(ptr)+abs(pti).lt.OH)goto 333

333   continue

      t  = ptr*cos(fim)-pti*sin(fim)
      pti= ptr*sin(fim)+pti*cos(fim)
      ptr=t

c     a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb):
      fcsumr=prs*ptr+pis*pti
      if(y.lt.fcsumr.and.y.gt.0)ip=ip+1.0d0
      if(y.gt.fcsumr.and.y.lt.0)im=im+1.0d0
      it=it+1.0d0

4000  continue

      fcf6=(ip-im)/it*(amax+amax)*dble(nk)**(2*(N-N7+1))

      return
      end
c     ------------------------------------------------------------
      subroutine fcf7(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box,
     1m1,fcsumr,OH,mm,NB0)
c     fcf3 with interpolation for k-coefficients
c     and with batch of m mm vectors
      implicit none
      integer*4 maxb,nk0,M3,maxe,NB0
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,km,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL,mm,ii,im
      real*8 w(M3,M3),u(*),fcsumr(*),xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1OH,t,xkrn,xkin,ak,dk,dko
      real*8,allocatable,dimension(:)::ptr,pti
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      allocate(ptr(mm),pti(mm))
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi

      do 1 i=1,mm
1     fcsumr(i)=0.0d0

c     distribute k-"excitations" on electronic excited coeff
c     start with ground as it may need less k-points
      do 3 Nexc=0,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.m1)goto 333
c
c     odd excitation - positive (N+1)/2
c     even excitation - neagtive N/2
      do 43 I=N7,N
      if(mod(state(I),2).eq.0)then
       kstate(I)=-state(I)/2
      else
       kstate(I)=(state(I)+1)/2
      endif
      if(abs(kstate(i)).gt.nk)goto 333
43    continue
c
c     calculate the product xl...ground state the same for all batch
      prs=box**(N-N7+1)
      pis=0.0d0
      do 11 i=N7,N
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      do 7 ii=1,mm
      ptr(ii)=1.0d0
7     pti(ii)=0.0d0

      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
c     contribution to the phase:
      fim=fim+u(i)*ak*tp/box
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 333
      if(abs(km+1).gt.nk)goto 333
      dk=ak-dble(km)
      dko=1.0d0-dk
      do 111 ii=1,mm
      im=mv(ii+(i-1)*NB0)
c     coefficients by interpolation:
      xkrn=xkr(im,nk+1+km)*dko+dk*xkr(im,nk+1+km+1)
      xkin=xki(im,nk+1+km)*dko+dk*xki(im,nk+1+km+1)
      t      = ptr(ii)*xkrn-pti(ii)*xkin
      pti(ii)= ptr(ii)*xkin+pti(ii)*xkrn
111   ptr(ii)= t

      do 6 ii=1,mm
      t      = ptr(ii)*cos(fim)-pti(ii)*sin(fim)
      pti(ii)= ptr(ii)*sin(fim)+pti(ii)*cos(fim)
      ptr(ii)= t
6     fcsumr(ii)=fcsumr(ii)+(prs*ptr(ii)+pis*pti(ii))

c     find index to be changed
333   if(Nexc.eq.0)goto 3
      do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      deallocate(ptr,pti)
      return
      end
c     ------------------------------------------------------------
      subroutine fcf8(nv,mv,nk,w,u,xkr,xki,N7,N,box,
     1fcsumr,OH,ie0,mm,NB0,RSCALE)
c     integration by random number generation
c     and batch of m  mm vectors
      implicit none
      integer*4 maxb,nk0,M3,NB0
      parameter (maxb=100,M3=100,nk0=200)
      integer*4 nv(*),mv(*),N7,i,nk,N,j,km,kstate(M3),ii,k,imm,
     1ie0,mm,iu
      real*8 w(M3,M3),u(*),fcsumr(*),sumr,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,OH,t,xkrn,xkin,ak,dk,dko,
     1amaxi,frandom,RSCALE
      real*8,allocatable::it(:),ip(:),im(:),amax(:),y(:),ptr(:),pti(:)
      integer*4 ix,iy,iz
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      allocate(it(NB0),ip(NB0),im(NB0),amax(NB0),y(NB0),
     1ptr(NB0),pti(NB0))
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      do 1 i=1,mm
1     fcsumr(i)=0.0d0
      do 200 ii=1,mm
      amax(ii)=1.0d0
      do 200 i=N7,N
      amaxi=0.0d0
      imm=mv(ii+(i-1)*NB0)
      do 201 k=-nk,nk
      if(dabs(xkr(nv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(nv(i),nk+1+k))
201   if(dabs(xkr(imm,nk+1+k)).gt.amaxi)amaxi=dabs(xkr(imm,nk+1+k))
200   amax(ii)=amax(ii)*amaxi**2*RSCALE

      do 6 ii=1,mm
      it(ii)=0.0d0
      ip(ii)=0.0d0
6     im(ii)=0.0d0

      do 4000 iu=1,ie0

      do 43 I=N7,N
43    kstate(I)=int(-nk+int(dble(nk+nk)*frandom()))
      do 8 ii=1,mm
8     y(ii)=-amax(ii)+(amax(ii)+amax(ii))*frandom()
c
      prs=1.0d0
      pis=0.0d0
      do 2 ii=1,mm
      ptr(ii)=1.0d0
2     pti(ii)=0.0d0

c     calculate the product xl...
      do 11 i=N7,N
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 4000
      if(abs(km+1).gt.nk)goto 4000
      dk=ak-dble(km)
      dko=1.0d0-dk

c     calculate the phase:
      fim=fim+u(i)*ak*tp/box
c
      do 111 ii=1,mm
      imm=mv(ii+(i-1)*NB0)
      xkrn=xkr(imm,nk+1+km)*dko+dk*xkr(imm,nk+1+km+1)
      xkin=xki(imm,nk+1+km)*dko+dk*xki(imm,nk+1+km+1)
      t      =ptr(ii)*xkrn-pti(ii)*xkin
      pti(ii)=ptr(ii)*xkin+pti(ii)*xkrn
111   ptr(ii)=t

333   continue

      do 4 ii=1,mm
      t      = ptr(ii)*cos(fim)-pti(ii)*sin(fim)
      pti(ii)= ptr(ii)*sin(fim)+pti(ii)*cos(fim)
      ptr(ii)= t
      sumr=prs*ptr(ii)+pis*pti(ii)
      if(y(ii).lt.sumr.and.y(ii).gt.0)ip(ii)=ip(ii)+1.0d0
      if(y(ii).gt.sumr.and.y(ii).lt.0)im(ii)=im(ii)+1.0d0
4     it(ii)=it(ii)+1.0d0

4000  continue

      do 5 ii=1,mm
5     fcsumr(ii)=
     1(ip(ii)-im(ii))/it(ii)*amax(ii)*2.0d0*dble(nk)**(2*(N-N7+1))
      deallocate(it,ip,im,amax,y,ptr,pti)

      return
      end
c     ------------------------------------------------------------
      subroutine fcf8ah(nk,w,u,xkr,xki,N7g,Ng,N7e,Ne,box,
     1fcsumr,OH,ie0,NST,RSCALE,INg,ISg,ISe,INe,CFIG,NHBIGg,LEXCLg,
     1NHBIGe,LEXCLe,I1g,I2g,I1e,I2e,CFIE,maxn,ICCg,ICCe)
c
c     integration by random number generation
c     NDIM excited AH states
c     for mm anharmonic excited states
c
c     nk - maximum k expansion
c     w u normal mode transformation matrices
c     xkr xki HO Fourier coefficients
c     N7g,Ng,N7e,Ne   ground and upper normal mode limits
c     box box size
c     fsumr the FCF buffer
c     OH precision limit
c     ie0  random number cycle generation limit
c     NST  number of excited states in the buffer
c     RSCALE scale the integration space
c     INg,ISg,INe,ISe,ICCg,ICCe ground HO state definition
c     CFIG ground state eigenvector
c     CFIE excited state eigenvectors (different dimension than CFIG!)
c     NHBIGg,LEXClg,NHBIGe,LEXCle  state dimensions
c     I1g,I2g,I1e,I2e interval of groudn and excited HO functions
c     maxn maximum HO excitation
      implicit none
      integer*4 maxb,nk0,M3
      parameter (maxb=100,M3=100,nk0=200)
      integer*4 nv(M3),mv(M3),N7g,i,nk,Ng,j,km,kstate(M3),ii,k,imm,
     1ie0,iu,NST,NHBIGg,LEXCLg,NHBIGe,LEXCLe,maxn,jj,Ne,N7e,NDIM,
     1INg(NHBIGg,LEXCLg),ISg(NHBIGg),INe(NHBIGe,LEXCLe),ISe(NHBIGe),
     2ICCg(NHBIGg,LEXCLg),ICCe(NHBIGe,LEXCLe),I1g,I2g,I1e,I2e
      real*8 w(M3,M3),u(*),fcsumr(*),sumr,xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,OH,t,xkrn,xkin,ak,dk,dko,
     1amaxi,frandom,RSCALE,amax,kr,ki,CFIG(*),CFIE(*),ptr,pti
      real*8,allocatable::it(:),ip(:),im(:),y(:),ker(:),kei(:)
      integer*4 ix,iy,iz
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
      NDIM=I2e-I1e+1
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      allocate(it(NST),ip(NST),im(NST),y(NST),ker(NST),kei(NST))
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi
      do 1 i=1,NST
1     fcsumr(i)=0.0d0
      amax=1.0d0
      do 200 i=N7g,Ng
      amaxi=0.0d0
      do 201 k=-nk,nk
      do 201 ii=1,maxn
201   if(dabs(xkr(ii,nk+1+k)).gt.amaxi)amaxi=dabs(xkr(ii,nk+1+k))
200   amax=amax*amaxi**2*RSCALE

      do 6 ii=1,NST
      it(ii)=0.0d0
      ip(ii)=0.0d0
6     im(ii)=0.0d0

      do 4000 iu=1,ie0

      do 43 I=N7g,Ng
43    kstate(I)=int(-nk+int(dble(nk+nk)*frandom()))
      do 8 ii=1,NST
8     y(ii)=-amax+(amax+amax)*frandom()
c
c     calculate ground state coeff of k=vector
      kr=0.0d0
      ki=0.0d0
      do 501 J=I1g,I2g

c     retrieve state J:
      DO 808 ii=N7g,Ng
808   nv(ii)=1
      DO 809 ii=1,ISg(J)
809   nv(ICCg(J,ii))=INg(J,ii)+1

      prs=1.0d0
      pis=0.0d0

c     calculate the product xl...
      do 11 i=N7g,Ng
      t  = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i))
      pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 501
      kr=kr+CFIG(J-I1G+1)*prs
      ki=ki+CFIG(J-I1G+1)*pis
501   continue
c
      do 2 ii=1,NST
      ker(ii)=0.0d0
2     kei(ii)=0.0d0

      do 502 J=I1e,I2e
c     retrieve state J:
      DO 818 ii=N7e,Ne
818   mv(ii)=1
      DO 819 ii=1,ISe(J)
819   mv(ICCe(J,ii))=INe(J,ii)+1
c     calculate the product xs...
      fim=0.0d0
      ptr=1.0d0
      pti=0.0d0
      do 111 i=N7e,Ne
c     calculate corresponding ak:
      ak=0.0d0
      do 12 jj=N7g,Ng
c     w-first index excited
12    ak=ak+dble(kstate(jj))*w(i,jj)
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 4000
      if(abs(km+1).gt.nk)goto 4000
      dk=ak-dble(km)
      dko=1.0d0-dk

c     calculate the phase:
      fim=fim+u(i)*ak*tp/box
c
      xkrn=xkr(mv(i),nk+1+km)*dko+dk*xkr(mv(i),nk+1+km+1)
      xkin=xki(mv(i),nk+1+km)*dko+dk*xki(mv(i),nk+1+km+1)
      t  =ptr*xkrn-pti*xkin
      pti=ptr*xkin+pti*xkrn
111   ptr=t

      do 503 ii=1,NST
      ker(ii)=ker(ii)+CFIE(J+(ii-1)*NDIM)*ptr
503   kei(ii)=kei(ii)+CFIE(J+(ii-1)*NDIM)*pti

502   continue

      do 4 ii=1,NST
      t      = ker(ii)*cos(fim)-kei(ii)*sin(fim)
      kei(ii)= ker(ii)*sin(fim)+kei(ii)*cos(fim)
      ker(ii)= t
      sumr=prs*ker(ii)+pis*kei(ii)
      if(y(ii).lt.sumr.and.y(ii).gt.0)ip(ii)=ip(ii)+1.0d0
      if(y(ii).gt.sumr.and.y(ii).lt.0)im(ii)=im(ii)+1.0d0
4     it(ii)=it(ii)+1.0d0

4000  continue

      do 5 ii=1,NST
5     fcsumr(ii)=
     1(ip(ii)-im(ii))/it(ii)*amax*2.0d0*dble(nk)**(2*(Ng-N7g+1))
      deallocate(it,ip,im,y,ker,kei)

      return
      end
c     ------------------------------------------------------------
      subroutine fcf9(nv,mv,nk,kmax,w,u,LEXCL,xkr,xki,N7,N,box,
     1fcsumr,OH,mm,ykr,yki,kkr,NB0)
c     fcf3 with interpolation for k-coefficients
c     and with batch of m mm vectors
c     and for largest coeffs only
      implicit none
      integer*4 maxb,nk0,M3,maxe,NB0,nk
      parameter (maxb=100,M3=100,nk0=200,maxe=200)
      integer*4 nv(*),mv(*),N7,i,kmax,N,Nexc,iex,j,km,
     1state(M3),si(maxe),kstate(M3),ic,LEXCL,mm,ii,im,kkr(maxb,2*nk0+1)
      real*8 w(M3,M3),u(*),fcsumr(*),xkr(maxb,2*nk0+1),
     1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp,
     1OH,t,xkrn,xkin,ak,dk,dko,
     1ykr(maxb,2*nk0+1),yki(maxb,2*nk0+1)
      real*8,allocatable::ptr(:),pti(:)
c     fcf Frank-Condon factor
c     mv,nv  excited and ground quantum numbers (1=ground state,nosc=0)
c     nk number of k-expansion points (k=-nk ... nk)
c     qs=wql+u
c     LEXCL maximum sum of k-numbers
c     xkr,xki real and imaginary parts of the exp. coefficients
c     N7..N number of oscillators, typically N7=1,N=3*nat-6
c     box=L the NM interval (-L/2...L/2)
c     m1 maximum excitation on one center
      allocate(ptr(NB0),pti(NB0))
      pi=4.0d0*atan(1.0d0)
      tp=2.0d0*pi

      do 1 i=1,mm
1     fcsumr(i)=0.0d0

c     distribute k-"excitations" on electronic excited coeff
c     start with ground as it may need less k-points
      do 3 Nexc=0,LEXCL
c     initial indices - all to the first center:
      do 41 iex=1,Nexc
41    si(iex)=N7
5     DO  1966 I=N7,N
1966  STATE(I)=0

      do 4 iex=1,Nexc
4     state(si(iex))=state(si(iex))+1

      do 42 iex=1,Nexc
42    if(state(si(iex)).gt.kmax)goto 333
c
c     calculate the product xl...ground state the same for all batch
      prs=1.0d0
      pis=0.0d0
      do 11 i=N7,N
c     0 .. largest, 1- second largest, ...
      kstate(I)=kkr(nv(i),state(I)+1)
      t  = prs*ykr(nv(i),1+state(i))-pis*yki(nv(i),1+state(i))
      pis= prs*yki(nv(i),1+state(i))+pis*ykr(nv(i),1+state(i))
      prs=t
11    if(abs(prs)+abs(pis).lt.OH)goto 333
c
c     calculate the product xs...
      do 7 ii=1,mm
      ptr(ii)=box**(N-N7+1)
7     pti(ii)=0.0d0

      fim=0.0d0
      do 111 i=N7,N
c     calculate corresponding ak:
      ak=0.0d0
      do 12 j=N7,N
c     w-first index excited
12    ak=ak+dble(kstate(j))*w(i,j)
c     contribution to the phase:
      fim=fim+u(i)*ak*tp/box
      if(ak.lt.0)then
c      nearest lower integer:
       km=int(ak)-1
      else
       km=int(ak)
      endif
      if(abs(km).gt.nk)goto 333
      if(abs(km+1).gt.nk)goto 333
      dk=ak-dble(km)
      dko=1.0d0-dk
      do 111 ii=1,mm
      im=mv(ii+(i-1)*NB0)
c     coefficients by interpolation:
      xkrn=xkr(im,nk+1+km)*dko+dk*xkr(im,nk+1+km+1)
      xkin=xki(im,nk+1+km)*dko+dk*xki(im,nk+1+km+1)
      t      = ptr(ii)*xkrn-pti(ii)*xkin
      pti(ii)= ptr(ii)*xkin+pti(ii)*xkrn
111   ptr(ii)= t

      do 6 ii=1,mm
      t      = ptr(ii)*cos(fim)-pti(ii)*sin(fim)
      pti(ii)= ptr(ii)*sin(fim)+pti(ii)*cos(fim)
      ptr(ii)= t
6     fcsumr(ii)=fcsumr(ii)+(prs*ptr(ii)+pis*pti(ii))

c     find index to be changed
333   if(Nexc.eq.0)goto 3
      do 8 ic=Nexc,1,-1
8     if(si(ic).lt.N)goto 9
      goto 3
9     do 102 i=ic+1,Nexc
102   si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
3     continue
c     loop end 22222222222222222222222222222222222222222222222222

      deallocate(ptr,pti)
      return
      end
c     ------------------------------------------------------------
      subroutine intq(n,k,r,i)
c     int(q^n exp(-q^2/2-ik.q)/sqrt(2*pi),q=-inf...inf)
c     r,i .. real and imag. parts
      implicit none
      integer*4 n,ni
      real*8 k,r,i,anir,anii,anipr,anipi,anippr,anippi
      ni=0
      anir=1.0d0
      anii=0.0d0
      if(n.eq.ni)then
       r=anir*exp(-k**2/2.0d0)
       i=anii
       return
      endif
      anipr=0.0d0
      anipi=-k
      if(n.eq.ni+1)then
       r=anipr
       i=anipi*exp(-k**2/2.0d0)
       return
      endif
c     I(n+2) = -i k I(n+1) + (n+1) I(n):
1     anippr= k*anipi+dble(ni+1)*anir
      anippi=-k*anipr+dble(ni+1)*anii
      if(n.eq.ni+2)then
       r=anippr*exp(-k**2/2.0d0)
       i=anippi*exp(-k**2/2.0d0)
       return
      else
       ni=ni+1
       anir=anipr
       anii=anipi
       anipr=anippr
       anipi=anippi
       goto 1
      endif
      end
c     ------------------------------------------------------------
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
c     ------------------------------------------------------------
      function fik(k,x)
      implicit none
      real*8 fik,x,w,hk,p,ank,ff
      integer maxb,M3
      parameter (maxb=100,M3=100)
      common/ho/p(maxb,maxb),ank(maxb),w(M3)
      integer*4 k,kk
      ff=ank(k)*exp(-x**2/2.0d0)
      hk=0.0d0
      do 1 kk=1,k
1     hk=hk+p(k,kk)*x**(kk-1)
      fik=ff*hk
      return
      end
c     ------------------------------------------------------------
      subroutine  ALLIGN(N7,N,W,UU,e1,e2,Wi,IWR,WT)
      implicit none
      INTEGER*4 NTA,MENDELEV,MXB,i,ia,ix,iy,iz,ii,NAT,
     1ierr,iter,j,MB,maxiter,itrial,
     2IX3,M3,NM,iq,idum,iix,N7,N,IWR
      PARAMETER (M3=100,NTA=M3/3,MENDELEV=89,MXB=25,maxiter=1000)
      real*8 pi,A,XS,XB,RTOL,am,WT(M3,M3),
     1oxx,amxx,ck,sk,dr1,dr2,fi,fk,pp,rat,wx,yy,z0,
     2tol,rr(NTA,3),c,s2r(M3,M3),jmat(M3,M3),kmat(M3),Wi(M3,M3),
     3X(NTA),r(NTA,3),rg(NTA,3),d(NTA,3),u(3),re(NTA,3)
      integer*4 ITY(NTA)
      real*4 amas
      dimension amas(MENDELEV)
      data amas/1.007825,4.003,
     2  6.941, 9.012, 10.810,12.0  ,14.003074,15.994915,18.9984,20.179,
     3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948,
     4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847,
     4 58.933,58.700,63.546,65.380,
     4                  69.720,72.590,74.922,78.960,79.904,83.800,
     5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070,
     5 102.906,106.400,107.868,112.410,
     5                 114.82,118.69,121.75,127.600,126.905,131.300,
     6 132.905,137.330,138.906,
     6                 140.120,140.908,144.240,145.000,150.400,
     6 151.960,157.250,158.925,162.500,164.930,167.260,168.934,
     6 173.040,174.970,
     6 178.490,180.948,183.850,186.207,190.200,192.220,195.090,
     6 196.967,207.590,204.370,207.200,208.980,210.000,210.001,
     6 222.02,
     7 223.000,226.025,227.028/     
      data z0/0.0d0/
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,MB,ITER
      real*8 s1(M3,M3),s2(M3,M3),e1(M3),e2(M3),UU(*),W(M3,M3),
     1tp,clight,const,hpl,nav
      logical lex
C     
      pi=4.0d0*atan(1.0d0)
      tol=1.0d-7
      tp=2.0d0*pi

      OPEN(4,FILE='fg.inp',STATUS='OLD')
      READ(4,*)NM,idum,NAT
      if(NAT.gt.NTA)call report('too many atoms')
      DO 7 I=1,NAT
      READ(4,*)ITY(I),(rg(I,ix),ix=1,3)
7     X(I)=dble(AMAS(ITY(I)))
      read(4,*)
      do 1 I=1,NAT
      do 1 iq=1,NM
1     read(4,*)idum,idum,(s1(3*(I-1)+ix,iq),ix=1,3)
      read(4,*)
      read(4,*)(e1(iq),iq=1,NM)
      close(4)
      CLOSE(4)
      call ordermodes(e1,s1,NAT,NM,M3)

      inquire(file='FTRY.INP',exist=lex)
      if(lex)then
       open(4,file='FTRY.INP')
       read(4,*)
       read(4,*)
       read(4,444)(X(I),I=1,NAT)
       read(4,*)
       read(4,444)(X(I),I=1,NAT)
444    format(6F12.6)
       close(4)
       write(6,*)'Masses read from FTRY.INP'
      endif
c
      AM=z0
      do 3 i=1,nat
3     AM=AM+X(i)
      WRITE(6,609)NAT,AM
609   format(' fg.inp read,',i4,' ATOMS, MOLECULAR MASS ',f10.3)

      OPEN(4,FILE='fe.inp',STATUS='OLD')
      READ(4,*)
      DO 72 I=1,NAT
72    READ(4,*)ITY(I),(re(I,ix),ix=1,3)
      read(4,*)
      do 2 I=1,NAT
      do 2 iq=1,NM
2     read(4,*)idum,idum,(s2(3*(I-1)+ix,iq),ix=1,3)
      read(4,*)
      read(4,*)(e2(iq),iq=1,NM)
      close(4)
      CLOSE(4)
      write(6,*)'fe.inp read'
      call ordermodes(e2,s2,NAT,NM,M3)


      do 303 ix=1,3
      c=z0
      do 300 ia=1,NAT
300   c = c+rg(ia,ix)*X(ia)
      do 303 ia=1,NAT
303   rg(ia,ix)=rg(ia,ix)-c/AM
      write(6,*)'fg put to origin'

      do 203 ix=1,3
      c=z0
      do 200 ia=1,NAT
200   c=c+re(ia,ix)*X(ia)
      do 203 ia=1,NAT
203   re(ia,ix)=re(ia,ix)-c/AM
      write(6,*)'fe put to origin'
c
c     rotate re to rg
      do 2031 ia=1,NAT
      do 2031 ix=1,3
2031  r(ia,ix)=re(ia,ix)

      itrial=1
88881 iter=0
8888  iter=iter+1
      if(iter.gt.maxiter)then
       write(6,*)'Maximum number of iterations exceeded'
       itrial=itrial+1
       if(itrial.lt.4)then
        write(6,*)'Trial ',itrial
        goto 88881
       endif
       write(6,*)'Iteration did not converge'
       goto 881
      endif

      do 210 ia=1,NAT
      do 210 ix=1,3
210   d(ia,ix)=r(ia,ix)-rg(ia,ix)
c
c     loop over the three axes,(try three different orderings):

      do 208 IX3=1,3
      IX=itrial+IX3-1
      if(IX.eq.4)IX=1
      if(IX.eq.5)IX=2
      IY=IX+1
      if(IY.eq.4)IY=1
      IZ=IY+1
      if(IZ.eq.4)IZ=1
      oxx=z0
      amxx=z0
      do 209 ia=1,NAT
      dr1=dsqrt(rg(ia,IY)**2+rg(ia,IZ)**2)
      dr2=dsqrt( r(ia,IY)**2+ r(ia,IZ)**2)
      pp=rg(ia,IY)*r(ia,IY)+rg(ia,IZ)*r(ia,IZ)
      if(dabs(dr1*dr2).gt.1.0d-8)then
       rat=pp/dr1/dr2
       if(rat.gt. 1.0d0-1.0d-13)then
        fi=z0
       else
        if(rat.lt.-1.0d0+1.0d-13)then
         fi=pi
        else
         fi=acos(rat)
        endif
       endif
       wx=rg(ia,IY)*d(ia,IZ)-rg(ia,IZ)*d(ia,IY)
       if(wx.lt.0.0d0)fi=-fi
c      amxx=amxx+X(ia)*fi*(5.0d0*dr1**2+2.0d0*dr2**2-dr1*dr2)/6.0d0
       amxx=amxx+X(ia)*fi*(dr1**2+dr2**2+dr1*dr2)/3.0d0
       oxx=oxx+X(ia)*dr2**2
      endif
209   continue

      if(dabs(oxx).gt.1.0d-8)then
       fk=amxx/oxx
      else
       fk=z0
      endif
      u(IX)=fk
c
c     rotate back by this angle:
      ck=cos(-fk)
      sk=sin(-fk)
      do 208 ia=1,NAT
      YY=r(ia,IY)
      r(ia,IY)=YY*ck-r(ia,IZ)*sk
208   r(ia,IZ)=YY*sk+r(ia,IZ)*ck

c     write(6,*)iter,u(1)**2+u(2)**2+u(3)**2
      if(u(1)**2+u(2)**2+u(3)**2.gt.tol)goto 8888

881   write(6,608)iter
608   format(i7,' iterations')

      open(41,file='tr.x')
      write(41,*)'rotated coordinates of exc'
      write(41,*)NAT
      do 2043 ia=1,NAT
2043  write(41,4141)ITY(ia),(r(ia,ix),ix=1,3)
      close(41)

      write(6,*)
      MB=NAT
      do 217 ii=1,MB
      do 217 ix=1,3
      XB(ii,ix)=r(ii,ix)
217   XS(ii,ix)=re(ii,ix)
      call DOMATRIX(IERR)
c     XB = A XS
c     rotated Xe = A original Xe
      write(6,6004)A
6004  format(' Transformation matrix:',/,3(3F10.4,/))

c     mass-weighed S-matrices
      do 214 ia=1,NAT
      do 214 ix=1,3
      ii=3*(ia-1)+ix
      do 214 iq=1,NM
      s1(ii,iq)=s1(ii,iq)*dsqrt(X(ia))
214   s2(ii,iq)=s2(ii,iq)*dsqrt(X(ia))
c     now transpose is inverted
      
      do 213 ia=1,NAT
      do 213 ix=1,3
      ii=3*(ia-1)+ix
      rr(ia,ix)=0.0d0
      do 2131 iix=1,3
2131  rr(ia,ix)=rr(ia,ix)+A(ix,iix)*r(ia,iix)
      do 213 iq=1,NM
      s2r(ii,iq)=0.0d0
      do 213 iix=1,3
213   s2r(ii,iq)=s2r(ii,iq)+A(ix,iix)*s2(3*(ia-1)+iix,iq)
      write(6,*)'fe transformed'

      open(41,file='xr.x')
      write(41,*)'rotated coordinates of exc'
      write(41,*)NAT
      do 2042 ia=1,NAT
2042  write(41,4141)ITY(ia),(rr(ia,ix),ix=1,3)
4141  format(i4,3f12.6)
      close(41)

      open(41,file='xg.x')
      write(41,*)'ground geometry in origin'
      write(41,*)NAT
      do 204 ia=1,NAT
204   write(41,4141)ITY(ia),(rg(ia,ix),ix=1,3)
      close(41)

      call fixphase(s1,s2r,NAT,NM,M3)

      do 215 i=1,NM
      do 215 j=1,NM
      w(i,j)=0.0d0
      wi(i,j)=0.0d0
      jmat(i,j)=0.0d0
      do 2151 ii=1,3*nat
c     w,wi - first index always excited
c     transformation matrix between dimensionless normal mode coords
      jmat(i,j)=jmat(i,j)+s2r(ii,i)*s1(ii,j)
      if(e1(j).gt.0.0d0.and.e2(j).gt.0.0d0)then
       wi(i,j)=wi(i,j)+s2r(ii,i)*dsqrt(e1(j)/e2(i))*s1(ii,j)
       w(i,j)=w(i,j)+s2r(ii,i)*dsqrt(e2(i)/e1(j))*s1(ii,j)
      endif
2151  continue
215   continue
      clight=3.0d10
      nav=6.022d23
      hpl=1.054d-34
      const=dsqrt(1.0d-3*tp*clight/nav/hpl)*1.0d-10
      do 218 i=1,NM
      uu(i)=0.0d0
      kmat(i)=0.0d0
      do 218 ia=1,nat
      do 218 ix=1,3
      ii=3*(ia-1)+ix
      kmat(i)=kmat(i)-s2r(ii,i)*(rg(ia,ix)-rr(ia,ix))*dsqrt(X(ia))
218   if(e2(i).gt.0)uu(i)=uu(i)+
     1dsqrt(e2(i)*X(ia))*s2r(ii,i)*const*(rg(ia,ix)-rr(ia,ix))
      if(IWR.gt.0)then
       write(6,6007)
6007   format(/,' J-matrix:')
       do 221 j=1,3*NAT-6
221    write(6,6008)(jmat(i,j),i=1,3*NAT-6)
       write(6,6009)
6009   format(/,' K-vector:')
       write(6,6008)(kmat(i),i=1,3*NAT-6)
6008   format(10f8.2)
       write(6,6005)
6005   format(' mode      Eg      Ee  shift',/,'----------------')
       do 219 i=1,NM
219    write(6,6006)i,e1(i),e2(i),uu(i)
6006   format(i5,2f8.1,f9.3)
      endif

      if(N7.eq.0)then
       N7=7
       N=NM
      endif

      do 8 i=N7,NM
      do 8 j=N7,NM
8     wt(i-N7+1,j-N7+1)=wi(i,j)

      return
      
      END

      SUBROUTINE DOMATRIX(IERR)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=25)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.0d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(ANG)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(ANG)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(ANG)
      ANG(1)=ANG(3)
      ANG(2)=ANG(3)
      ANG(3)=0.0d0
      DO 4 I=1,3
4     P(4,I)=ANG(I)
      Y(4)=FU(ANG)
      FTOL=1.0d-7
      ITMAX=1500
      ITER=0
99999 ILO=1
      IF(Y(1).GT.Y(2))THEN
       IHI=1
       INHI=2
      ELSE
       IHI=2
       INHI=1
      ENDIF
      DO 5 I=1,4
      IF(Y(I).LT.Y(ILO))ILO=I
      IF(Y(I).GT.Y(IHI))THEN
       INHI=IHI
       IHI=I
      ELSE
       IF(Y(I).GT.Y(INHI))THEN
        IF(I.NE.IHI)INHI=I
       ENDIF
      ENDIF
5     CONTINUE
      IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-11)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX)THEN
       WRITE(3,*)' Rotation has not converged !'
       IERR=1
       RETURN
      ENDIF
      ITER=ITER+1
      DO 55 I=1,3
55    PBAR(I)=0.0d0
      DO 6 I=1,4
      IF(I.NE.IHI)THEN
       DO 7 J=1,3
7      PBAR(J)=PBAR(J)+P(I,J)
      ENDIF
6     CONTINUE
      DO 8 J=1,3
      PBAR(J)=PBAR(J)/3.0d0
8     PR(J)=2.0d0*PBAR(J)-P(IHI,J)
      YPR=FU(PR)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(PRR)
       IF(YPRR.LT.Y(ILO))THEN
        DO 10 J=1,3
10      P(IHI,J)=PRR(J)
        Y(IHI)=YPRR
       ELSE
        DO 11 J=1,3
11      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ELSE
       IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
         DO 12 J=1,3
12       P(IHI,J)=PR(J)
         Y(IHI)=YPR
        ENDIF
        DO 13 J=1,3
13      PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J)
        YPRR=FU(PRR)
        IF(YPRR.LT.Y(IHI))THEN
         DO 14 J=1,3
14       P(IHI,J)=PRR(J)
         Y(IHI)=YPRR
        ELSE
         DO 15 I=1,4
         IF(I.NE.ILO)THEN
          DO 16 J=1,3
          PR(J)=0.5d0*(P(I,J)+P(ILO,J))
16        P(I,J)=PR(J)
          Y(I)=FU(PR)
         ENDIF
15       CONTINUE
        ENDIF
       ELSE
        DO 17 J=1,3
17      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ENDIF
      GOTO 99999
      END

      FUNCTION FU(ANG)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=25)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      DIMENSION ANG(3)
      S1=SIN(ANG(1))
      S2=SIN(ANG(2))
      S3=SIN(ANG(3))
      C1=COS(ANG(1))
      C2=COS(ANG(2))
      C3=COS(ANG(3))
C     ang1 = theta
C     ang2 = phi
C     ang3 = cappa
      A(1,1)=C1*C2*C3-S2*S3
      A(1,2)=C1*S2*C3+C2*S3
      A(1,3)=-S1*C3
      A(2,1)=-C1*C2*S3-S2*C3
      A(2,2)=-C1*S2*S3+C2*C3
      A(2,3)=S1*S3
      A(3,1)=S1*C2
      A(3,2)=S1*S2
      A(3,3)=C1
      R=0.0d0
      DO 1 I=1,NAT
      TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3)
      TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3)
      TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3)
1     R=R+(TX-XB(I,1))**2+(TY-XB(I,2))**2+(TZ-XB(I,3))**2
      FU=R
      RETURN
      END
c
      SUBROUTINE DTRM(A,M,N,D,INDX,C)
C
C Subroutine for evaluating the determinant of a matrix using 
C the partial-pivoting Gaussian elimination scheme.
C
      implicit none
      integer*4 M,N,INDX(M),I,MSGN,J
      real*8 A(M,M),D,C(M)
C
      CALL ELGS(A,M,N,INDX,C)
      D    = 1.0d0
      DO     100 I = 1, N
100   D = D*A(INDX(I),I)
      MSGN = 1
      DO     200 I = 1, N
      DO   150 WHILE (I.NE.INDX(I))
      MSGN = -MSGN
      J = INDX(I)
      INDX(I) = INDX(J)
      INDX(J) = J
  150 END DO
  200 CONTINUE
      D = dble(MSGN)*D
      RETURN
      END
C
      SUBROUTINE ELGS(A,M,N,INDX,C)
C
C Subroutine to perform the partial-pivoting Gaussian elimination.
C A(N,N) is the original matrix in the input and transformed
C matrix plus the pivoting element ratios below the diagonal in
C the output.  INDX(N) records the pivoting order.

C
      implicit none
      integer*4 M,N,INDX(M),I,J,K,ITMP
      real*8 A(M,M),C(M),C1,PI1,PI,PJ
C
      DO     50    I = 1, N
50    INDX(I) = I
C
C Find the rescaling factors, one from each row
C
      DO     100   I = 1, N
      C1= 0.0d0
      DO    90   J = 1, N
90    C1 = MAX(C1,DABS(A(I,J)))
100   C(I) = C1
C
C Search the pivoting (largest) element from each column
C
      DO     200   J = 1, N-1
        PI1 = 0.0d0
        DO   150   I = J, N
          PI = ABS(A(INDX(I),J))/C(INDX(I))
          IF (PI.GT.PI1) THEN
            PI1 = PI
            K   = I
          ELSE
          ENDIF
  150   CONTINUE
C
C Interchange the rows via INDX(N) to record pivoting order
C
        ITMP    = INDX(J)
        INDX(J) = INDX(K)
        INDX(K) = ITMP
        DO   170   I = J+1, N
          PJ  = A(INDX(I),J)/A(INDX(J),J)
C
C Record pivoting ratios below the diagonal
C
          A(INDX(I),J) = PJ
C
C Modify other elements accordingly
C
          DO 160   K = J+1, N
            A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
  160     CONTINUE
  170   CONTINUE
  200 CONTINUE
C
      RETURN
      END
c     ------------------------------------------------------------
      subroutine ordermodes(e,s,NAT,NM,M3)
      integer*4 NAT,NM,i,j,k,M3
      real*8 s(M3,M3),e(*),t
      do 1 i=1,NM
      do 1 j=i+1,NM
      if(e(j).lt.e(i))then
       t=e(j)
       e(j)=e(i)
       e(i)=t
       do 2 k=1,3*NAT
       t=s(k,j)
       s(k,j)=s(k,i)
2      s(k,i)=t
      endif
1     continue
      return
      end
c     ------------------------------------------------------------
      subroutine fixphase(u,s,NAT,NM,M3)
      integer*4 NAT,NM,i,k,M3
      real*8 s(M3,M3),u(M3,M3),t
      do 1 i=1,NM
      t=0.0d0
      do 2 k=1,3*NAT
2     t=t+s(k,i)*u(k,i)
      if(t.lt.0)then
       do 3 k=1,3*NAT
3      s(k,i)=-s(k,i)
      endif
1     continue
      return
      end
c     function mclock()
c     mclock=0
c     return
c     end
c     ------------------------------------------------------------
      function frandom()
      implicit none
      real*8 frandom
c
c     Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2
c
c     Returns a pseudo-random number rectangularly distributed
c     between 0 and 1.   The cycle length is 6.95E+12 (See page 123
c     of Applied Statistics (1984) vol.33), not as claimed in the
c     original article.
c
c     IX, IY and IZ should be set to integer values between 1 and
c     30000 before the first entry.
c
      integer*4 ix, iy, iz
      common /randc/ ix, iy, iz
c
      ix = mod(171 * ix, 30269)
      iy = mod(172 * iy, 30307)
      iz = mod(170 * iz, 30323)
      frandom = mod(dble(ix) / 30269.0d0 + dble(iy) / 30307.0d0 +
     1                        dble(iz) / 30323.0d0, 1.0d0)
      return
      end
c     ------------------------------------------------------------
      SUBROUTINE LOADC(s,J,CFJ,N)
      INTEGER*4 J,I,N 
      REAL*8 CFJ(*)
      character*(*) s
      OPEN(45,FILE=s,FORM='UNFORMATTED',STATUS='OLD')
      DO 1 K=1,J-1
1     READ(45)
      READ(45,err=999)(CFJ(I),I=1,N)
      CLOSE(45)
      RETURN 
999   write(6,*)J,N,I
      call report('reading error')
      END