PROGRAM NEW4FTR5
c     with Raman extension
c     with dynamic allocation (F90)
      IMPLICIT none
      integer*4 nsim,it0,mclock
c
c     local variables:
      integer*4 i9,i7,io,iti,itj,itkl,l,ik,NA,NM,NOAT,NOSC,I,J,IA,iel,
     1nproc,omp_get_num_procs,omp_get_max_threads,OMP_GET_THREAD_NUM,
     1K,it,ii,NP,nmd,nhalf,ix,is,nff,inn,jn,nd,nrest,ipr,it00,seed,i8

      real*8  ur(51),ui(51),tr(51),ti(51),rtr,s1,s2,s3,fi,uix,urx,
     1ei,wi,t,ANORM,sf,CONST,c1,CHATOT,BM,frandom,ce,se,ann,
     5WMIN,WMAX,DELW,dt,dd,d,con,wj,FC,FD,c2,pi,SAL0,SAL1,SAG0,SAG1,
     1SA1,stemp,s4,di,ri
      logical lex,LG,LDW,LEB,LVEL,LSEC,LT3,lrtr
      character*5 s5

      integer*4,allocatable::jbb(:),j33(:),nstart(:),nend(:),
     1mm(:),nn(:)
      real*8,allocatable::abb(:),a33(:),sff(:),sar(:),st(:),
     1sti(:),str(:),AXR(:),AXI(:),BT(:),P(:,:),A(:,:),sai(:),
     1ALPHA(:,:),G(:,:),AT(:,:),ZIM(:),CHAR(:),
     1st2r(:),st2i(:),st3r(:),st3i(:),smmr(:),smmi(:),
     2smr(:),smi(:),sm3r(:),sm3i(:),ZM(:),AMULT(:)

      integer*4 irx,iry,irz
      common /randc/ irx, iry, irz
      irx=1
      iry=1000
      irz=99

      CONST=4.359828d0/0.5291772d0/0.5291772d0
      c1=1302.828d0
      pi=4.0d0*atan(1.0d0)
      it0=mclock()
      it00=-1
C
      WRITE(*,60000)
60000 FORMAT('      SPECTRA GENERATION BY TIME-PROPAGATION OF',/,
     1       '        CONTINUEOUS FOURIER TRANSFORMATION',/,
     1       '       PC version - in cartesian coordinates',/,
     2       '                 Scaling of masses!',/,/,
     2       '                For giant force fields',/,/,
     3       '            Petr Bour, Prague 2008',/,/)
C
      OPEN(17,FILE='FILE.X')
      READ(17,*)
      READ(17,*)NOAT
      CLOSE(17)

      NA=3*NOAT

      allocate(ZIM(NOAT),CHAR(NOAT),ZM(NA),AMULT(NA))
C
C     Read-in scalefactors:
      OPEN(15,FILE='FTRY.INP')
      READ(15,*)
      READ(15,*)NOSC
      DO 222 I=1,NOSC
222   READ(15,*) J,NM,(ix,J=1,NM)
1040  FORMAT(6G12.6)
C
C     Read-in atomic charges for FPC:
      READ(15,*) (CHAR(I),I=1,NOAT)
      CHATOT=0.0d0
      DO 77 I=1,NOAT
77    CHATOT=CHATOT+CHAR(I)
      WRITE(*,777)CHATOT,NOAT
777   FORMAT(' TOTAL CHARGE',F20.6,/I6,' atoms')
C
C     Read-in atomic masses:
      READ(15,*)
      READ(15,1040)(ZIM(I),I=1,NOAT)
      BM=0.0d0
      DO 102 IA=1,NOAT
102   BM=BM+ZIM(IA)
      DO 220 I=1,NOAT
      DO 220 J=1,3
220   ZM(3*(I-1)+J)=1/SQRT(ZIM(I))
      WRITE(*,778)BM
778   FORMAT(' TOTAL MASS  ',F20.6)
      CLOSE(15)
      WRITE(*,*)' FTRY.INP READ IN'
      write(6,6009)mclock()-it0
6009  format(' Time: ',i10,' msec')
C
      DO 201 I=1,NA
