program ftcpmdx
      character*80 fn
      parameter (np0=10000,nat0=200)
      dimension ss(np0),sc(np0),x(3*nat0),x0(3*nat0),xm(3*nat0),
     1ise(nat0),ls(nat0)
      logical sel,ls
      if(iargc().eq.0)then
       write(6,*)'Filename:'
       read(5,*)fn
      else
       call getarg(1,fn)
      endif
      inquire(file='SEL.PAR',exist=sel)
      if(sel)then
       write(6,*)'Selected atoms'
       open(67,file='SEL.PAR')
       read(67,*)nsel
       read(67,*)(ise(io),io=1,nsel)
       write(6,*)(ise(io),io=1,nsel)
       close(67)
      endif
      wmin=0.0
      wmax=4000.0
      np=4000
      do 4 i=1,np
      ss(i)=0.0
4     sc(i)=0.0
      t=0.0
c     time step in atomic units:
      dt=200.0
      dt=dt/50.0
c     time step in seconds:
      dt=dt*2.418e-17
      open(9,file=fn)
      open(91,file='time.prn')
      ii=0
1     read(9,*,end=99,err=99)
      read(9,*)nat
      if(nat.gt.nat0)stop
      do 5 ia=1,nat
5     read(9,*)idumm,(x(ix+3*(ia-1)),ix=1,3)
      if(ii.eq.0)then
       do 8 ia=1,nat
8      ls(ia)=.true.
       if(sel)then
        do 9 ia=1,nat
9       ls(ia)=.false.
        do 10 ia=1,nsel
10      ls(ise(ia))=.true.
       endif
      endif
      ii=ii+1
      t=t+dt
      if(ii.eq.1)then
       do 6 ia=1,3*nat
6      x0(ia)=x(ia)
      endif
      a=0.0
      do 7 ia=1,nat
      do 7 ix=1,3
      ind=ix+3*(ia-1)
7     if(ls(ia))a=a+(x(ind)-xm(ind))*x0(ind)
      do 62 ia=1,3*nat
62    xm(ia)=x(ia)
      if(ii.gt.2)then
       write(91,*)t,a
       do 3 i=1,np
c      frequency in cm-1:
       w=wmin+real(i-1)/real(np-1)*(wmax-wmin)
c      frequency in au:
       w=2.0*3.141592*3.0e10*w
       ss(i)=ss(i)+a*sin(w*t)
3      sc(i)=sc(i)+a*cos(w*t)
      endif
      goto 1
99    close(9)
      close(91)
      write(6,*)ii,' points'
      open(9,file='S.PRN')
      do 2 i=1,np
      w=wmin+real(i-1)/real(np-1)*(wmax-wmin)
2     write(9,900)w,ss(i)**2+sc(i)**2
900   format(f10.2,f25.9)
      close(9)
      stop
      end