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