Octopus
propagator_etrs.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
22 use accel_oct_m
23 use batch_oct_m
24 use debug_oct_m
29 use grid_oct_m
31 use global_oct_m
36 use ions_oct_m
39 use lda_u_oct_m
41 use math_oct_m
46 use parser_oct_m
49 use space_oct_m
52 use types_oct_m
53 use v_ks_oct_m
56 use xc_oct_m
57
58 implicit none
59
60 private
61
62 public :: &
63 td_etrs, &
64 td_etrs_sc, &
65 td_aetrs, &
67
68contains
69
70 ! ---------------------------------------------------------
72 subroutine td_etrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
73 ions_dyn, ions, mc)
74 type(v_ks_t), intent(inout) :: ks
75 type(namespace_t), intent(in) :: namespace
76 type(electron_space_t), intent(in) :: space
77 type(hamiltonian_elec_t), intent(inout) :: hm
78 type(partner_list_t), intent(in) :: ext_partners
79 type(grid_t), intent(inout) :: gr
80 type(states_elec_t), intent(inout) :: st
81 type(propagator_base_t), intent(inout) :: tr
82 real(real64), intent(in) :: time
83 real(real64), intent(in) :: dt
84 type(ion_dynamics_t), intent(inout) :: ions_dyn
85 type(ions_t), intent(inout) :: ions
86 type(multicomm_t), intent(inout) :: mc
87
88 type(xc_copied_potentials_t) :: vhxc_t1, vhxc_t2
89 type(gauge_field_t), pointer :: gfield
90
91 push_sub(td_etrs)
92
93 if (hm%theory_level /= independent_particles) then
94
95 call hm%ks_pot%store_potentials(vhxc_t1)
96
97 call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt, dt)
98
99 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
100 calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
101
102 call hm%ks_pot%store_potentials(vhxc_t2)
103 call hm%ks_pot%restore_potentials(vhxc_t1)
104 call hm%update(gr, namespace, space, ext_partners, time = time - dt)
105
106 else
107
108 call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
109
110 end if
111
112 ! propagate dt/2 with H(t)
113
114 ! first move the ions to time t
115 call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
116 ext_partners, mc, time, dt)
117
118 gfield => list_get_gauge_field(ext_partners)
119 if(associated(gfield)) then
120 call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
121 end if
122
123 call hm%ks_pot%restore_potentials(vhxc_t2)
124
125 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
126
127 ! propagate dt/2 with H(time - dt)
128 call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
129
130 pop_sub(td_etrs)
131 end subroutine td_etrs
132
133 ! ---------------------------------------------------------
135 subroutine td_etrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
136 ions_dyn, ions, mc, sctol, scsteps)
137 type(v_ks_t), intent(inout) :: ks
138 type(namespace_t), intent(in) :: namespace
139 type(electron_space_t), intent(in) :: space
140 type(hamiltonian_elec_t), intent(inout) :: hm
141 type(partner_list_t), intent(in) :: ext_partners
142 type(grid_t), intent(inout) :: gr
143 type(states_elec_t), intent(inout) :: st
144 type(propagator_base_t), intent(inout) :: tr
145 real(real64), intent(in) :: time
146 real(real64), intent(in) :: dt
147 type(ion_dynamics_t), intent(inout) :: ions_dyn
148 type(ions_t), intent(inout) :: ions
149 type(multicomm_t), intent(inout) :: mc
150 real(real64), intent(in) :: sctol
151 integer, optional, intent(out) :: scsteps
152
153 real(real64) :: diff
154 integer :: ik, ib, iter
155 class(wfs_elec_t), allocatable :: psi2(:, :)
156 ! these are hardcoded for the moment
157 integer, parameter :: niter = 10
158 type(gauge_field_t), pointer :: gfield
159 type(xc_copied_potentials_t) :: vhxc_t1, vhxc_t2
160
161 push_sub(td_etrs_sc)
162
163 assert(hm%theory_level /= independent_particles)
164
165 call hm%ks_pot%store_potentials(vhxc_t1)
166
167 call messages_new_line()
168 call messages_write(' Self-consistency iteration:')
169 call messages_info(namespace=namespace)
170
171 !Propagate the states to t+dt/2 and compute the density at t+dt
172 call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt, dt)
173
174 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
175 calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
176
177
178 call hm%ks_pot%store_potentials(vhxc_t2)
179 call hm%ks_pot%restore_potentials(vhxc_t1)
180
181 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt)
182
183 ! propagate dt/2 with H(t)
184
185 ! first move the ions to time t
186 call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
187 ext_partners, mc, time, dt)
188
189 gfield => list_get_gauge_field(ext_partners)
190 if(associated(gfield)) then
191 call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
192 end if
193
194 call hm%ks_pot%restore_potentials(vhxc_t2)
195
196 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
197
198 safe_allocate_type_array(wfs_elec_t, psi2, (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
199
200 ! store the state at half iteration
201 do ik = st%d%kpt%start, st%d%kpt%end
202 do ib = st%group%block_start, st%group%block_end
203 call st%group%psib(ib, ik)%copy_to(psi2(ib, ik), copy_data=.true.)
204 end do
205 end do
206
207 do iter = 1, niter
208
209 call hm%ks_pot%store_potentials(vhxc_t2)
210
211 call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half * dt)
212
213 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
214 time = time, calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
215 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
216
217 ! now check how much the potential changed
218 diff = hm%ks_pot%check_convergence(vhxc_t2, gr, st%rho, st%qtot)
219
220 call messages_write(' step ')
221 call messages_write(iter)
222 call messages_write(', residue = ')
223 call messages_write(abs(diff), fmt = '(1x,es9.2)')
224 call messages_info(namespace=namespace)
225
226 if (diff <= sctol) exit
227
228 if (iter /= niter) then
229 ! we are not converged, restore the states
230 do ik = st%d%kpt%start, st%d%kpt%end
231 do ib = st%group%block_start, st%group%block_end
232 call psi2(ib, ik)%copy_data_to(gr%np, st%group%psib(ib, ik))
233 end do
234 end do
235 end if
236
237 end do
238
239 if (hm%lda_u_level /= dft_u_none) then
240 call lda_u_write_u(hm%lda_u, namespace=namespace)
241 call lda_u_write_v(hm%lda_u, namespace=namespace)
242 end if
243
244 ! print an empty line
245 call messages_info(namespace=namespace)
246
247 if (present(scsteps)) scsteps = iter
248
249 do ik = st%d%kpt%start, st%d%kpt%end
250 do ib = st%group%block_start, st%group%block_end
251 call psi2(ib, ik)%end()
252 end do
253 end do
254
255 safe_deallocate_a(psi2)
256
257 pop_sub(td_etrs_sc)
258 end subroutine td_etrs_sc
259
260 ! ---------------------------------------------------------
262 subroutine td_aetrs(namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
263 type(namespace_t), intent(in) :: namespace
264 type(electron_space_t), intent(in) :: space
265 type(hamiltonian_elec_t), intent(inout) :: hm
266 type(grid_t), intent(inout) :: gr
267 type(states_elec_t), intent(inout) :: st
268 type(propagator_base_t), intent(inout) :: tr
269 real(real64), intent(in) :: time
270 real(real64), intent(in) :: dt
271 type(ion_dynamics_t), intent(inout) :: ions_dyn
272 type(ions_t), intent(inout) :: ions
273 type(partner_list_t), intent(in) :: ext_partners
274 type(multicomm_t), intent(inout) :: mc
275
276 type(gauge_field_t), pointer :: gfield
277
278 push_sub(td_aetrs)
279
280 ! propagate half of the time step with H(time - dt)
281 call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
282
283 !Get the potentials from the interpolation
284 call propagation_ops_elec_interpolate_get(hm, gr, tr%vks_old)
285
286 ! move the ions to time t
287 call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, &
288 ions, ext_partners, mc, time, dt)
289
290 !Propagate gauge field
291 gfield => list_get_gauge_field(ext_partners)
292 if(associated(gfield)) then
293 call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
294 end if
295
296 !Update Hamiltonian
297 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
298
299 !Do the time propagation for the second half of the time step
300 call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
301
302 pop_sub(td_aetrs)
303 end subroutine td_aetrs
304
305 ! ---------------------------------------------------------
307 subroutine td_caetrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
308 ions_dyn, ions, mc)
309 type(v_ks_t), intent(inout) :: ks
310 type(namespace_t), intent(in) :: namespace
311 type(electron_space_t), intent(in) :: space
312 type(hamiltonian_elec_t), intent(inout) :: hm
313 type(partner_list_t), intent(in) :: ext_partners
314 type(grid_t), intent(inout) :: gr
315 type(states_elec_t), intent(inout) :: st
316 type(propagator_base_t), intent(inout) :: tr
317 real(real64), intent(in) :: time
318 real(real64), intent(in) :: dt
319 type(ion_dynamics_t), intent(inout) :: ions_dyn
320 type(ions_t), intent(inout) :: ions
321 type(multicomm_t), intent(inout) :: mc
322
323 integer :: ik, ispin, ip, ist, ib
324 real(real64) :: vv
325 complex(real64) :: phase
326 type(density_calc_t) :: dens_calc
327 integer(int64) :: pnp, wgsize, dim2, dim3
328 type(accel_mem_t) :: phase_buff
329 type(gauge_field_t), pointer :: gfield
330 type(xc_copied_potentials_t) :: vold
331
332 push_sub(td_caetrs)
333
334 ! Get the interpolated KS potentials into vold
335 call hm%ks_pot%get_interpolated_potentials(tr%vks_old, 2, storage=vold)
336 ! And set it to the Hamiltonian
337 call hm%ks_pot%restore_potentials(vold)
338
339 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt)
340
341 call v_ks_calc_start(ks, namespace, space, hm, st, ions, ions%latt, ext_partners, &
342 time = time - dt, calc_energy = .false., calc_current = .false.)
343
344 ! propagate half of the time step with H(time - dt)
345 call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
346
347 call v_ks_calc_finish(ks, hm, namespace, space, ions%latt, st, ext_partners)
348
349 call hm%ks_pot%set_interpolated_potentials(tr%vks_old, 1)
350
351 call hm%ks_pot%perform_interpolation(tr%vks_old, (/time - dt, time - m_two*dt, time - m_three*dt/), time)
352
353 ! Replace vold by 0.5(vhxc+vold)
354 call hm%ks_pot%mix_potentials(vold, dt)
356 ! copy vold to a cl buffer
357 if (accel_is_enabled() .and. hm%apply_packed()) then
358 if (family_is_mgga_with_exc(hm%xc)) then
359 call messages_not_implemented('CAETRS propagator with accel and MGGA with energy functionals', namespace=namespace)
360 end if
361 pnp = accel_padded_size(gr%np)
362 call accel_create_buffer(phase_buff, accel_mem_read_only, type_float, pnp*st%d%nspin)
363 call vold%copy_vhxc_to_buffer(int(gr%np, int64), st%d%nspin, pnp, phase_buff)
364 end if
365
366 !Get the potentials from the interpolator
367 call propagation_ops_elec_interpolate_get(hm, gr, tr%vks_old)
368
369 ! move the ions to time t
370 call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, &
371 ions, ext_partners, mc, time, dt)
372
373 gfield => list_get_gauge_field(ext_partners)
374 if(associated(gfield)) then
375 call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
376 end if
377
378 call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
379
380 call density_calc_init(dens_calc, st, gr, st%rho)
381
382 ! propagate the other half with H(t)
383 do ik = st%d%kpt%start, st%d%kpt%end
384 ispin = st%d%get_spin_index(ik)
385
386 do ib = st%group%block_start, st%group%block_end
387 if (hm%apply_packed()) then
388 call st%group%psib(ib, ik)%do_pack()
389 if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_pack()
390 end if
391
392 call profiling_in("CAETRS_PHASE")
393 select case (st%group%psib(ib, ik)%status())
394 case (batch_not_packed)
395 do ip = 1, gr%np
396 vv = vold%vhxc(ip, ispin)
397 phase = cmplx(cos(vv), -sin(vv), real64)
398 do ist = 1, st%group%psib(ib, ik)%nst_linear
399 st%group%psib(ib, ik)%zff_linear(ip, ist) = st%group%psib(ib, ik)%zff_linear(ip, ist)*phase
400 end do
401 end do
402 case (batch_packed)
403 do ip = 1, gr%np
404 vv = vold%vhxc(ip, ispin)
405 phase = cmplx(cos(vv), -sin(vv), real64)
406 do ist = 1, st%group%psib(ib, ik)%nst_linear
407 st%group%psib(ib, ik)%zff_pack(ist, ip) = st%group%psib(ib, ik)%zff_pack(ist, ip)*phase
408 end do
409 end do
411 call accel_set_kernel_arg(kernel_phase, 0, pnp*(ispin - 1))
412 call accel_set_kernel_arg(kernel_phase, 1, phase_buff)
413 call accel_set_kernel_arg(kernel_phase, 2, st%group%psib(ib, ik)%ff_device)
414 call accel_set_kernel_arg(kernel_phase, 3, log2(st%group%psib(ib, ik)%pack_size(1)))
415
416 wgsize = accel_max_workgroup_size()/st%group%psib(ib, ik)%pack_size(1)
417 dim3 = pnp/(accel_max_size_per_dim(2)*wgsize) + 1
418 dim2 = min(accel_max_size_per_dim(2)*wgsize, pad(pnp, wgsize))
419
420 call accel_kernel_run(kernel_phase, (/st%group%psib(ib, ik)%pack_size(1), dim2, dim3/), &
421 (/st%group%psib(ib, ik)%pack_size(1), wgsize, 1_int64/))
422
423 end select
424 call profiling_out("CAETRS_PHASE")
425
426 call hm%phase%set_phase_corr(gr, st%group%psib(ib, ik))
427 if (hamiltonian_elec_inh_term(hm)) then
428 call tr%te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), m_half*dt, &
429 inh_psib = hm%inh_st%group%psib(ib, ik))
430 else
431 call tr%te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), m_half*dt)
432 end if
433 call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik))
434
435 call density_calc_accumulate(dens_calc, st%group%psib(ib, ik))
436
437 if (hm%apply_packed()) then
438 call st%group%psib(ib, ik)%do_unpack()
439 if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_unpack()
440 end if
441 end do
442 end do
443
444 call density_calc_end(dens_calc)
445
446 if (accel_is_enabled() .and. hm%apply_packed()) then
447 call accel_release_buffer(phase_buff)
448 end if
449
450 pop_sub(td_caetrs)
451 end subroutine td_caetrs
452
453end module propagator_etrs_oct_m
454
455!! Local Variables:
456!! mode: f90
457!! coding: utf-8
458!! End:
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
integer pure function, public accel_max_size_per_dim(dim)
Definition: accel.F90:2179
subroutine, public accel_release_buffer(this)
Definition: accel.F90:1250
type(accel_kernel_t), target, save, public kernel_phase
Definition: accel.F90:290
pure logical function, public accel_is_enabled()
Definition: accel.F90:402
integer, parameter, public accel_mem_read_only
Definition: accel.F90:183
integer pure function, public accel_max_workgroup_size()
Definition: accel.F90:1476
This module implements batches of mesh functions.
Definition: batch.F90:133
integer, parameter, public batch_not_packed
functions are stored in CPU memory, unpacked order
Definition: batch.F90:282
integer, parameter, public batch_device_packed
functions are stored in device memory in packed order
Definition: batch.F90:282
integer, parameter, public batch_packed
functions are stored in CPU memory, in transposed (packed) order
Definition: batch.F90:282
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public density_calc_end(this, allreduce, symmetrize)
Finalize the density calculation.
Definition: density.F90:559
subroutine, public density_calc_accumulate(this, psib, async)
Accumulate weighted orbital densities for a batch psib.
Definition: density.F90:392
subroutine, public density_calc_init(this, st, gr, density)
initialize the density calculator
Definition: density.F90:188
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_three
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:117
pure logical function, public hamiltonian_elec_inh_term(hm)
This module defines classes and functions for interaction partners.
A module to handle KS potential, without the external potential.
integer, parameter, public independent_particles
subroutine, public lda_u_write_u(this, iunit, namespace)
Definition: lda_u_io.F90:530
subroutine, public lda_u_write_v(this, iunit, namespace)
Definition: lda_u_io.F90:579
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, hm_base, phase, energy)
Definition: lda_u.F90:798
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_new_line()
Definition: messages.F90:1134
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
subroutine, public propagation_ops_elec_propagate_gauge_field(wo, gfield, dt, time, save_gf)
subroutine, public propagation_ops_elec_update_hamiltonian(namespace, space, st, mesh, hm, ext_partners, time)
subroutine, public propagation_ops_elec_exp_apply(te, namespace, st, mesh, hm, dt)
subroutine, public propagation_ops_elec_interpolate_get(hm, mesh, vks_old)
subroutine, public propagation_ops_elec_move_ions(wo, gr, hm, st, namespace, space, ions_dyn, ions, ext_partners, mc, time, dt, save_pos)
subroutine, public propagation_ops_elec_fuse_density_exp_apply(te, namespace, st, gr, hm, dt, dt2, vmagnus)
subroutine, public td_aetrs(namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
Propagator with approximate enforced time-reversal symmetry.
subroutine, public td_etrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, ions_dyn, ions, mc, sctol, scsteps)
Propagator with enforced time-reversal symmetry and self-consistency.
subroutine, public td_etrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, ions_dyn, ions, mc)
Propagator with enforced time-reversal symmetry.
subroutine, public td_caetrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, ions_dyn, ions, mc)
Propagator with approximate enforced time-reversal symmetry.
This module handles spin dimensions of the states and the k-point distribution.
type(type_t), public type_float
Definition: types.F90:133
subroutine, public v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
Definition: v_ks.F90:1106
subroutine, public v_ks_calc_start(ks, namespace, space, hm, st, ions, latt, ext_partners, time, calc_energy, calc_current, force_semilocal)
This routine starts the calculation of the Kohn-Sham potential. The routine v_ks_calc_finish must be ...
Definition: v_ks.F90:783
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:736
Definition: xc.F90:114
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
Definition: xc.F90:595
The calculator object.
Definition: density.F90:169
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:169
Stores all communicators and groups.
Definition: multicomm.F90:206
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:139
int true(void)