Octopus
propagation_spectrum.F90
Go to the documentation of this file.
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
23 use debug_oct_m
24 use global_oct_m
25 use io_oct_m
26 use kick_oct_m
29 use parser_oct_m
33
34 implicit none
35
36 integer :: in_file(3), out_file(3), ref_file, eq_axes, nspin, &
37 lmax, time_steps, ierr
38 logical :: calculate_tensor, reference_multipoles
39 type(spectrum_t) :: spectrum
40 type(unit_system_t) :: file_units
41 character(len=80) :: refmultipoles
42
43 ! Initialize stuff
45
46 call getopt_init(ierr)
47 refmultipoles = ""
48 if (ierr == 0) call getopt_propagation_spectrum(refmultipoles)
49 call getopt_end()
50
51 call parser_init()
52
53 call messages_init()
54
55 call io_init()
57
59
60 call spectrum_init(spectrum, global_namespace)
61
62 select case (spectrum%spectype)
64 call read_files('multipoles', refmultipoles)
65 call calculate_absorption('cross_section', global_namespace)
66 case (spectrum_p_power)
67 call calculate_dipole_power("multipoles", 'dipole_power')
69 call calculate_ftchd('ftchd', 'dynamic_structure_factor')
71 call calculate_rotatory_strength("angular", "rotatory_strength")
72 case default
73 write(message(1), '(a)') 'No PropagationSpectrumType defined,'
74 write(message(2), '(a)') 'cannot calculate the spectrum.'
75 call messages_warning(2)
76 end select
77
79 call io_end()
80 call messages_end()
81
82 call parser_end()
83
84 call global_end()
85
86contains
87
120 subroutine read_files(fname, reffname)
121 character(len=*), intent(in) :: fname
122 character(len=*), intent(in) :: reffname
123
124 type(kick_t) :: kick
125 real(real64) :: dt
126
127 push_sub(read_files)
128
129 in_file(1) = io_open(trim(fname), global_namespace, action='read', status='old', die=.false.)
130 if (in_file(1) == -1) in_file(1) = io_open('td.general/'//trim(fname), global_namespace, &
131 action='read', status='old', die=.false.)
132 if (in_file(1) /= -1) then
133 write(message(1),'(3a)') 'File "', trim(fname), '" found. This will be the only file to be processed.'
134 write(message(2),'(a)') '(If more than one file is to be used, the files should be called'
135 write(message(3),'(5a)') '"', trim(fname), '.1", "', trim(fname), '.2", etc.)'
136 write(message(4),'(a)')
137 call messages_info(4)
138
139 ! OK, so we only have one file. Now we have to see what it has inside.
140 call spectrum_mult_info(global_namespace, in_file(1), nspin, kick, time_steps, dt, file_units, lmax=lmax)
141 eq_axes = kick%pol_equiv_axes
142 if (eq_axes == 3) then
143 calculate_tensor = .true.
144 write(message(1),'(3a)') 'The file "', trim(fname), '" tells me that the system has three equivalent axes.'
145 write(message(2),'(a)') 'I will calculate the full tensor, written in file "XXXX_tensor".'
146 call messages_info(2)
147 else if (eq_axes == 2) then
148 write(message(1),'(3a)') 'The file "', trim(fname), '" tells me that the system has two equivalent axes.'
149 write(message(2),'(a)') 'However, I am only using this file; cannot calculate the full tensor.'
150 write(message(3),'(a)') 'A file "XXXX_vector" will be generated instead.'
151 call messages_warning(3)
152 calculate_tensor = .false.
153 else
154 write(message(1),'(3a)') 'The file "', trim(fname), '" tells me that the system has no usable symmetry. '
155 write(message(2),'(a)') 'However, I am only using this file; cannot calculate the full tensor.'
156 write(message(3),'(a)') 'A file "XXXX_vector" will be generated instead.'
157 call messages_warning(3)
158 calculate_tensor = .false.
159 end if
160
161 else ! We will try to load more fname.1 files...
162
163 ! In this case, we will always want the full tensor
164 calculate_tensor = .true.
165
166 in_file(1) = io_open(trim(fname)//'.1', global_namespace, action='read', &
167 status='old', die=.false.)
168 if (in_file(1) == -1) in_file(1) = io_open('td.general/'//trim(fname)//'.1', global_namespace, &
169 action='read', status='old', die=.false.)
170 if (in_file(1) == -1) then ! Could not find proper files. Die and complain.
171 write(message(1),'(5a)') 'No "', trim(fname), '" or "', trim(fname), '.1" file found. At least one of those'
172 write(message(2),'(a)') 'should be visible.'
173 call messages_fatal(2)
174 end if
175
176 call spectrum_mult_info(global_namespace, in_file(1), nspin, kick, time_steps, dt, file_units, lmax=lmax)
177 eq_axes = kick%pol_equiv_axes
178
179 if (eq_axes == 3) then
180 write(message(1),'(3a)') 'The file "', trim(fname), '.1" tells me that the system has three equivalent axes.'
181 write(message(2),'(a)') 'I will calculate the full tensor, written in file "cross_section_tensor".'
182 call messages_info(2)
183
184 else if (eq_axes == 2) then
185 in_file(2) = io_open(trim(fname)//'.2', global_namespace, action='read', status='old', die=.false.)
186 if (in_file(2) == -1) in_file(2) = io_open('td.general/'//trim(fname)//'.2', global_namespace, &
187 action='read', status='old', die=.false.)
188 if (in_file(2) == -1) then
189 write(message(1),'(3a)') 'The file "', trim(fname), '.1" tells me that the system has two equivalent axes,'
190 write(message(2),'(3a)') 'but I cannot find a "', trim(fname), '.2".'
191 call messages_fatal(2)
192 end if
193 write(message(1),'(5a)') 'Found two files, "', trim(fname), '.1" and "', trim(fname), '.2".'
194 write(message(2),'(a)') 'Two polarization axes are equivalent. I will generate the full tensor.'
195 call messages_info(2)
196
197 else ! No equivalent axes
198 in_file(2) = io_open(trim(fname)//'.2', global_namespace, action='read', &
199 status='old', die=.false.)
200 if (in_file(2) == -1) in_file(2) = io_open('td.general/'//trim(fname)//'.2', global_namespace, &
201 action='read', status='old', die=.false.)
202 if (in_file(2) == -1) then
203 write(message(1),'(3a)') 'The file "', trim(fname), '.1" tells me that the system has three inequivalent axes,'
204 write(message(2),'(3a)') 'but I cannot find a "', trim(fname), '.2".'
205 call messages_fatal(2)
206 end if
207 in_file(3) = io_open(trim(fname)//'.3', global_namespace, action='read', &
208 status='old', die=.false.)
209 if (in_file(3) == -1) in_file(3) = io_open('td.general/'//trim(fname)//'.3', global_namespace, &
210 action='read', status='old', die=.false.)
211 if (in_file(3) == -1) then
212 write(message(1),'(3a)') 'The file "', trim(fname), '.1" tells me that the system has three inequivalent axes,'
213 write(message(2),'(3a)') 'but I cannot find a "', trim(fname), '.3".'
214 call messages_fatal(2)
215 end if
216 write(message(1),'(7a)') 'Found three files, "', trim(fname), '.1", "', trim(fname), '.2" and "', trim(fname), '.3".'
217 write(message(2),'(a)') 'No symmetry information will be used.'
218 call messages_info(2)
219 end if
220
221 end if
222
223 if (reffname == "") then
224 reference_multipoles = .false.
225 else
226 reference_multipoles = .true.
227 ref_file = io_open(trim(reffname), global_namespace, action='read', status='old')
228 end if
229
230 pop_sub(read_files)
231 end subroutine read_files
232
233
234 !----------------------------------------------------------------------------
235 subroutine calculate_absorption(fname, namespace)
236 character(len=*), intent(in) :: fname
237 type(namespace_t), intent(in) :: namespace
238
239 integer :: ii, jj
240 character(len=150), allocatable :: filename(:)
241
242 push_sub(calculate_absorption)
243
244 if (.not. calculate_tensor) then
245
246 out_file(1) = io_open(trim(fname)//'_vector', global_namespace, action='write')
247 if (.not. reference_multipoles) then
248 call spectrum_cross_section(spectrum, namespace, in_file(1), out_file(1))
249 else
250 call spectrum_cross_section(spectrum, namespace, in_file(1), out_file(1), ref_file)
251 end if
252 call io_close(in_file(1))
253 call io_close(out_file(1))
254
255 else
256 select case (eq_axes)
257 case (0, 1)
258 jj = 3
259 case (2)
260 jj = 2
261 case (3)
262 jj = 1
263 end select
264
265 safe_allocate(filename(1:jj))
266 do ii = 1, jj
267 write(filename(ii),'(2a,i1)') trim(fname), '_vector.',ii
268 out_file(ii) = io_open(trim(filename(ii)), global_namespace, action='write')
269 if (.not. reference_multipoles) then
270 call spectrum_cross_section(spectrum, namespace, in_file(ii), out_file(ii))
271 else
272 call spectrum_cross_section(spectrum, namespace, in_file(ii), out_file(ii), ref_file)
273 end if
274 call io_close(in_file(ii))
275 call io_close(out_file(ii))
276 in_file(ii) = io_open(trim(filename(ii)), global_namespace, action='read', status='old')
277 end do
278
279 out_file(1) = io_open(trim(fname)//'_tensor', global_namespace, action='write')
280 call spectrum_cross_section_tensor(spectrum, namespace, out_file(1), in_file(1:jj))
281 do ii = 1, jj
282 call io_close(in_file(ii))
283 end do
284 call io_close(out_file(1))
285
286 end if
287
288 pop_sub(calculate_absorption)
289 end subroutine calculate_absorption
290
291 !----------------------------------------------------------------------------
292 subroutine calculate_dipole_power(fname_in, fname_out)
293 character(len=*), intent(in) :: fname_in, fname_out
294
295 push_sub(calculate_dipole_power)
296
297 in_file(1) = io_open(trim(fname_in), global_namespace, action='read', status='old', die=.false.)
298 if (in_file(1) < 0) in_file(1) = io_open('td.general/'//trim(fname_in), global_namespace, &
299 action='read', status='old', die=.false.)
300 if (in_file(1) >= 0) then
301 write(message(1),'(3a)') 'File "', trim(fname_in), '" found.'
302 write(message(2),'(a)')
303 call messages_info(2)
304 end if
305
306 out_file(1) = io_open(trim(fname_out), global_namespace, action='write')
307 call spectrum_dipole_power(spectrum, global_namespace, in_file(1), out_file(1))
308
309 call io_close(in_file(1))
310 call io_close(out_file(1))
311
313 end subroutine calculate_dipole_power
314
315 !----------------------------------------------------------------------------
316 subroutine calculate_rotatory_strength(fname_in, fname_out)
317 character(len=*), intent(in) :: fname_in, fname_out
318
320
321 in_file(1) = io_open(trim(fname_in), global_namespace, action='read', status='old', die=.false.)
322 if (in_file(1) < 0) in_file(1) = io_open('td.general/'//trim(fname_in), global_namespace, &
323 action='read', status='old', die=.false.)
324 if (in_file(1) >= 0) then
325 write(message(1),'(3a)') 'File "', trim(fname_in), '" found.'
326 write(message(2),'(a)')
327 call messages_info(2)
328 end if
329
330 out_file(1) = io_open(trim(fname_out), global_namespace, action='write')
331 call spectrum_rotatory_strength(spectrum, global_namespace, in_file(1), out_file(1))
332
333 call io_close(in_file(1))
334 call io_close(out_file(1))
335
337 end subroutine calculate_rotatory_strength
338
339
340 !----------------------------------------------------------------------------
341 subroutine calculate_ftchd(fname_in, fname_out)
342 character(len=*), intent(in) :: fname_in, fname_out
343
344 push_sub(calculate_ftchd)
345
346 ! read files
347 in_file(1) = io_open(trim(fname_in) // '.sin', global_namespace, action='read', status='old', die=.false.)
348 in_file(2) = io_open(trim(fname_in) // '.cos', global_namespace, action='read', status='old', die=.false.)
349
350 out_file(1) = io_open(trim(fname_out), global_namespace, action='write')
351 call spectrum_dyn_structure_factor(spectrum, global_namespace, in_file(1), in_file(2), out_file(1))
352 call io_close(in_file(1))
353 call io_close(out_file(1))
354
355 pop_sub(calculate_ftchd)
356 end subroutine calculate_ftchd
357
358end program propagation_spectrum
359
360!! Local Variables:
361!! mode: f90
362!! coding: utf-8
363!! End:
subroutine, public getopt_init(ierr)
Initializes the getopt machinery. Must be called before attempting to parse the options....
subroutine, public getopt_end
subroutine, public global_end()
Finalise parser varinfo file, and MPI.
Definition: global.F90:382
type(mpi_comm), parameter, public serial_dummy_comm
Alias MPI_COMM_UNDEFINED for the specific use case of initialising Octopus utilities with no MPI supp...
Definition: global.F90:265
subroutine, public init_octopus_globals(comm)
Initialise Octopus-specific global constants and files. This routine performs no initialisation calls...
Definition: global.F90:354
Definition: io.F90:114
subroutine, public io_init(defaults)
If the argument defaults is present and set to true, then the routine will not try to read anything f...
Definition: io.F90:161
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
subroutine, public io_end()
Definition: io.F90:258
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
subroutine, public messages_end()
Definition: messages.F90:277
subroutine, public messages_init(output_dir)
Definition: messages.F90:224
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
type(namespace_t), public global_namespace
Definition: namespace.F90:132
subroutine, public parser_init()
Initialise the Octopus parser.
Definition: parser.F90:450
subroutine, public parser_end()
End the Octopus parser.
Definition: parser.F90:481
subroutine, public profiling_end(namespace)
Definition: profiling.F90:413
subroutine, public profiling_init(namespace)
Create profiling subdirectory.
Definition: profiling.F90:255
subroutine, public spectrum_cross_section(spectrum, namespace, in_file, out_file, ref_file)
Definition: spectrum.F90:710
subroutine, public spectrum_cross_section_tensor(spectrum, namespace, out_file, in_file)
Definition: spectrum.F90:430
subroutine, public spectrum_dyn_structure_factor(spectrum, namespace, in_file_sin, in_file_cos, out_file)
Definition: spectrum.F90:1262
subroutine, public spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
Definition: spectrum.F90:213
subroutine, public spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax)
Definition: spectrum.F90:2339
subroutine, public spectrum_dipole_power(spectrum, namespace, in_file, out_file)
Definition: spectrum.F90:1154
integer, parameter, public spectrum_energyloss
Definition: spectrum.F90:176
integer, parameter, public spectrum_rotatory
Definition: spectrum.F90:176
subroutine, public spectrum_rotatory_strength(spectrum, namespace, in_file, out_file)
Definition: spectrum.F90:1415
integer, parameter, public spectrum_absorption
Definition: spectrum.F90:176
integer, parameter, public spectrum_p_power
Definition: spectrum.F90:176
This module defines the unit system, used for input and output.
subroutine, public unit_system_init(namespace)
subroutine calculate_rotatory_strength(fname_in, fname_out)
subroutine read_files(fname, reffname)
subroutine calculate_ftchd(fname_in, fname_out)
subroutine calculate_absorption(fname, namespace)
subroutine calculate_dipole_power(fname_in, fname_out)
program propagation_spectrum
int true(void)