1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21  program vibrational
22    use batch_oct_m
23    use command_line_oct_m
24    use geometry_oct_m
25    use global_oct_m
26    use io_oct_m
27    use messages_oct_m
28    use namespace_oct_m
29    use parser_oct_m
30    use profiling_oct_m
31    use simul_box_oct_m
32    use space_oct_m
33    use spectrum_oct_m
34    use unit_oct_m
35    use unit_system_oct_m
36
37    implicit none
38
39    integer :: iunit, ierr, ii, jj, iter, read_iter, ntime, nvaf, nvel, ivel
40    FLOAT, allocatable :: vaf(:), time(:), velocities(:, :), ftvaf(:)
41    type(geometry_t)  :: geo
42    type(space_t)     :: space
43    type(simul_box_t) :: sb
44    type(spectrum_t) :: spectrum
45    type(batch_t) :: vafb, ftvafb
46    FLOAT :: ww, curtime, vaftime, deltat
47    integer :: ifreq, max_freq
48    integer :: skip
49
50    ! Initialize stuff
51    call global_init(is_serial = .true.)
52
53    call getopt_init(ierr)
54    call getopt_end()
55
56    call parser_init()
57
58    call messages_init()
59
60    call io_init()
61
62    call unit_system_init(global_namespace)
63
64    call spectrum_init(spectrum, global_namespace, &
65      default_energy_step = units_to_atomic(unit_invcm, CNST(0.2)), &
66      default_max_energy  = units_to_atomic(unit_invcm, CNST(5000.0)))
67
68    !%Variable VibrationalSpectrumTimeStepFactor
69    !%Type integer
70    !%Default 10
71    !%Section Utilities::oct-vibrational_spectrum
72    !%Description
73    !% In the calculation of the vibrational spectrum, it is not necessary
74    !% to read the velocity at every time step. This variable controls
75    !% the integer factor between the simulation time step and the
76    !% time step used to calculate the vibrational spectrum.
77    !%End
78
79    call messages_obsolete_variable(global_namespace, 'PropagationSpectrumTimeStepFactor', 'VibrationalSpectrumTimeStepFactor')
80    call parse_variable(global_namespace, 'VibrationalSpectrumTimeStepFactor', 10, skip)
81    if(skip <= 0) call messages_input_error(global_namespace, 'VibrationalSpectrumTimeStepFactor')
82
83    max_freq = spectrum_nenergy_steps(spectrum)
84
85    if (spectrum%end_time < M_ZERO) spectrum%end_time = huge(spectrum%end_time)
86
87    call space_init(space, global_namespace)
88    call geometry_init(geo, global_namespace, space)
89    call simul_box_init(sb, global_namespace, geo, space)
90
91    ! Opens the coordinates files.
92    iunit = io_open('td.general/coordinates', global_namespace, action='read')
93
94    call io_skip_header(iunit)
95
96    ntime = 1
97    iter = 1
98
99    ! check the number of time steps we will read
100    do
101      read(unit = iunit, iostat = ierr, fmt = *) read_iter, curtime, &
102        ((geo%atom(ii)%x(jj), jj = 1, 3), ii = 1, geo%natoms), &
103        ((geo%atom(ii)%v(jj), jj = 1, 3), ii = 1, geo%natoms)
104
105      curtime = units_to_atomic(units_out%time, curtime)
106
107      if(ierr /= 0 .or. curtime >= spectrum%end_time) then
108        iter = iter - 1	! last iteration is not valid
109        ntime = ntime - 1
110        exit
111      end if
112
113      if(iter /= read_iter + 1) then
114        call messages_write("Error while reading file 'td.general/coordinates',", new_line = .true.)
115        call messages_write('expected iteration ')
116        call messages_write(iter - 1)
117        call messages_write(', got iteration ')
118        call messages_write(read_iter)
119        call messages_write('.')
120        call messages_fatal()
121      end if
122
123      ! ntime counts how many steps are gonna be used
124      if (curtime >= spectrum%start_time .and. mod(iter, skip) == 0) ntime = ntime + 1
125
126      iter = iter + 1 !counts number of timesteps (with time larger than zero up to SpecEndTime)
127    end do
128
129    call io_close(iunit)
130
131    nvel = geo%natoms*space%dim
132
133    SAFE_ALLOCATE(time(1:ntime))
134    SAFE_ALLOCATE(velocities(1:nvel, 1:ntime))
135
136    ! Opens the coordinates files.
137    iunit = io_open('td.general/coordinates', global_namespace, action='read')
138
139    call io_skip_header(iunit)
140
141    ntime = 1
142    iter = 1
143
144    do
145      read(unit = iunit, iostat = ierr, fmt = *) read_iter, curtime, &
146        ((geo%atom(ii)%x(jj), jj = 1, 3), ii = 1, geo%natoms), &
147        ((geo%atom(ii)%v(jj), jj = 1, 3), ii = 1, geo%natoms)
148
149      curtime = units_to_atomic(units_out%time, curtime)
150
151      if(ierr /= 0 .or. curtime >= spectrum%end_time) then
152        iter = iter - 1	! last iteration is not valid
153        ntime = ntime - 1
154        exit
155      end if
156
157      ASSERT(iter == read_iter + 1)
158
159      if (curtime >= spectrum%start_time .and. mod(iter, skip) == 0) then
160
161        time(ntime) = curtime
162        ivel = 1
163        do ii = 1, geo%natoms
164          do jj = 1, space%dim
165            velocities(ivel, ntime) = units_to_atomic(units_out%velocity, geo%atom(ii)%v(jj))
166            ivel = ivel + 1
167          end do
168        end do
169
170        ntime = ntime + 1
171      end if
172
173      iter = iter + 1 !counts number of timesteps (with time larger than zero up to SpecEndTime)
174    end do
175
176    call io_close(iunit)
177
178    deltat = time(2) - time(1)
179
180    !%Variable VibrationalSpectrumTime
181    !%Type integer
182    !%Section Utilities::oct-vibrational_spectrum
183    !%Description
184    !% This variable controls the maximum time for the calculation of
185    !% the velocity autocorrelation function. The default is the total
186    !% propagation time.
187    !%End
188    call parse_variable(global_namespace, 'VibrationalSpectrumTime', ntime*deltat, vaftime)
189
190    nvaf = int(vaftime/deltat)
191
192    SAFE_ALLOCATE(vaf(1:ntime))
193
194    call messages_write('Time step = ')
195    call messages_write(deltat, units = units_out%time)
196    call messages_info()
197
198    call messages_new_line()
199    call messages_write('Calculating the velocity autocorrelation function')
200    call messages_info()
201
202    call calculate_vaf(vaf)
203
204   !print the vaf
205    iunit = io_open('td.general/velocity_autocorrelation', global_namespace, action='write')
206
207800 FORMAT(80('#'))
208    write(unit = iunit, iostat = ierr, fmt = 800)
209    write(unit = iunit, iostat = ierr, fmt = '(8a)')  '# HEADER'
210    write(unit = iunit, iostat = ierr, fmt = '(a,4x,a6,a7,a1,10x,a10)') &
211      '#       Iter', 'time [',units_out%time%abbrev,']', 'VAF [a.u.]'
212    write(unit = iunit, iostat = ierr, fmt = 800)
213
214    do jj = 1, nvaf
215      write(unit = iunit, iostat = ierr, fmt = *) jj, units_from_atomic(units_out%time, (jj - 1)*deltat), vaf(jj)
216    end do
217
218    call io_close(iunit)
219
220
221    SAFE_ALLOCATE(ftvaf(1:max_freq))
222
223    ftvaf = M_ONE
224
225    call batch_init(vafb, vaf)
226
227    call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, 1, nvaf, M_ZERO, deltat, vafb)
228
229    call batch_init(ftvafb, ftvaf)
230
231    call spectrum_fourier_transform(spectrum%method, SPECTRUM_TRANSFORM_COS, spectrum%noise, &
232      1, nvaf, M_ZERO, deltat, vafb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, ftvafb)
233
234    call vafb%end()
235    call ftvafb%end()
236
237
238    !and print the spectrum
239    iunit = io_open('td.general/vibrational_spectrum', global_namespace, action='write')
240
241    write(unit = iunit, iostat = ierr, fmt = 800)
242    write(unit = iunit, iostat = ierr, fmt = '(8a)')  '# HEADER'
243    write(unit = iunit, iostat = ierr, fmt = '(a17,6x,a15)') '#   Energy [1/cm]', 'Spectrum [a.u.]'
244    write(unit = iunit, iostat = ierr, fmt = 800 )
245
246    do ifreq = 1, max_freq
247      ww = spectrum%energy_step*(ifreq - 1) + spectrum%min_energy
248      write(unit = iunit, iostat = ierr, fmt = '(2e20.10)') units_from_atomic(unit_invcm, ww), ftvaf(ifreq)
249    end do
250
251    call io_close(iunit)
252
253    SAFE_DEALLOCATE_A(vaf)
254
255    SAFE_DEALLOCATE_A(ftvaf)
256
257    call simul_box_end(sb)
258    call geometry_end(geo)
259    call space_end(space)
260
261    SAFE_DEALLOCATE_A(time)
262
263    call io_end()
264    call messages_end()
265
266    call parser_end()
267    call global_end()
268
269  contains
270
271    subroutine calculate_vaf(vaf)
272      FLOAT, intent(out) :: vaf(:)
273      integer :: itm, itn
274
275      PUSH_SUB(calculate_vaf)
276
277      write (message(1), '(a)') "Read velocities from '"// &
278        trim(io_workpath('td.general/coordinates', global_namespace))//"'"
279      call messages_info(1)
280
281      !calculating the vaf, formula from
282      !
283      ! http://www.timteatro.net/2010/09/29/velocity-autocorrelation-and-vibrational-spectrum-calculation/
284
285      vaf = M_ZERO
286
287      do itm = 1, ntime
288        vaf(itm) = M_ZERO
289
290        do itn = 1, ntime - itm + 1
291
292          do ivel = 1, nvel
293            vaf(itm) = vaf(itm) + velocities(ivel, itm + itn - 1)*velocities(ivel, itn)
294          end do
295        end do
296
297        vaf(itm) = vaf(itm)/TOFLOAT(ntime - itm + 1)
298
299      end do
300
301      do itm = 2, ntime
302        vaf(itm) = vaf(itm)/vaf(1)
303      end do
304      vaf(1) = M_ONE
305
306      POP_SUB(calculate_vaf)
307    end subroutine calculate_vaf
308
309  end program vibrational
310
311!! Local Variables:
312!! mode: f90
313!! coding: utf-8
314!! End:
315