PROGRAM stsim
      implicit none
      integer*4 N0,A,B,K,NOAT,NQ,NAC,I,J,L,IAT,IMD,IMDMAX,
     1LL,NM0,n3,iu,i1,i2,IX,IY,NX,NY,iargc
      PARAMETER (N0=15,NM0=3*N0)
      REAL*8 SV(N0,NM0,3),PP(N0,3,3),II(N0,3,3),JJ(N0,3,3),
     1 PV(N0,3,3),I0(N0,3,3),R(3,N0),FRQ(NM0),ELD(3),MGD(3),MDN(3),
     2MDD(3),MDT(3),BOHR,FD,FC,WLIM,W,WR,ABSW,x,y,RDO(NM0),
     3DS(NM0),ymax,ymin,p,Q,PEF,T,pds,rds,
     8d(NM0**3),c(NM0**3),dx,dy,xmax,xmin,Emm,E
      character*80 sa
c     
      FD =SQRT(388.9219d0)
      FC=131.14D-8
      BOHR=0.52917715d0
      WLIM=0.000001d0
      WRITE(*,*)' Spectra simulation - anharmonic and temperature'
      WRITE(*,*)'         corrections - experimental code'
      WRITE(*,6666)
6666  FORMAT(/,
     1'  INPUT FILES: F   .INP ... VIBRATIONAL PARAMETERS',/,
     1'               FILE.TEN ... TENSORS',/,
     1'               STSIM.OPT ... options',/,
     3'  OUTPUT FILE: DOG.TAB',/,
     4'               DIP.MOM',/)

      open(9,file='STSIM.OPT')
      read(9,*)i1,i2,nx,ny,xmin,xmax,ymin,ymax
c     i1,i2 ... modes of interest
c     nx,ny .. number of points in each direction
c     dx,dy .. increments in dimensionless coordinates
c     xmin,xmax,ymin,ymax ... integration limits
      dx=(xmax-xmin)/dble(nx)
      dy=(ymax-ymin)/dble(ny)
      close(9)
      if(iargc().gt.0)then
       call getarg(1,sa)
       read(sa,*)T
      else
       T=298.0D0
      endif

      write(6,444)T
 444  format('Temperature ',f6.1,' K')
       
      OPEN(15,FILE='FILE.TEN',STATUS='OLD')
      READ (15,*)NOAT
      OPEN(17,FILE='F.INP',STATUS='OLD')
      READ(17,*)NQ,NAC
      IF(NAC.NE.3*NOAT)call 
     1report(' Number of atoms if FILE.TEN and F.INP is different !')
      IF(NOAT.gt.N0)call report('Too many atoms')
      DO 35 K=1,NOAT
      READ(17,*)R(1,K),(R(I,K),I=1,3)
      DO 35 I=1,3
35    R(I,K)=R(I,K)/BOHR
      WRITE(*,*)' Geometry read in'
      IMDMAX=0
      READ(17,*)
      if(NQ.gt.NM0)call report('too many modes')
      DO 45 LL=1,NOAT
      DO 45 I=1,NQ
      READ (17,*) IAT,IMD,SV(IAT,I,1),SV(IAT,I,2),SV(IAT,I,3)
      IF (IMD.GT.IMDMAX)IMDMAX=IMD
45    CONTINUE
      READ(17,*)
      READ(17,*)(FRQ(I),I=1,NQ)
      WRITE(*,*)' S-matrix read in'
C
C  READ DIPOLE DERIVATIVES = TENSOR P(A,B)
C
      DO 10 L=1,NOAT
      DO 10 J=1,3
10    READ (15,*) (PP(L,J,I),I=1,3)
C
C  READ TENSORS I(A,B) AND J(A,B) FOR THE 0 0 0        ORIGIN:
C    (=ELECTRONIC & NUCLEAR TERMS)
C
      DO 220 L=1,NOAT
      DO 220 J=1,3
 220  READ (15,*) (I0(L,J,I),I=1,3)
      DO 230 L=1,NOAT
      DO 230 J=1,3
230   READ (15,*) (JJ(L,J,I),I=1,3)
C
C  READ VELOCITY REPRESENTATION OF DIPOLE DERIVATIVES
C
      DO 100 L=1,NOAT
      DO 100 J=1,3