201   AMULT(I)=1.0d0
c
      allocate(j33(NA),a33(NA))
      write(6,*)'j33,a33 allocated'

      WRITE(*,*)' Making scaled mass-weighed FF from AFILE.FC'
c     number of elements:
      iel=0
      OPEN(20,FILE='AFILE.FC',STATUS='OLD',form='unformatted')
      OPEN(21,FILE='A3.SCR',form='unformatted')
      do 20 i=1,NA
      fi=ZM(i)*AMULT(i)*CONST
      read(20)ii,(j33(j),j=1,ii),(a33(j),j=1,ii)
      iel=iel+ii
20    write(21)ii,(j33(j),j=1,ii),(a33(j)*fi*ZM(j33(j))
     1*AMULT(j33(j)) ,j=1,ii)
      close(20)
      close(21)
      WRITE(*,*)' Cartesian FF rewritten into A3.SCR'
      write(6,6009)mclock()-it0
c
      write(6,6010)iel,dble(iel)/dble(NA*(NA+1)/2)*100.0d0
6010  format(i20,' elements (',f9.4, '%)')
      allocate(jbb(iel),abb(iel),mm(NA),nn(NA))
      write(6,*)'Space for the martix allocated'

      iel=0
      open(49,file='A3.SCR',form='unformatted')
      do 441 i=1,NA
      read(49)nn(i),(j33(ii),ii=1,nn(i)),(a33(ii),ii=1,nn(i))
      mm(i)=iel
      do 441 ii=1,nn(i)
      iel=iel+1
      jbb(iel)=j33(ii)   
441   abb(iel)=a33(ii)   
      close(49)
      write(6,4900)(iel*8+iel*4+NA*8)/1000
4900  format(' Matrix loaded into memory',i12,' kb')
      write(6,6009)mclock()-it0
c
c     load the matrix into memory - non-zero elements only
c     where the line starts   : mm
c     how many elements it has: nn
c
c     set the options:
      WMIN=100.0d0
      WMAX=4000.0d0
      DELW=10.0d0
      NP=4000
      nmd=10000
      dt=0.01d0
      LG=.true.
      LDW=.false.
      LEB=.false.
      LVEL=.false.
      LSEC=.true.
      LT3=.false.
      lrtr=.false.
      nhalf=1
      nsim=10
      nproc=1
      seed=0
      stemp=298.0d0
      inquire(file='NEW4FT.OPT',exist=lex)
      if(lex)then
       open(77,file='NEW4FT.OPT')
111    read(77,770,end=112,err=112)s5
770    format(a5)
       if(s5(1:4).eq.'WMIN')read(77,*)WMIN
       if(s5(1:4).eq.'WMAX')read(77,*)WMAX
       if(s5(1:4).eq.'DELW')read(77,*)DELW
       if(s5(1:3).eq.'NMD')read(77,*)NMD
       if(s5(1:2).eq.'DT')read(77,*)dt
       if(s5(1:2).eq.'NP')read(77,*)NP
       if(s5(1:2).eq.'LG')read(77,*)LG
       if(s5(1:3).eq.'LDW')read(77,*)LDW
       if(s5(1:3).eq.'LEB')read(77,*)LEB
       if(s5(1:3).eq.'LVE')read(77,*)LVEL
       if(s5(1:3).eq.'LSE')read(77,*)LSEC
       if(s5(1:3).eq.'LT3')read(77,*)LT3
       if(s5(1:4).eq.'NHAL')read(77,*)nhalf
       if(s5(1:4).eq.'NSIM')read(77,*)nsim
       if(s5(1:4).eq.'#PRO')read(77,*)nproc
       if(s5(1:4).eq.'RAND')read(77,*)seed
       if(s5(1:4).eq.'LRTR')read(77,*)lrtr
       if(s5(1:4).eq.'TEMP')read(77,*)stemp
       goto 111
112    close(77)
      endif
      irx=irx+seed
      iry=iry+seed+seed
      irz=irz+seed*5
      if(nhalf.eq.0)nhalf=1
      write(6,6001)WMIN,WMAX,DELW,NP,nmd,dt,nsim,nhalf,
     1omp_get_num_procs(),omp_get_max_threads(),nproc,seed
