1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  2000  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #################################################################
9c     ##                                                             ##
10c     ##  program spectrum  --  power spectrum from autocorrelation  ##
11c     ##                                                             ##
12c     #################################################################
13c
14c
15c     "spectrum" computes a power spectrum over a wavelength range
16c     from the velocity autocorrelation as a function of time
17c
18c
19      program spectrum
20      use files
21      use iounit
22      use math
23      use units
24      implicit none
25      integer i,k
26      integer next,nsamp
27      integer ivel,nvel
28      integer maxvel
29      integer maxfreq
30      integer freeunit
31      real*8 factor,aver
32      real*8 norm,step
33      real*8 freq,time
34      real*8, allocatable :: vel(:)
35      real*8, allocatable :: intense(:)
36      logical exist,done
37      character*240 velfile
38      character*240 record
39      character*240 string
40c
41c
42c     perform the standard initialization functions
43c
44      call initial
45c
46c     try to get a filename from the command line arguments
47c
48      call nextarg (velfile,exist)
49      if (exist) then
50         call basefile (velfile)
51         call suffix (velfile,'vel','old')
52         inquire (file=velfile,exist=exist)
53      end if
54c
55c     ask for the velocity autocorrelation data filename
56c
57      do while (.not. exist)
58         write (iout,10)
59   10    format (/,' Enter Name of Velocity Autocorrelation',
60     &              ' File :  ',$)
61         read (input,20)  velfile
62   20    format (a240)
63         call basefile (velfile)
64         call suffix (velfile,'vel','old')
65         inquire (file=velfile,exist=exist)
66      end do
67c
68c     get the time step between autocorrelation data points
69c
70      step = -1.0d0
71      call nextarg (string,exist)
72      if (exist)  read (string,*,err=30,end=30)  step
73   30 continue
74      if (step .le. 0.0d0) then
75         write (iout,40)
76   40    format (/,' Enter Time Step for Autocorrelation Data',
77     &              ' [1.0 fs] :  ',$)
78         read (input,50)  step
79   50    format (f20.0)
80      end if
81      if (step .le. 0.0d0)  step = 1.0d0
82      step = 0.001d0 * step
83c
84c     open the velocity autocorrelation data file
85c
86      ivel = freeunit ()
87      open (unit=ivel,file=velfile,status='old')
88      rewind (unit=ivel)
89c
90c     read through file headers to the start of the data
91c
92      done = .false.
93      do while (.not. done)
94         read (ivel,60)  record
95   60    format (a240)
96         next = 1
97         call getword (record,string,next)
98         if (string(1:10) .eq. 'Separation') then
99            done = .true.
100            read (ivel,70)
101   70       format ()
102         end if
103      end do
104c
105c     perform dynamic allocation of some local arrays
106c
107      maxvel = 100000
108      maxfreq = 5000
109      allocate (vel(0:maxvel))
110      allocate (intense(maxfreq))
111c
112c     read the velocity autocorrelation as a function of time
113c
114      do i = 1, maxvel
115         read (ivel,80,err=90,end=90)  record
116   80    format (a240)
117         read (record,*)  k,nsamp,aver,norm
118         nvel = k
119         vel(k) = norm
120      end do
121   90 continue
122c
123c     compute the power spectrum via discrete Fourier transform
124c
125      factor = 2.0d0 * pi * lightspd
126      do i = 1, maxfreq
127         freq = factor * dble(i)
128         intense(i) = 0.0d0
129         do k = 0, nvel
130            time = step * dble(k)
131            intense(i) = intense(i) + vel(k)*cos(freq*time)
132         end do
133         intense(i) = 1000.0d0 * step * intense(i)
134      end do
135c
136c     print the power spectrum intensity at each wavelength
137c
138      write (iout,100)
139  100 format (/,' Power Spectrum from Velocity Autocorrelation :',
140     &        //,4x,'Frequency (cm-1)',10x,'Intensity',/)
141      do i = 1, maxfreq
142         write (iout,110)  dble(i),intense(i)
143  110    format (3x,f12.2,8x,f16.6)
144      end do
145c
146c     perform deallocation of some local arrays
147c
148      deallocate (vel)
149      deallocate (intense)
150c
151c     perform any final tasks before program exit
152c
153      call final
154      end
155