100   READ (15,*) (PV(L,J,I),I=1,3)
      CLOSE(15)
      WRITE(*,*)' Dipole derivatives read in'
      DO 140 L=1,NOAT
      DO 140 A=1,3
      DO 140 B=1,3
140   II (L,A,B) = I0(L,A,B)
      CALL VCDD0(II,PV,PP,R,NOAT,N0)
      WRITE(*,*)' Tensors transformed'

      call inpcd('CQQ.SCR.TXT',c,NOAT,FRQ)
      call inpcd('DQQ.SCR.TXT',d,NOAT,FRQ)

C
      DO 413 I=1,NQ
      DS(I)=0.0d0
      RDO(I)=0.0d0
413   continue

      DO 41 I=1,NQ
      W=FRQ(I)
      ABSW=DABS(W)
      WR=DSQRT(ABSW)
      DO 56 B=1,3
      MGD(B)=0.0d0
      ELD(B)=0.0d0
      MDN(B)=0.0d0
      MDD(B)=0.0d0
      MDT(B)=0.0d0
      DO 56 L=1,NOAT
      DO 56 A=1,3
C     ELECTRONIC MOMENT TOTAL:
      ELD(B) = ELD(B) +             PP(L,A,B)*SV(L,I,A)
C     MAGNETIC DIPOLE NUCLEAR    SAME FOR DO/GO
      MDN(B) = MDN(B) +            JJ(L,A,B)  *SV(L,I,A)
C     MAGNETIC DIPOLE ELECTRONIC DO
      MDD(B) = MDD(B) +             II(L,A,B)*SV(L,I,A)
C     MAGNETIC DIPOLE TOTAL      DO
      MDT(B) = MDT(B) + (JJ(L,A,B)+ II(L,A,B))*SV(L,I,A)
C     MAGNETIC DIPOLE TOTAL  GO
56    MGD(B) = MGD(B) + (I0(L,A,B)+JJ(L,A,B)) *SV(L,I,A)
      DO 51 B=1,3
      IF(ABSW.GT.WLIM)ELD(B)=ELD(B)/WR*FD
      MDT(B)=MDT(B)*FC*WR
      MDN(B)=MDN(B)*FC*WR
      MDD(B)=MDD(B)*FC*WR
51    MGD(B)=MGD(B)*FC*WR
      DS(I) = ELD(1)*ELD(1)+ELD(2)*ELD(2)+ELD(3)*ELD(3)
      RDO(I)= MDT(1)*ELD(1)+MDT(2)*ELD(2)+MDT(3)*ELD(3)
41    IF(ABSW.LT.WLIM)DS(I)=0.0d0
c
      n3=3*NOAT

c     precalculate probability sum
      Q=0.0d0
      E=0.0d0
      x=xmin-dx
      do 414 ix=1,nx
      x=x+dx
      y=ymin-dy
      do 414 iy=1,ny
      y=y+dy
414   Q=Q+PEF(x,y,FRQ,NQ,c,d,n3,E,T,i1,i2)
c
c     write probability on disk
      open(93,file='p.prn')
      x=xmin-dx
      do 415 ix=1,nx
      x=x+dx
      y=ymin-dy
      do 415 iy=1,ny
      y=y+dy
      p=PEF(x,y,FRQ,NQ,c,d,n3,E,T,i1,i2)/Q
415   write(93,9393)x,y,E,p
9393  format(2F15.5,f10.1,f15.8)
      close(93)

      OPEN(18,FILE='DOG.TAB')
      WRITE (18,501)
  501 FORMAT (/,'          FREQ.     ',
     1'  D [D^2]     ','  R [D^-2]DO         R - CO      I[km/mol]',
     2/,'-----------------------------------------------------------')
      DO 412 I=1,NQ