6001  format(' Spectral interval ',f10.1,' - ',f10.1,' cm-1',/,
     1       ' bandwidth         ',f10.2,' cm-1',/,
     2       ' number of spectral points ',i15,/,
     2       ' number of MD points       ',i15,/,
     2       ' integration time step     ',f15.3,/,
     2       ' number of simultanous spec',i15,/,
     2       ' pure propagation MD steps ',i15,/,
     2       ' available processors      ',i15,/,
     2       ' maximum of threads        ',i15,/,
     2       ' requested processors      ',i15,/,
     2       ' random seed               ',i15,/)

      allocate(sar(51*nproc*NP),st(4*NP),sff(NP),str(NA*nproc),
     1sti(NA*nproc),sai(51*nproc*NP))
      write(6,*)'Spectral space allocated'

      if(nproc.gt.omp_get_num_procs())
     1call report('not enough processors')
      if(nproc.gt.omp_get_max_threads())then
       write(6,*)'Number of threads will increase'
       call omp_set_num_threads(nproc)
      endif
      if(nproc.lt.omp_get_max_threads())then
       write(6,*)'Number of threads will decrease'
       call omp_set_num_threads(nproc)
      endif
      call omp_set_dynamic(0)
      call omp_set_nested(0)
c
c     precalculate parallel distribution:
      allocate(nstart(nproc),nend(nproc))
      nd=nsim/nproc
      nrest=nsim-nd*nproc
      nstart(1)=1
      nend(1)=nd
      if(nrest.ge.1)nend(1)=nend(1)+1
      do 7 ipr=2,nproc
      nstart(ipr)=nend(ipr-1)+1
      nend(ipr)=nstart(ipr)+nd-1
7     if(nrest.ge.ipr)nend(ipr)=nend(ipr)+1
      write(6,6007)
6007  format(' Parallel distribution:')
      do 72 ipr=1,nproc
72    write(6,60071)ipr,nstart(ipr),nend(ipr),nend(ipr)-nstart(ipr)+1
60071 format(i2,i5,' - ',i5,i10)


      write(6,*)' Velocity integration   : ',LVEL
      write(6,*)' Second derivs on fly   : ',LSEC
      write(6,*)' Third derivs on fly    : ',LT3
      write(6,*)' Double-ener integration: ',LDW
      write(6,*)' Debug DEBUG.TXT        : ',LEB
      write(6,*)' SQRT(Random)           : ',lrtr
      if(LG)then
       write(6,*)'Gaussian peaks'
       d=DELW/2.0d0/dsqrt(log(2.0d0))
      else
       write(6,*)'Lorentzian peaks'
       d=DELW/2.0d0
      endif
c
c     pre-calculate convolution points
      dd=(WMAX-WMIN)/dble(NP-1)
      if(LDW)then
       nff=0
       wj=-dd
       do 11 j=1,NP
       wj=wj+dd
       con=sf(d,LG,0.0d0,wj)
       if(con.gt.1.0d-3)then
        nff=nff+1
        sff(nff)=con
       else
        goto 12
       endif
11     continue
12     continue
       write(6,*)nff,' convolution points'
      endif
c
c     read dipole derivatives
      allocate(P(NA,3),A(NA,3),ALPHA(NA,9),G(NA,9),AT(NA,27))
      call READTEN(NA,P,A,NOAT)
c
c     read POLARIZABILITY derivatives
      call READTTT(NA,ALPHA,AT,G,NOAT)
c
c     Make mass-weighed derivatives for the propagation
      FD=dsqrt(388.9219d0)
      FC=131.14d-8
      do 10 i=1,NA
      do 101 ix=1,3
      P(i,ix)=P(i,ix)*ZM(i)*FD
101   A(i,ix)=A(i,ix)*ZM(i)*FC
      do 103 ix=1,9
      ALPHA(i,ix)=ALPHA(i,ix)*ZM(i)
