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