c     spread each frequency over n1,n2 boltzmann surface
c     E=sum(i)wi xi^2/2 + sum(i,j,k)cijk xi xj xk/6 + sum (i,j,k,l)dijkl
c     xi xj xk xl /24
c      =sum(i)wi xi^2/2 + sum(i,j,k,all different)cijk xi xj xk/6 
c     +sum(i<>j)ciij xi^2 xj/2+ sum(i)ciii xi^3/6
c     +sum (i,j,k,l,all different (not considered))dijkl/24 xi xj xk xl
c     +sum (i,j,k,all different)diijk xi^2 xj xk/2
c     +sum (i<>j)diiij xi^3 xj/6 + sum (i<>j)diijj xi^2 xj^2/4
c     +sum (i)diiii xi^4 /24
c     Em=wm xm + sum(i<>m,j<>m)cmjk xi xj/2 
c     + sum(i<>m)cmii xi^2/2+ sum(i<>m)cmmi xm xi+ cmmm xm^2/2+
c     +sum (i,j,k,all different (not considered))dijkm xi xj xk /6
c     +sum (i<>j)dmmij xm xi xj +sum (i<>j)diimj xi^2 xj
c     +sum (i<>m)dmmmi xm^2 xi/2 +sum (i<>m)diiim xi^3 /6 
c     +sum (i<>m)diimm xm xi^2 +dmmmm xm^3/6
c     n<>m:
c     Emn= cmnn xn+ cmmn xm
c     +sum (i,j,all different (not considered))dijmn xi xj /2
c     +2 sum (i<>n)dmmin xm xi +sum (i<>n)diimn xi^2 +2sum (i<>n)dnnmi xn xi
c     +dmmmn xm^2 /2 +dnnnm xn^2 /2  +2 dnnmm xm xn 
c     mm:
c     Emm= wm+sum(j<>m)cmmj xj+ cmmm xm
c     +sum (i<>j)dmmij xi xj+sum (i<>m)(dmmmi xm xi +diimm xi^2)+dmmmm xm^2/2
c     x... mode 7, y ... mode 8
c     numbering according to AH constants:
      iu=NQ-i+1+6
      x=xmin-dx
      do 412 ix=1,nx
      x=x+dx
      y=ymin-dy
      do 412 iy=1,ny
      y=y+dy
      p=PEF(x,y,FRQ,NQ,c,d,n3,E,T,i1,i2)/Q
      Emm=FRQ(I)
      if(iu.ne.i1)Emm=Emm+c(i1+n3*(iu-1)+n3**2*(iu-1))*x
      if(iu.ne.i2)Emm=Emm+c(i2+n3*(iu-1)+n3**2*(iu-1))*y
      if(iu.eq.i2)Emm=Emm+c(i2+n3*(i2-1)+n3**2*(i2-1))*y
      if(iu.eq.i1)Emm=Emm+c(i1+n3*(i1-1)+n3**2*(i1-1))*x
      if(iu.ne.i1.and.iu.ne.i2)Emm=Emm+d(iu+n3*(i1-1)+n3**2*(i2-1))*y*x
      if(iu.eq.i1)Emm=Emm+d(iu+n3*(iu-1)+n3**2*(i2-1))*x*y
      if(iu.eq.i1)Emm=Emm+d(iu+n3*(i2-1)+n3**2*(i2-1))*y*y
      if(iu.eq.i2)Emm=Emm+d(iu+n3*(iu-1)+n3**2*(i1-1))*x*y
      if(iu.eq.i2)Emm=Emm+d(iu+n3*(i1-1)+n3**2*(i1-1))*x*x
      if(iu.eq.i1)Emm=Emm+d(i1+n3*(i1-1)+n3**2*(i1-1))*x*x/2.0d0
      if(iu.eq.i2)Emm=Emm+d(i2+n3*(i2-1)+n3**2*(i2-1))*x*x/2.0d0
      pds=p*DS(I)
      rds=p*RDO(I)
      if(pds.lt.1.0d-20)pds=0.0d0
      if(rds.lt.1.0d-20)rds=0.0d0
412   if(iu.eq.i2)WRITE (18,523) I,Emm,pds,rds
523   FORMAT(I5,f16.7,4G15.6)
      WRITE(18,1819)
1819  FORMAT('-------------------------------------------------')
      CLOSE(18)
C
      CLOSE(17)
      WRITE(*,*)' >> PROGRAM TERMINATED <<'
      STOP
      END

      SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT,N00)
C     STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI
      IMPLICIT none
      INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA,N00,NAT
      real*8 VCD(N00,3,3),VEL(N00,3,3),DIP(N00,3,3),
     1C(3,N00),SKEW,EPS
      DO 90 BETA =1,3
      DO 90 GAMMA=1,3
      DO 90 DELTA=1,3
      SKEW=0.25D0*EPS(BETA,GAMMA,DELTA)
      IF (DABS(SKEW).GT.1.0D-10)THEN
      DO 80 ALPHA=1,3
      DO 80 LAMBDA=1,NAT