103   G(i,ix)=G(i,ix)*ZM(i)
      do 10 ix=1,27
10    AT(i,ix)=AT(i,ix)*ZM(i)
c
c     create bit tensor
c       1 ...  3 P
c       4 ...  6 A
c       7 ... 15 ALPHA
c      16 ... 24 G
c      25 ... 51 A

      allocate(BT(NA*51))
      do 78 i=1,NA
      do 781 ix=1,3
      BT(i+(ix   -1)*NA)=P(i,ix)
781   BT(i+(ix+ 3-1)*NA)=A(i,ix)
      do 782 ix=1,9
      BT(i+(ix+ 6-1)*NA)=ALPHA(i,ix)
782   BT(i+(ix+15-1)*NA)=G(i,ix)
      do 78  ix=1,27
78    BT(i+(ix+24-1)*NA)=AT(i,ix)


      write(6,6009)mclock()-it0
c
c     initialize spectra   
      do 2994 i=1,4*NP
2994  st(i)=0.0d0
c
c     initialize random function
      do 2993 i=1,13*seed
2993  if(frandom(0).gt.100.0d0)write(6,*)'blbost'
c
      allocate(AXI(NA*nproc),AXR(NA*nproc),st3r(NA*nproc),
     1st3i(NA*nproc),smmr(NA*nproc),
     1smmi(NA*nproc),smr(NA*nproc),smi(NA*nproc),sm3r(NA*nproc),
     2st2r(NA*nproc),st2i(NA*nproc),sm3i(NA*nproc))
      write(6,*)'AXI, AXR allocated'
      it00=mclock()
c
C$OMP PARALLEL DEFAULT(SHARED)
C$OMP+ PRIVATE(ipr,is,io,i,ANORM,rtr,t,it,ii,inn,jn,ann,ui,ur,uix,urx,
C$OMP+ ix,wi,ei,ce,se,tr,ti,j,ik,di,ri,SAL0,SAL1,SAG0,SAG1,SA1,
C$OMP+ i9,i8,i7,iti,itj,k,l,itkl)

c     pseudo-private variables:
c     51*nproc*NP:
c     sai,sar  
c     nproc*NA:
c     sti,str,st3r,st3i,st2r,st2i,sm3r,sm3i,smmr,smmi,smr,smi,AXR,AXI

      ipr=OMP_GET_THREAD_NUM()+1
      write(6,*)' Integration started, thread ',ipr
c     loop over spectra:
      do 899 is=nstart(ipr),nend(ipr)
      if(ipr.eq.1.and.is.eq.2)then
       it00=mclock()-it00
       write(6,6005)it00,it00*(nend(ipr)-nstart(ipr)+1)
6005   format('One spectrum time:',i10,' total estimate:',i11)
      endif

c     initialize vectors u m alpha g A:
      io=(ipr-1)*51*NP
      do 8991 i=io+1,io+51*NP
      sai(i)=0.0d0
8991  sar(i)=0.0d0
c
c     Initialize trial vector-real and imaginary part,str, sti
      do 2492 i=1,(ipr-1)*(131+13*is)
2492  if(frandom(0).gt.100.0d0+dble(i-i))write(6,*)'blbost'
      io=(ipr-1)*NA
      do 2992 i=io+1,io+NA
      sti(i)=frandom(0)
      if(sti(i).gt.100.0d0)write(6,*)'blbost'
2992  sti(i)=0.0d0

      ANORM=0.0d0
      do 89921 i=io+1,io+NA
      rtr=frandom(0)-0.5d0
      if(lrtr)then
       str(i)=dsqrt(dabs(rtr))
       if(rtr.lt.0.0d0)str(i)=-str(i)
      else
       str(i)=rtr
      endif
89921 ANORM=ANORM+str(i)**2
      ANORM=1.0d0/dsqrt(ANORM)
      do 81 i=io+1,io+NA
81    str(i)=str(i)*ANORM
c  
c     Initialize second,third derivatives and old vectors
      do 8993 i=io+1,io+NA
      st3r(i)=0.0d0
      st3i(i)=0.0d0
      st2r(i)=0.0d0
      st2i(i)=0.0d0
      sm3r(i)=str(i) 
      sm3i(i)=sti(i) 
      smmr(i)=str(i) 
      smmi(i)=sti(i) 
      smr(i)=str(i) 
