program powersp
      implicit none
c     spectrum from autocorrelation function,
c     form tinker MD positions and velocities
      integer*4 nat0,i,ie,is,ie2,np0,ix,ia,nat,ii,is2,
     1it
      parameter (nat0=2000,np0=10000)
      character*80 s80,fn,fl,lines
      character*1 ok
      logical lex
      real*8 v(3*nat0),v1(3*nat0),st(np0),sum,t,dt,w,dw,wmax,
     1sums,sumc,tp,sol,wmin

      write(6,*)' File name (without extension):'
      read(5,500)fn
 500  format(a80)
      write(6,*)' Coordinates or velocity (c/v):'
      read(5,509)ok
509   format(a1)
      if(ok.eq.'c')ok='C'
      call trim(fn,is,ie)
      i=1
1     write(s80,609)i
609   format(i6)
      if(i.lt.100)s80(4:4)='0'
      if(i.lt.10)s80(5:5)='0'
      call trim(s80,is2,ie2)
      if(ok.eq.'C')then
       fl=fn(is:ie)//'.'//s80(is2:ie2)
      else
       fl=fn(is:ie)//'.'//s80(is2:ie2)//'v'
      endif
      write(6,506)fl(1:50)
506   format(a50)
      inquire(file=fl,exist=lex)
      if(lex)then
       open(8,file=fl)
       read(8,801)nat
801    format(i6)

       if(nat.gt.nat0)then
        write(6,*)'Too many atoms'
        stop
       endif

       do 3 ia=1,nat
       read(8,500)lines
       lines(1:11)='           '
3      read(lines,*)(v(ix+3*(ia-1)),ix=1,3)

       if(i.eq.1)then
	do 2 ii=1,3*nat
2       v1(ii)=v(ii)
       endif
c
c      Calculate the correlated property:
       sum=0.0d0
       do 4 ii=1,3*nat
4      sum=sum+v1(ii)*v(ii)
       st(i)=sum

       close(8)
       i=i+1
       if(i.le.np0)then
        goto 1
       else
        write(6,*)'Number of points exceeded'
	goto 1001
       endif
      else
       i=i-1
      endif
1001  write(6,6001)i
 6001 format(i7,' points read')
      open(8,file='T.PRN')
      do 5 ii=1,i
5     write(8,804)ii,st(ii)
804   format(i7,f22.6)
      close(8)
      write(6,*)'T.PRN written'

      write(6,*)'delta t (fs):'
      read(5,*)dt

      dt=dt*1.0d-15
      tp=2.0d0*3.141592653589d0
      sol=3.0d10
      dw=tp/dt/dble(i)
      write(6,6004)dw/tp/sol,dw/tp/sol*dble(i)
6004  format(' dw = ',f12.2,' cm-1, wmax = ',f15.2,' cm-1',/,/,
     1' Input new dw, wmin, wmax:')
      read(5,*)dw,wmin,wmax
      dw=dw*tp*sol
      wmax=wmax*tp*sol
c
      open(8,file='W.PRN')

      w=wmin-dw
6     w=w+dw
      sums=0.0d0
      sumc=0.0d0
      t=-dt/2
      do 62 it=1,i
      t=t+dt
      sums=sums+sin(t*w)*st(it)
62    sumc=sumc+cos(t*w)*st(it)
      write(8,805)w/2.0d0/3.141592653589d0/3.0d10,
     1(sums**2+sumc**2)/dble(i)
805   format(f12.2,f25.6)
      if(w+dw.le.wmax)goto 6

      close(8)
      write(6,*)'W.PRN written'
      end
      
      subroutine trim(fn,is,ie)
      implicit none
      character*80 fn
      integer*4 is,ie
      do 1 is=1,len(fn)
1     if(fn(is:is).ne.' ')goto 2
2     do 3 ie=len(fn),1,-1
3     if(fn(ie:ie).ne.' ')goto 4
4     return
      end