80    VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA)
     1 -SKEW*C(GAMMA,LAMBDA)
     2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA))
      ENDIF
90    CONTINUE
      RETURN
      END

      FUNCTION EPS(I,J,K)
      integer*4 I,J,K
      REAL*8 EPS
      EPS=0.0D0
      IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0D0
      IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0D0
      IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0D0
      IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0D0
      IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0D0
      IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0D0
      RETURN
      END


      SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (N0=15000)
      DIMENSION A(N,N),Z(N,N),D(N),E(3*N0)
      COMMON/DIAG/E
      IF (IEIG .LT. 0) GO TO 110
      DO 100 I = 1, N
      DO 100 J = 1, I
      Z(I,J) = A(I,J)
  100 CONTINUE
  110 ICODE = IABS(IEIG)
      IF (N .EQ. 1) GO TO 320
      DO 300 II = 2, N
      I = N + 2 - II
      L = I - 1
      H = 0.0D0
      SCALE = 0.0D0
      IF (L .LT. 2) GO TO 130
      DO 120 K = 1, L
  120 SCALE = SCALE + DABS(Z(I,K))
      IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140
  130 E(I) = Z(I,L)
      GO TO 290
  140 DO 150 K = 1, L
      F = Z(I,K) / SCALE
      Z(I,K) = F
      H = H + F * F
  150 CONTINUE
      G = -DSIGN(DSQRT(H),F)
      E(I) = SCALE * G
      H = H - F * G
      Z(I,L) = F - G
      F = 0.0D0
      DO 240 J = 1, L
      IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H)
      G = 0.0D0
      DO 180 K = 1, J
  180 G = G + Z(J,K) * Z(I,K)
      JP1 = J + 1
      IF (L .LT. JP1) GO TO 220
      DO 200 K = JP1, L
  200 G = G + Z(K,J) * Z(I,K)
  220 E(J) = G / H
      F = F + E(J) * Z(I,J)
  240 CONTINUE
      HH = F / (H + H)
      DO 260 J = 1, L
      F = Z(I,J)
      G = E(J) - HH * F
      E(J) = G
      DO 260 K = 1, J
      Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260 CONTINUE
      DO 280 K = 1, L
  280 Z(I,K) = SCALE * Z(I,K)
  290 D(I) = H
  300 CONTINUE
  320 D(1) = 0.0D0
      E(1) = 0.0D0
      IF (ICODE .NE. 2) GO TO 600
      DO 500 I = 1, N
      L = I - 1
      IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380
      DO 360 J = 1, L
      G = 0.0D0
      DO 340 K = 1, L
  340 G = G + Z(I,K) * Z(K,J)
      DO 360 K = 1, L
      Z(K,J) = Z(K,J) - G * Z(K,I)
  360 CONTINUE
  380 D(I) = Z(I,I)
      Z(I,I) = 1.0D0
      IF (L .LT. 1) GO TO 500
      DO 400 J = 1, L
      Z(I,J) = 0.0D0
      Z(J,I) = 0.0D0
  400 CONTINUE
  500 CONTINUE
  620 CALL TQL12 (N,Z,D,IERR,ICODE)
      RETURN
  600 DO 610 I=1,N
  610 D(I) = Z(I,I)
      GO TO 620
      END
      SUBROUTINE TQL12(N,Z,D,IERR,ICODE)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (N0=15000)
      DIMENSION Z(N,N),D(N),E(3*N0)
      REAL*8 MACHEP
      COMMON/DIAG/E
      EPS = 1.0D0
   10 EPS = 0.50D0*EPS
      TOL1 = EPS + 1.0D0
      IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10
      IF(TOL1-EPS.EQ.1.0D0)  EPS = EPS + EPS
      MACHEP = EPS
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
      DO 100 I = 2, N
  100 E(I-1) = E(I)
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
      DO 240 L = 1, N
      J = 0
      H = MACHEP * (DABS(D(L)) + DABS(E(L)))
      IF (B .LT. H) B = H
      DO 110 M = L, N
      IF (DABS(E(M)) .LE. B) GO TO 120
  110 CONTINUE
  120 IF (M .EQ. L) GO TO 220
  130 IF (J .EQ. 30) GO TO 1000
      J = J + 1
      L1 = L + 1
      G = D(L)
      P = (D(L1) - G) / (2.0D0 * E(L))
      R = DSQRT(P*P+1.0D0)
      D(L) = E(L) / (P + DSIGN(R,P))
      H = G - D(L)
      DO 140 I = L1, N
  140 D(I) = D(I) - H
      F = F + H
      P = D(M)
      C = 1.0D0
      S = 0.0D0
      MML = M - L
      DO 200 II = 1, MML
      I = M - II
      G = C * E(I)
      H = C * P
      IF (DABS(P) .LT. DABS(E(I))) GO TO 150
      C = E(I) / P
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * P * R
      S = C / R
      C = 1.0D0 / R
      GO TO 160
  150 C = P / E(I)
      R = DSQRT(C*C+1.0D0)
      E(I+1) = S * E(I) * R
      S = 1.0D0 / R
      C = C * S
  160 P = C * D(I) - S * G
      D(I+1) = H + S * (C * G + S * D(I))
      IF (ICODE .NE. 2) GO TO 200
      DO 180 K = 1, N
      H = Z(K,I+1)
      Z(K,I+1) = S * Z(K,I) + C * H
      Z(K,I) = C * Z(K,I) - S * H
  180 CONTINUE
  200 CONTINUE
      E(L) = S * P
      D(L) = C * P
      IF (DABS(E(L)) .GT. B) GO TO 130
  220 D(L) = D(L) + F
  240 CONTINUE
      GO TO 1001
 1000 IERR = L
 1001 RETURN
      END

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end
      