8993  smi(i)=sti(i) 
c
c     do the propagation:
      t=0.0d0
      do 3 it=1,nmd
      io=(ipr-1)*NA
      t=t+dt
c
c     mutliply st by FF, store in AX:
      do 4993 i=io+1,io+NA
      AXR(i)=0.0d0
4993  AXI(i)=0.0d0

      DO 4 I=1,NA
      DO 4 ii=1,nn(I)
      inn=mm(I)+ii
      jn=jbb(inn)
      ann=abb(inn)
      AXR(io+jn)=AXR(io+jn)+ann*str(io+I)
      AXR(io+I )=AXR(io+I )+ann*str(io+jn)
      AXI(io+jn)=AXI(io+jn)+ann*sti(io+I)
4     AXI(io+I )=AXI(io+I )+ann*sti(io+jn)

      do 5 i=io+1,io+NA
      str(i)=str(i)-AXI(i)*dt
5     sti(i)=sti(i)+AXR(i)*dt

      if(LSEC)then
       do 52 i=io+1,io+NA
       str(i)=str(i)+st2r(i)
52     sti(i)=sti(i)+st2i(i)
      endif
      if(LT3)then
       do 54 i=io+1,io+NA
       str(i)=str(i)+st3r(i)
54     sti(i)=sti(i)+st3i(i)
      endif
c
c     normalize
      ANORM=0.0d0
      do 4995 i=io+1,io+NA
4995  ANORM=ANORM+str(i)**2+sti(i)**2
      if(dabs(ANORM-1.0d0).gt.0.3d0)write(6,*)it,'norm too high',ANORM
      ANORM=1.0d0/dsqrt(ANORM)
      do 8995 i=io+1,io+NA
      sti(i)=sti(i)*ANORM
8995  str(i)=str(i)*ANORM
c
c     calculate second derivatives*dt^2*0.5
c     calculate third derivatives*dt^3/6
      do 84 i=io+1,io+NA
      st3r(i)=(3.0d0*(smmr(i)-smr(i))+str(i)-sm3r(i))/6.0d0
      st3i(i)=(3.0d0*(smmi(i)-smi(i))+sti(i)-sm3i(i))/6.0d0
      st2r(i)=(str(i)+smmr(i))*0.5d0-smr(i)
84    st2i(i)=(sti(i)+smmi(i))*0.5d0-smi(i)
c
c     remember this and previous vectors
      do 85 i=io+1,io+NA
      sm3r(i)=smmr(i) 
      sm3i(i)=smmi(i) 
      smmr(i)=smr(i) 
      smmi(i)=smi(i) 
      smr(i)=str(i) 
85    smi(i)=sti(i) 

c     time-dependent dipole moments and autocorrelation function: 
c          ur ui:
c       1 ...  3 u
c       4 ...  6 m
c       7 ... 15 al
c      16 ... 24 g
c      25 ... 51 a

      do 71 ix=1,51
      urx=0.0d0
      uix=0.0d0
      ii=(ix-1)*NA
      do 73 i=1,NA
      urx=urx+BT(i+ii)*str(io+i)
73    uix=uix+BT(i+ii)*sti(io+i)
      ur(ix)=urx*dt
71    ui(ix)=uix*dt

c     accumulate FT - double w-integration to disperse it:
      wi=WMIN-dd
      do 8 i=1,NP
      wi=wi+dd
      ei=(wi/c1)**2
      ce=cos(ei*t)
      se=sin(ei*t)
      
      do 74 ix=1,51
      tr(ix)=ce*ur(ix)+se*ui(ix)
74    ti(ix)=ce*ui(ix)-se*ur(ix)

      io=(ipr-1)*51*NP
      if(LDW)then
       do 82 j=i,min(NP,i+nff-1)
82     call sadd(sff(j-i+1),tr,ti,sar,sai,io+51*(j-1))
       do 76 j=max(1,i-nff+1),i-1
