Octopus
spin_susceptibility.F90
Go to the documentation of this file.
1!! Copyright (C) 2018 N. Tancogne-Dejean
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
22 use batch_oct_m
24 use global_oct_m
25 use io_oct_m
26 use ions_oct_m
27 use kick_oct_m
29 use math_oct_m
32 use parser_oct_m
34 use space_oct_m
36 use unit_oct_m
38
39 implicit none
40
41 integer :: in_file, out_file, ref_file, ii, jj, kk, ierr, ib, num_col, num_col_cart
42 integer :: time_steps, time_steps_ref, energy_steps, istart, iend, ntiter, iq
43 real(real64) :: dt, tt, ww, dt_ref
44 real(real64), allocatable :: ftreal(:,:), ftimag(:,:)
45 real(real64), allocatable :: m_cart(:,:), m_cart_ref(:,:), magnetization(:,:,:)
46 type(spectrum_t) :: spectrum
47 type(batch_t) :: vecpotb, ftrealb, ftimagb
48 type(kick_t) :: kick
49 character(len=256) :: header
50 character(len=MAX_PATH_LEN) :: fname, ref_filename
51 complex(real64), allocatable :: chi(:,:)
52
53 ! Initialize stuff
55
56 call getopt_init(ierr)
57 if (ierr == 0) call getopt_dielectric_function()
58 call getopt_end()
59
60 call parser_init()
61
62 call messages_init()
63
64 call io_init()
66
68
69 call spectrum_init(spectrum, global_namespace)
70
71 in_file = io_open('td.general/total_magnetization', global_namespace, action='read', status='old')
72 rewind(in_file)
73 read(in_file,*)
74 read(in_file,*)
75 call kick_read(kick, in_file, global_namespace)
76 call io_skip_header(in_file)
77 call spectrum_count_time_steps(global_namespace, in_file, time_steps, dt)
78
79 time_steps = time_steps + 1
80
81 num_col_cart = 12*kick%nqvec
82 safe_allocate(m_cart(1:time_steps, 1:num_col_cart))
83
84 call io_skip_header(in_file)
85
86 do ii = 1, time_steps
87 read(in_file, *) jj, tt, (m_cart(ii,ib),ib = 1, num_col_cart)
88 end do
89
90 call io_close(in_file)
91
92 if (parse_is_defined(global_namespace, 'TransientMagnetizationReference')) then
93 !%Variable TransientMagnetizationReference
94 !%Type string
95 !%Default "."
96 !%Section Utilities::oct-spin_susceptibility
97 !%Description
98 !% In case of delayed kick, the calculation of the transient spin susceptibility requires
99 !% to substract a reference calculation, containing dynamics of the magnetization without the kick
100 !% This reference must be computed having
101 !% TDOutput = total_magnetization.
102 !% This variables defined the directory in which the reference total_magnetization file is,
103 !% relative to the current folder
104 !%End
105
106 call parse_variable(global_namespace, 'TransientMagnetizationReference', '.', ref_filename)
107 ref_file = io_open(trim(ref_filename)//'/total_magnetization', global_namespace, action='read', status='old')
108 call io_skip_header(ref_file)
109 call spectrum_count_time_steps(global_namespace, ref_file, time_steps_ref, dt_ref)
110 time_steps_ref = time_steps_ref + 1
111
112 if (time_steps_ref < time_steps) then
113 message(1) = "The reference calculation does not contain enought time steps"
115 end if
116
117 if (.not. is_close(dt_ref, dt)) then
118 message(1) = "The time step of the reference calculation is different from the current calculation"
119 call messages_fatal(1)
120 end if
121
122 !We remove the reference
123 safe_allocate(m_cart_ref(1:time_steps_ref, 1:num_col_cart))
124 call io_skip_header(ref_file)
125 do ii = 1, time_steps_ref
126 read(ref_file, *) jj, tt, (m_cart_ref(ii, kk), kk = 1, num_col_cart)
127 end do
128 call io_close(ref_file)
129 do ii = 1, time_steps
130 do kk = 1, num_col_cart
131 m_cart(ii, kk) = m_cart(ii, kk) - m_cart_ref(ii, kk)
132 end do
133 end do
134 safe_deallocate_a(m_cart_ref)
135 end if
136
137 !We now perform the change of basis to the rotating basis
138 !In this basis we have only m_+(q), m_-(q), and m_z(+/-q)
139 !where z means here along the easy axis
140 num_col = 8
141 safe_allocate(magnetization(1:time_steps, 1:num_col, 1:kick%nqvec))
142
143 do iq = 1, kick%nqvec
144 !Real part of m_x
145 magnetization(:,1,iq) = m_cart(:,(iq-1)*12+1)*kick%trans_vec(1,1) &
146 + m_cart(:,(iq-1)*12+3)*kick%trans_vec(2,1) &
147 + m_cart(:,(iq-1)*12+5)*kick%trans_vec(3,1)
148 !We add -Im(m_y)
149 magnetization(:,1,iq) = magnetization(:,1,iq) -(m_cart(:,(iq-1)*12+2)*kick%trans_vec(1,2) &
150 + m_cart(:,(iq-1)*12+4)*kick%trans_vec(2,2) &
151 + m_cart(:,(iq-1)*12+6)*kick%trans_vec(3,2))
152
153 !Im part of m_x
154 magnetization(:,2,iq) = m_cart(:,(iq-1)*12+2)*kick%trans_vec(1,1) &
155 + m_cart(:,(iq-1)*12+4)*kick%trans_vec(2,1) &
156 + m_cart(:,(iq-1)*12+6)*kick%trans_vec(3,1)
157 !We add +Re(m_y)
158 magnetization(:,2,iq) = magnetization(:,2,iq) +(m_cart(:,(iq-1)*12+1)*kick%trans_vec(1,2) &
159 + m_cart(:,(iq-1)*12+3)*kick%trans_vec(2,2) &
160 + m_cart(:,(iq-1)*12+5)*kick%trans_vec(3,2))
161
162 !Real part of m_x
163 magnetization(:,3,iq) = m_cart(:,(iq-1)*12+7)*kick%trans_vec(1,1) &
164 + m_cart(:,(iq-1)*12+9)*kick%trans_vec(2,1) &
165 + m_cart(:,(iq-1)*12+11)*kick%trans_vec(3,1)
166 !We add +Im(m_y)
167 magnetization(:,3,iq) = magnetization(:,3,iq) +(m_cart(:,(iq-1)*12+8)*kick%trans_vec(1,2) &
168 + m_cart(:,(iq-1)*12+10)*kick%trans_vec(2,2) &
169 + m_cart(:,(iq-1)*12+12)*kick%trans_vec(3,2))
170
171 !Im part of m_x
172 magnetization(:,4,iq) = m_cart(:,(iq-1)*12+8)*kick%trans_vec(1,1) &
173 + m_cart(:,(iq-1)*12+10)*kick%trans_vec(2,1) &
174 + m_cart(:,(iq-1)*12+12)*kick%trans_vec(3,1)
175 !We add -Re(m_y)
176 magnetization(:,4,iq) = magnetization(:,4,iq) -(m_cart(:,(iq-1)*12+7)*kick%trans_vec(1,2) &
177 + m_cart(:,(iq-1)*12+9)*kick%trans_vec(2,2) &
178 + m_cart(:,(iq-1)*12+11)*kick%trans_vec(3,2))
179
180 !Real and Im part of m_z
181 magnetization(:,5,iq) = m_cart(:,(iq-1)*12+1)*kick%easy_axis(1) &
182 + m_cart(:,(iq-1)*12+3)*kick%easy_axis(2) &
183 + m_cart(:,(iq-1)*12+5)*kick%easy_axis(3)
184 magnetization(:,6,iq) = m_cart(:,(iq-1)*12+2)*kick%easy_axis(1) &
185 + m_cart(:,(iq-1)*12+4)*kick%easy_axis(2) &
186 + m_cart(:,(iq-1)*12+6)*kick%easy_axis(3)
187 magnetization(:,7,iq) = m_cart(:,(iq-1)*12+7)*kick%easy_axis(1) &
188 + m_cart(:,(iq-1)*12+9)*kick%easy_axis(2) &
189 + m_cart(:,(iq-1)*12+11)*kick%easy_axis(3)
190 magnetization(:,8,iq) = m_cart(:,(iq-1)*12+8)*kick%easy_axis(1) &
191 + m_cart(:,(iq-1)*12+10)*kick%easy_axis(2) &
192 + m_cart(:,(iq-1)*12+12)*kick%easy_axis(3)
193 end do
194
195 safe_deallocate_a(m_cart)
196
197 write(message(1), '(a, i7, a)') "Info: Read ", time_steps, " steps from file '"// &
198 trim(io_workpath('td.general/total_magnetization', global_namespace))//"'"
199 call messages_info(1)
200
201 write(header, '(9a15)') '# time', 'Re[m_+(q,t)]', 'Im[m_+(q,t)]', &
202 'Re[m_-(-q, t)]', 'Im[m_-(-q,t)]', &
203 'Re[m_z(q,t)]', 'Im[m_z(q,t)]', 'Re[m_z(-q,t)]', 'Im[m_z(-q,t)]'
204
205 do iq = 1, kick%nqvec
206
207 write(fname, '(a,i3.3)') 'td.general/transverse_magnetization_q', iq
208 out_file = io_open(trim(fname), global_namespace, action='write')
209 write(out_file,'(a)') trim(header)
210 do kk = 1, time_steps
211 write(out_file, '(9e15.6)') (kk - 1)*dt, (magnetization(kk,ii,iq), ii = 1,8)
212 end do
213 call io_close(out_file)
214
215
216 ! Find out the iteration numbers corresponding to the time limits.
217 call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
218
219 istart = max(1, istart)
220
221 energy_steps = spectrum_nenergy_steps(spectrum)
222
223 safe_allocate(ftreal(1:energy_steps, 1:num_col))
224 safe_allocate(ftimag(1:energy_steps, 1:num_col))
225
226 call batch_init(vecpotb, 1, 1, num_col, magnetization(:, :, iq))
227 call batch_init(ftrealb, 1, 1, num_col, ftreal)
228 call batch_init(ftimagb, 1, 1, num_col, ftimag)
229
230
231 call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart, iend, spectrum%start_time, dt, vecpotb)
232
233 call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
234 istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
235 spectrum%max_energy, spectrum%energy_step, ftrealb)
236
237 call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
238 istart, iend, spectrum%start_time, dt, vecpotb, spectrum%min_energy, &
239 spectrum%max_energy, spectrum%energy_step, ftimagb)
240
241
242 call vecpotb%end()
243 call ftrealb%end()
244 call ftimagb%end()
245
246 assert(abs(anint(num_col*m_half) - num_col*m_half) <= m_epsilon)
247 safe_allocate(chi(1:energy_steps, 1:num_col/2))
248 do ii = 1, num_col/2
249 do kk = 1, energy_steps
250 chi(kk,ii) = (ftreal(kk,(ii-1)*2+1) + m_zi*ftimag(kk, (ii-1)*2+1)&
251 -ftimag(kk, (ii-1)*2+2) + m_zi*ftreal(kk, (ii-1)*2+2))/kick%delta_strength
252 end do
253 end do
254
255 safe_deallocate_a(ftreal)
256 safe_deallocate_a(ftimag)
257
258
259 write(header, '(9a18)') '# energy', 'Re[\chi_{+-}(q)]', 'Im[\chi_{+-}(q)]', &
260 'Re[\chi_{-+}(-q)]', 'Im[\chi_{-+}(-q)]', &
261 'Re[\chi_{zz}(q)]', 'Im[\chi_{zz}(q)]', 'Re[\chi_{zz}(-q)]', 'Im[\chi_{zz}(-q)]'
262
263 write(fname, '(a,i3.3)') 'td.general/spin_susceptibility_q', iq
264 out_file = io_open(trim(fname), global_namespace, action='write')
265 write(out_file,'(a)') trim(header)
266 do kk = 1, energy_steps
267 ww = (kk-1)*spectrum%energy_step + spectrum%min_energy
268 write(out_file, '(13e15.6)') ww, &
269 (real(chi(kk,ii), real64), aimag(chi(kk,ii)), ii = 1, num_col/2)
270 end do
271 call io_close(out_file)
272
273 safe_deallocate_a(chi)
274 end do !iq
275
276 safe_deallocate_a(magnetization)
277
278 call kick_end(kick)
279
281 call io_end()
282 call messages_end()
283 call parser_end()
284 call global_end()
285
286end program spin_susceptibility
287
288!! Local Variables:
289!! mode: f90
290!! coding: utf-8
291!! End:
initialize a batch with existing memory
Definition: batch.F90:267
This module implements batches of mesh functions.
Definition: batch.F90:133
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
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
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_skip_header(iunit)
Definition: io.F90:597
subroutine, public io_end()
Definition: io.F90:258
character(len=max_path_len) function, public io_workpath(path, namespace)
Definition: io.F90:270
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
subroutine, public kick_read(kick, iunit, namespace)
Definition: kick.F90:819
subroutine, public kick_end(kick)
Definition: kick.F90:795
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public messages_end()
Definition: messages.F90:277
subroutine, public messages_init(output_dir)
Definition: messages.F90:224
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
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
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_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
Definition: spectrum.F90:2507
subroutine, public spectrum_fourier_transform(method, transform, noise, time_start, time_end, t0, time_step, time_function, energy_start, energy_end, energy_step, energy_function)
Computes the sine, cosine, (or "exponential") Fourier transform of the real function given in the tim...
Definition: spectrum.F90:2637
subroutine, public spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
Definition: spectrum.F90:213
subroutine, public spectrum_signal_damp(damp_type, damp_factor, time_start, time_end, t0, time_step, time_function)
Definition: spectrum.F90:2551
integer, parameter, public spectrum_transform_cos
Definition: spectrum.F90:171
integer, parameter, public spectrum_transform_sin
Definition: spectrum.F90:171
subroutine, public spectrum_count_time_steps(namespace, iunit, time_steps, dt)
Definition: spectrum.F90:2385
pure integer function, public spectrum_nenergy_steps(spectrum)
Definition: spectrum.F90:2946
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
This module defines the unit system, used for input and output.
subroutine, public unit_system_init(namespace)
program spin_susceptibility