c     ========================================
      subroutine inpcd(fn,c,nat,w)
      implicit none
      real*8 c(*),a,w(*),cm
      integer*4 nat,i,j,k,n,n3
      character*(*) fn
      n3=nat*3
      cm=219470.0d0
      do 2 i=1,n3**3
2     c(i)=0.0d0
      open(88,file=fn)
      n=0
1     read(88,*,end=999,err=999)i,j,k,a
      if(fn.eq.'CQQ.SCR.TXT')then
       a=a/dsqrt(abs(w(n3-i+1)*w(n3-j+1)*w(n3-k+1)/cm**3))
       write(6,*)'c',i,j,k,a
      else
       a=a/dsqrt(abs(w(n3-i+1)**2*w(n3-j+1)*w(n3-k+1)/cm**4))
       write(6,*)'d',i,i,j,k,a
      endif
      n=n+1
      c(i+(j-1)*n3+(k-1)*n3**2)=a
      c(i+(k-1)*n3+(j-1)*n3**2)=a
      if(fn.eq.'CQQ.SCR.TXT')then
       c(j+(k-1)*n3+(i-1)*n3**2)=a
       c(j+(i-1)*n3+(k-1)*n3**2)=a
       c(k+(i-1)*n3+(j-1)*n3**2)=a
       c(k+(j-1)*n3+(i-1)*n3**2)=a
      endif
      goto 1
  999 close(88)
      write(6,900)n,fn
900   format(i10,' constants read from ',A11)
      return
      end
c     ========================================
      function PEF(x,y,FRQ,NQ,c,d,n3,E,T,i1,i2)
      implicit none
      real*8 x,y,FRQ(*),c(*),d(*),PEF,E,T
      integer*4 n3,NQ,i1,i2
      E=FRQ(NQ)*x**2+FRQ(NQ-1)*y**2+
     1c(i1+n3*(i1-1)+n3**2*(i2-1))*y*x**2/2.0d0+
     1c(i2+n3*(i2-1)+n3**2*(i1-1))*y**2*x/2.0d0+
     1c(i1+n3*(i1-1)+n3**2*(i1-1))*x**3/6.0d0+
     1c(i2+n3*(i2-1)+n3**2*(i2-1))*y**3/6.0d0+
     1d(i1+n3*(i1-1)+n3**2*(i2-1))*x**3*y/6.0d0+
     1d(i2+n3*(i2-1)+n3**2*(i1-1))*y**3*x/6.0d0+
     1d(i2+n3*(i1-1)+n3**2*(i1-1))*y**2*x**2/2.0d0+
     1d(i1+n3*(i1-1)+n3**2*(i1-1))*x**4/24.0d0+
     1d(i2+n3*(i2-1)+n3**2*(i2-1))*y**4/24.0d0
c              h        *E*c     /R        /T*Nav
      PEF=exp(-6.626d-34*E*3.0d10/8.31441d0/T*6.022d23)
      return
      end