76     call sadd(sff(i-j+1),tr,ti,sar,sai,io+51*(j-1))
      else
       call sadd(1.0d0,tr,ti,sar,sai,io+51*(i-1))
      endif
8     continue
c
3     continue
c     end of propagation
c
c     add the spectrum to the total:
C$OMP CRITICAL
      wi=WMIN-dd
      do 8997 i=1,NP
      wi=wi+dd
c     1 abs 
c     2 vcd 
c     3 Raman-180
c     4 ROA-ICP-180

      ik=51*((ipr-1)*NP+i-1)
      di=sar(ik+1)*sar(ik+1)+sar(ik+2)*sar(ik+2)+sar(ik+3)*sar(ik+3)
      ri=sar(ik+1)*sar(ik+4)+sar(ik+2)*sar(ik+5)+sar(ik+3)*sar(ik+6)
      SAL0=0.0d0
      SAL1=0.0d0
      SAG0=0.0d0
      SAG1=0.0d0
      SA1=0.0d0
      i9=ik+6
      i8=ik+15
      i7=ik+24
      DO 6 ii=1,3
      DO 6 J=1,3
      iti=ii+3*(ii-1)
      itj=J+3*(J-1)
      io=ii+3*(J-1)
c     SAL0=SAL0+ALPHA(IQ,I,I)*    ALPHA(IQ,J,J)
      SAL0=SAL0+sar(i9+iti)*sar(i9+itj)+
     1          sai(i9+iti)*sai(i9+itj)
c     SAL1=SAL1+ALPHA(IQ,I,J)*    ALPHA(IQ,I,J)
      SAL1=SAL1+sar(i9+io )*sar(i9+io )+
     1          sai(i9+io )*sai(i9+io )
c     SAG0=SAG0+ALPHA(IQ,I,I)*    G(IQ,J,J)
      SAG0=SAG0+sar(i9+iti )*sar(i8+itj )+
     1          sai(i9+iti )*sai(i8+itj )
c     SAG1=SAG1+ALPHA(IQ,I,J)*    G(IQ,I,J)
      SAG1=SAG1+sar(i9+io )*sar(i8+io )+
     1          sai(i9+io )*sai(i8+io )
c     DO 6 K=1,3
c     DO 6 L=1,3
c     SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0
c     itkl=K+3*(L-1)+9*(J-1)
c     SA1=SA1+EPS(ii,K,L)*(sar(i9+io)*sar(i7+itkl)+
c    1                     sai(i9+io)*sai(i7+itkl))/3.0d0
      K=ii+1
      if(K.gt.3)K=1
      L=K+1
      if(L.gt.3)L=1
      itkl=i7+K+3*(L-1)+9*(J-1)
      SA1=SA1+(sar(i9+io)*sar(itkl)+sai(i9+io)*sai(itkl))/3.0d0
      itkl=i7+L+3*(K-1)+9*(J-1)
6     SA1=SA1-(sar(i9+io)*sar(itkl)+sai(i9+io)*sai(itkl))/3.0d0
c
      st(     i)=st(     i)+di
      st(  NP+i)=st(  NP+i)+ri
      st(2*NP+i)=st(2*NP+i)+7.0d0*SAL1+SAL0
8997  st(3*NP+i)=st(3*NP+i)+8.0d0*(3.0d0*SAG1-SAG0+SA1)
C$OMP END CRITICAL

      if(is.eq.nend(ipr))then
       write(6,91)ipr,nstart(ipr),nend(ipr)
91     format('Thread',i3,' spectra ',i4,' - ',i4,' added.')
      endif
899   continue
C$OMP END PARALLEL 
c     
      if(LDW)then
       c2=1.0d0/dble(NA)*2.0d0/pi*2.0d0*d/c1**2/dble(nsim)
      else
       c2=1.0d0/dble(NA)*2.0d0/pi/dble(nmd)/dt/dble(nsim)
      endif
      open(91,file='D.PRN')
      open(92,file='R.PRN')
      open(93,file='RAM.PRN')
      open(94,file='ROA.PRN')
      wi=WMIN-dd
      do 64 i=1,NP
      wi=wi+dd
      s1=st(     i)
      s2=st(  NP+i)
      s3=st(2*NP+i)
      s4=st(3*NP+i)
      if(LDW)then
       s1=s1*108.7d0*wi**2*c2
       s2=s2*435.0d0*wi**3*c2
       s3=s3*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp)))
       s4=s4*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp)))
      else
       s1=s1*108.7d0*wi*c2
       s2=s2*435.0d0*wi**2*c2
       s3=s3*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp)))
       s4=s4*dsqrt(dabs(wi))/(1.0d0-exp(-(1.44d0*wi/stemp)))
      endif
      write(91,494)wi,s1
      write(92,494)wi,s2
      write(93,494)wi,s3
64    write(94,494)wi,s4
494   format(g15.6,' ',g15.6)
      close(91)
      close(92)
      close(93)
      close(94)
      write(6,*)' Spectra written'
      write(6,6009)mclock()-it0
      END
c
      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      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.1000.0d0)then
        sf=1.0d0/d/(dd+1.0d0)/pi
       else
        sf=0.0d0
       endif
      endif
      return
      end

      SUBROUTINE READTEN(N0,P,A,NAT)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION P(N0,3),A(N0,3)
      open(15,FILE='FILE.TEN')
      READ (15,*) NAT
      WRITE(*,*)NAT,' atoms in FILE.TEN'
      DO 10 L=1,NAT
      DO 10 J=1,3                                                     
10    READ (15,*) (P(3*(L-1)+J,I),I=1,3)
      WRITE(6,*)' Electric dipole derivatives read in '
      DO 220 L=1,NAT                                                 
      DO 220 J=1,3                                                    
220   READ (15,*) (A(3*(L-1)+J,I),I=1,3)
      WRITE(6,*)' MAgnetic dipole derivatives read in '
      DO 230 L=1,NAT                                                 
      DO 230 J=1,3                                                    
230   READ (15,*) 
      DO 100 L=1,NAT                                                 
      DO 100 J=1,3                                                    
100   READ (15,*)
      close(15)
      RETURN
      END

      subroutine sadd(con,tr,ti,sar,sai,i33)
      implicit none
      real*8 sar(*),sai(*),con,tr(*),ti(*)
      integer*4 i33,ix
      do 1 ix=1,51
      sar(i33+ix)=sar(i33+ix)+con*tr(ix)
1     sai(i33+ix)=sai(i33+ix)+con*ti(ix)
      return
      end
C
      SUBROUTINE READTTT(M,ALPHA,A,G,NAT)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION ALPHA(M,9),G(M,9),A(M,27)
      logical lex
      inquire(file='FILE.TTT',exist=lex)
      if(.not.lex)return
      open(15,file='FILE.TTT')
      READ(15,*)
      READ(15,*)NATC
      if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file')
      READ(15,*)
      READ(15,*)
      DO 1 I=1,3
      READ(15,*)
      DO 1 L=1,NAT
      DO 1 IX=1,3
      IIND=3*(L-1)+IX
1     READ(15,*)(ALPHA(IIND,I+3*(J-1)),J=1,2),
     1(ALPHA(IIND,I+3*(J-1)),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 2 I=1,3
      READ(15,*)
      DO 2 L=1,NAT
      DO 2 IX=1,3
      IIND=3*(L-1)+IX
2     READ(15,*)(G(IIND,I+3*(J-1)),J=1,2),
     1(G(IIND,I+3*(J-1)),J=1,3)
      READ(15,*)
      READ(15,*)
      DO 3 I=1,3
      DO 3 J=1,3
      READ(15,*)
      DO 3 L=1,NAT
      DO 3 IX=1,3
      IIND=3*(L-1)+IX
3     READ(15,*)(A(IIND,I+3*(J-1)+9*(K-1)),K=1,2),
     1(A(IIND,I+3*(J-1)+9*(K-1)),K=1,3)
      CLOSE(15)
      RETURN
      END
      
      function frandom(i)
      implicit none
      real*8 frandom
      integer*4 i
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