Octopus
propagation_ops_elec.F90
Go to the documentation of this file.
1!! Copyright (C) 2019 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 accel_oct_m
23 use batch_oct_m
24 use debug_oct_m
30 use global_oct_m
31 use grid_oct_m
36 use ions_oct_m
37 use, intrinsic :: iso_fortran_env
38 use lasers_oct_m
39 use lda_u_oct_m
41 use mesh_oct_m
49 use space_oct_m
53 use xc_oct_m
54
55 implicit none
56
57 private
58 public :: &
72
74 private
75 type(ion_state_t) :: ions_state
76 real(real64), allocatable :: vecpot(:), vecpot_vel(:)
78
79
80contains
81
82 ! ---------------------------------------------------------
83 subroutine propagation_ops_elec_update_hamiltonian(namespace, space, st, mesh, hm, ext_partners, time)
84 type(namespace_t), intent(in) :: namespace
85 class(space_t), intent(in) :: space
86 type(states_elec_t), intent(inout) :: st
87 class(mesh_t), intent(in) :: mesh
88 type(hamiltonian_elec_t), intent(inout) :: hm
89 type(partner_list_t), intent(in) :: ext_partners
90 real(real64), intent(in) :: time
91
92
94
95 call profiling_in('ELEC_UPDATE_H')
96
97 call calculate_mxll_dipole_field(hm, mesh, st)
98 call hm%update(mesh, namespace, space, ext_partners, time = time)
99 call lda_u_update_occ_matrices(hm%lda_u, namespace, mesh, st, hm%phase, hm%energy)
100
101 call profiling_out('ELEC_UPDATE_H')
102
105
106
107 ! ---------------------------------------------------------
108 subroutine propagation_ops_elec_move_ions(wo, gr, hm, st, namespace, space, ions_dyn, ions, &
109 ext_partners, mc, time, dt, save_pos)
110 class(propagation_ops_elec_t), intent(inout) :: wo
111 type(grid_t), intent(inout) :: gr
112 type(hamiltonian_elec_t), intent(inout) :: hm
113 type(states_elec_t), intent(inout) :: st
114 type(namespace_t), intent(in) :: namespace
115 type(electron_space_t), intent(in) :: space
116 type(ion_dynamics_t), intent(inout) :: ions_dyn
117 type(ions_t), intent(inout) :: ions
118 type(partner_list_t), intent(in) :: ext_partners
119 type(multicomm_t), intent(inout) :: mc
120 real(real64), intent(in) :: time
121 real(real64), intent(in) :: dt
122 logical, optional, intent(in) :: save_pos
123
124 real(real64) :: dt_ions
125
127
128 call profiling_in('ELEC_MOVE_IONS')
129
130
131 if (ions_dyn%is_active()) then
132 dt_ions = dt * ions_dyn%ionic_scale
133 if (optional_default(save_pos, .false.)) then
134 call ion_dynamics_save_state(ions_dyn, ions, wo%ions_state)
135 end if
136
137 call propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, ions_dyn, ions, &
138 mc, time, dt_ions)
139
140 call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = time)
141 end if
142
143 call profiling_out('ELEC_MOVE_IONS')
144
146 end subroutine propagation_ops_elec_move_ions
147
148 ! ---------------------------------------------------------
149 subroutine propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, ions_dyn, ions, &
150 mc, time, dt_ions)
151 type(grid_t), intent(inout) :: gr
152 type(hamiltonian_elec_t), intent(inout) :: hm
153 type(states_elec_t), intent(inout) :: st
154 type(namespace_t), intent(in) :: namespace
155 type(electron_space_t), intent(in) :: space
156 type(ion_dynamics_t), intent(inout) :: ions_dyn
157 type(ions_t), intent(inout) :: ions
158 type(multicomm_t), intent(inout) :: mc
159 real(real64), intent(in) :: time
160 real(real64), intent(in) :: dt_ions
161
163
164 call ion_dynamics_propagate(ions_dyn, ions, time, dt_ions, namespace)
165
166 ! Update lattice vectors and regenerate grid
167 if (ions_dyn%cell_relax()) then
168 call ion_dynamics_box_update(namespace, gr, space, ions%latt)
169 call electrons_lattice_vectors_update(namespace, gr, space, hm%psolver, hm%kpoints, &
170 mc, st%qtot, ions%latt)
171 end if
172
175
176 ! ---------------------------------------------------------
177 subroutine propagation_ops_elec_restore_ions(wo, ions_dyn, ions)
178 class(propagation_ops_elec_t), intent(inout) :: wo
179 type(ion_dynamics_t), intent(inout) :: ions_dyn
180 type(ions_t), intent(inout) :: ions
181
182
184
185 call profiling_in('ELEC_RESTORE_IONS')
186
187 if (ions_dyn%is_active()) then
188 call ion_dynamics_restore_state(ions_dyn, ions, wo%ions_state)
189 end if
190
191 call profiling_out('ELEC_RESTORE_IONS')
192
195
196 ! ---------------------------------------------------------
197 subroutine propagation_ops_elec_propagate_gauge_field(wo, gfield, dt, time, save_gf)
198 class(propagation_ops_elec_t), intent(inout) :: wo
199 type(gauge_field_t), intent(inout) :: gfield
200 real(real64), intent(in) :: dt
201 real(real64), intent(in) :: time
202 logical, optional, intent(in) :: save_gf
204
206
207 call profiling_in('ELEC_MOVE_GAUGE')
208
209 if (gauge_field_is_propagated(gfield)) then
210 if (optional_default(save_gf, .false.)) then
211 safe_allocate(wo%vecpot(1:gfield%space%dim))
212 safe_allocate(wo%vecpot_vel(1:gfield%space%dim))
213 call gauge_field_get_vec_pot(gfield, wo%vecpot)
214 call gauge_field_get_vec_pot_vel(gfield, wo%vecpot_vel)
215 end if
217 end if
218
219 call profiling_out('ELEC_MOVE_GAUGE')
220
223
224 ! ---------------------------------------------------------
225 subroutine propagation_ops_elec_restore_gauge_field(wo, namespace, space, hm, mesh, ext_partners)
226 class(propagation_ops_elec_t), intent(inout) :: wo
227 type(namespace_t), intent(in) :: namespace
228 class(space_t), intent(in) :: space
229 type(hamiltonian_elec_t), intent(inout) :: hm
230 class(mesh_t), intent(in) :: mesh
231 type(partner_list_t), intent(in) :: ext_partners
232
233 type(gauge_field_t), pointer :: gfield
234
236
237 call profiling_in('ELEC_RESTORE_GAUGE')
238
239 gfield => list_get_gauge_field(ext_partners)
240 if (associated(gfield)) then
241 call gauge_field_set_vec_pot(gfield, wo%vecpot)
242 call gauge_field_set_vec_pot_vel(gfield, wo%vecpot_vel)
243 safe_deallocate_a(wo%vecpot)
244 safe_deallocate_a(wo%vecpot_vel)
245 call hm%update(mesh, namespace, space, ext_partners)
246 end if
247
248 call profiling_out('ELEC_RESTORE_GAUGE')
249
252
253 ! ---------------------------------------------------------
254 subroutine propagation_ops_elec_exp_apply(te, namespace, st, mesh, hm, dt)
255 type(exponential_t), intent(inout) :: te
256 type(namespace_t), intent(in) :: namespace
257 type(states_elec_t), intent(inout) :: st
258 class(mesh_t), intent(in) :: mesh
259 type(hamiltonian_elec_t), intent(inout) :: hm
260 real(real64), intent(in) :: dt
261
262 integer :: ik, ib
263
265
266 call profiling_in('ELEC_EXP_APPLY')
267
268 do ik = st%d%kpt%start, st%d%kpt%end
269 call propagation_ops_do_pack(st, hm, st%group%block_start, ik)
270 do ib = st%group%block_start, st%group%block_end
271 if (ib + 1 <= st%group%block_end) call propagation_ops_do_pack(st, hm, ib+1, ik)
273
274 call hm%phase%set_phase_corr(mesh, st%group%psib(ib, ik), async=.true.)
275 if (hamiltonian_elec_inh_term(hm)) then
276 call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt, &
277 inh_psib = hm%inh_st%group%psib(ib, ik))
278 else
279 call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt)
280 end if
281 call hm%phase%unset_phase_corr(mesh, st%group%psib(ib, ik))
282
283 call propagation_ops_do_unpack(st, hm, ib, ik)
284 if (ib-1 >= st%group%block_start) call propagation_ops_finish_unpack(st, hm, ib-1, ik)
285 end do
286 call propagation_ops_finish_unpack(st, hm, st%group%block_end, ik)
287 end do
288
289 call accel_finish()
290
291 call profiling_out('ELEC_EXP_APPLY')
294
295 end subroutine propagation_ops_elec_exp_apply
296
297 ! ---------------------------------------------------------
298 subroutine propagation_ops_elec_fuse_density_exp_apply(te, namespace, st, gr, hm, dt, dt2, vmagnus)
299 type(exponential_t), intent(inout) :: te
300 type(namespace_t), intent(in) :: namespace
301 type(states_elec_t), intent(inout) :: st
302 type(grid_t), intent(in) :: gr
303 type(hamiltonian_elec_t), intent(inout) :: hm
304 real(real64), intent(in) :: dt
305 real(real64), optional, intent(in) :: dt2
306 real(real64), optional, intent(in) :: vmagnus(:,:,:)
307
308 integer :: ik, ib
309 type(wfs_elec_t) :: zpsib_dt
310 type(density_calc_t) :: dens_calc
311
313
314 call profiling_in('ELEC_FUSE_DENS_EXP_APPLY')
315
316 call density_calc_init(dens_calc, st, gr, st%rho)
317
318 do ik = st%d%kpt%start, st%d%kpt%end
319 call propagation_ops_do_pack(st, hm, st%group%block_start, ik)
320 do ib = st%group%block_start, st%group%block_end
321 if (ib + 1 <= st%group%block_end) call propagation_ops_do_pack(st, hm, ib+1, ik)
322 call accel_set_stream(ib)
323
324 call hm%phase%set_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
325 if (present(dt2)) then
326 call st%group%psib(ib, ik)%copy_to(zpsib_dt)
327 if (st%group%psib(ib, ik)%is_packed()) call zpsib_dt%do_pack(copy = .false.)
328
329 !propagate the state to dt/2 and dt, simultaneously, with H(time - dt)
330 if (hamiltonian_elec_inh_term(hm)) then
331 call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, &
332 deltat2 = dt2, inh_psib = hm%inh_st%group%psib(ib, ik))
333 else
334 call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, &
335 deltat2 = dt2)
336 end if
337 call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
338
339 !use the dt propagation to calculate the density
340 call density_calc_accumulate(dens_calc, zpsib_dt)
341
342 call zpsib_dt%end()
343 else
344 !propagate the state to dt with H(time - dt)
345 if (hamiltonian_elec_inh_term(hm)) then
346 call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, vmagnus=vmagnus, &
347 inh_psib = hm%inh_st%group%psib(ib, ik))
348 else
349 call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, vmagnus=vmagnus)
350 end if
351 call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.)
352
353 !use the dt propagation to calculate the density
354 call density_calc_accumulate(dens_calc, st%group%psib(ib, ik), async=.true.)
355 end if
356
357 call propagation_ops_do_unpack(st, hm, ib, ik)
358 if (ib-1 >= st%group%block_start) call propagation_ops_finish_unpack(st, hm, ib-1, ik)
359 end do
360 call propagation_ops_finish_unpack(st, hm, st%group%block_end, ik)
361 end do
362
363 call density_calc_end(dens_calc)
364
365 call accel_finish()
366
367 call profiling_out('ELEC_FUSE_DENS_EXP_APPLY')
368
371
372 ! ---------------------------------------------------------
373 subroutine propagation_ops_do_pack(st, hm, ib, ik)
374 type(states_elec_t), intent(inout) :: st
375 type(hamiltonian_elec_t), intent(inout) :: hm
376 integer, intent(in) :: ib
377 integer, intent(in) :: ik
378
380 if (hm%apply_packed()) then
381 call accel_set_stream(ib)
382 call st%group%psib(ib, ik)%do_pack(async=.true.)
383 if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_pack(async=.true.)
384 end if
386 end subroutine propagation_ops_do_pack
387
388 ! ---------------------------------------------------------
389 subroutine propagation_ops_do_unpack(st, hm, ib, ik)
390 type(states_elec_t), intent(inout) :: st
391 type(hamiltonian_elec_t), intent(inout) :: hm
392 integer, intent(in) :: ib
393 integer, intent(in) :: ik
394
396 if (hm%apply_packed()) then
397 call accel_set_stream(ib)
398 call st%group%psib(ib, ik)%do_unpack(async=.true.)
399 if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_unpack(async=.true.)
400 end if
402 end subroutine propagation_ops_do_unpack
403
404 ! ---------------------------------------------------------
405 subroutine propagation_ops_finish_unpack(st, hm, ib, ik)
406 type(states_elec_t), intent(inout) :: st
407 type(hamiltonian_elec_t), intent(inout) :: hm
408 integer, intent(in) :: ib
409 integer, intent(in) :: ik
410
412 if (hm%apply_packed()) then
413 call accel_set_stream(ib)
414 call st%group%psib(ib, ik)%finish_unpack()
415 if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%finish_unpack()
416 end if
418 end subroutine propagation_ops_finish_unpack
419
420 ! ---------------------------------------------------------
421 subroutine propagation_ops_elec_interpolate_get(hm, mesh, vks_old)
422 type(hamiltonian_elec_t), intent(inout) :: hm
423 class(mesh_t), intent(in) :: mesh
424 type(potential_interpolation_t), intent(inout) :: vks_old
425
427
428 call hm%ks_pot%get_interpolated_potentials(vks_old, 0)
429
431
433
434 ! ---------------------------------------------------------
435 subroutine calculate_mxll_dipole_field(hm, mesh, st)
436 type(hamiltonian_elec_t), target,intent(inout) :: hm
437 class(mesh_t), intent(in) :: mesh
438 type(states_elec_t), intent(in) :: st
439
440 real(real64), parameter :: density_threshold = 1.0e-8_real64
441 integer :: idir
442
443 if (.not. hm%mxll%calc_field_dip) return
444
445 select case (hm%mxll%coupling_mode)
447 if (hm%mxll%add_electric_dip) then
448 hm%mxll%e_field_dip = get_field_dip(hm%mxll%e_field)
449 end if
450 if (hm%mxll%add_magnetic_dip) then
451 hm%mxll%b_field_dip = get_field_dip(hm%mxll%b_field)
452 end if
454 hm%mxll%vec_pot_dip = get_field_dip(hm%mxll%vec_pot)
455 case default
456 return
457 end select
458
459 contains
460
461 function get_field_dip(field_full) result(field_dip)
462 real(real64) :: field_dip(3)
463 real(real64), intent(in) :: field_full(:,:)
464
465 real(real64), allocatable :: total_density(:), mask_density(:)
466 real(real64) :: integral_mask
467
468 field_dip(1:3) = m_zero
469
470 select case (hm%mxll%dipole_field)
471 case (dipole_average)
472 safe_allocate(total_density(1:mesh%np))
473 safe_allocate(mask_density(size(total_density)))
474 total_density = sum(st%rho(1:mesh%np,:), dim=2)
475 ! This mask defines a region in which the fields are averaged
476 ! The region is determined by a non-zero electronic density (density > density_threshold)
477 ! The purpose of defining a region is to average the EM fields in regions overlapping with the
478 ! electronic system, and not including values of the EM fields in the rest of the box.
479 mask_density = merge(m_one, m_zero, total_density > density_threshold)
480 integral_mask = dmf_integrate(mesh, mask_density)
481 do idir = 1, mesh%box%dim
482 field_dip(idir) = dmf_integrate(mesh, mask_density*field_full(:,idir))/integral_mask
483 end do
484 safe_deallocate_a(total_density)
485 safe_deallocate_a(mask_density)
486
487 case (dipole_at_com)
488 if (mesh%mpi_grp%rank == hm%mxll%center_of_mass_rankmin) then
489 field_dip(1:3) = field_full(hm%mxll%center_of_mass_ip,1:3)
490 end if
491 call mesh%allreduce(field_dip(:))
492 end select
493
494 end function get_field_dip
495
496 end subroutine calculate_mxll_dipole_field
497
498
501
502!! Local Variables:
503!! mode: f90
504!! coding: utf-8
505!! End:
subroutine, public accel_finish()
Definition: accel.F90:985
subroutine, public accel_set_stream(stream_number)
Definition: accel.F90:1420
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
subroutine, public density_calc_end(this, allreduce, symmetrize)
Finalize the density calculation.
Definition: density.F90:561
subroutine, public density_calc_accumulate(this, psib, async)
Accumulate weighted orbital densities for a batch psib.
Definition: density.F90:394
subroutine, public density_calc_init(this, st, gr, density)
initialize the density calculator
Definition: density.F90:190
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
subroutine, public gauge_field_set_vec_pot_vel(this, vec_pot_vel)
subroutine, public gauge_field_set_vec_pot(this, vec_pot)
subroutine, public gauge_field_get_vec_pot(this, vec_pot)
subroutine, public gauge_field_do_algorithmic_operation(this, operation, dt, time)
subroutine, public gauge_field_get_vec_pot_vel(this, vec_pot_vel)
logical pure function, public gauge_field_is_propagated(this)
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_one
Definition: global.F90:192
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
pure logical function, public hamiltonian_elec_inh_term(hm)
This module defines classes and functions for interaction partners.
subroutine, public ion_dynamics_restore_state(this, ions, state)
subroutine, public ion_dynamics_propagate(this, ions, time, dt, namespace)
Interface for the ion/cell dynamics.
subroutine, public ion_dynamics_save_state(this, ions, state)
subroutine, public electrons_lattice_vectors_update(namespace, gr, space, psolver, kpoints, mc, qtot, new_latt)
subroutine, public ion_dynamics_box_update(namespace, gr, space, new_latt)
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, phase, energy)
Definition: lda_u.F90:798
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
integer, parameter, public length_gauge_dipole
integer, parameter, public dipole_average
integer, parameter, public velocity_gauge_dipole
integer, parameter, public multipolar_expansion
integer, parameter, public dipole_at_com
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public propagation_ops_elec_restore_ions(wo, ions_dyn, ions)
subroutine, public propagation_ops_elec_propagate_gauge_field(wo, gfield, dt, time, save_gf)
subroutine, public propagation_ops_do_unpack(st, hm, ib, ik)
subroutine, public propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, ions_dyn, ions, mc, time, dt_ions)
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_do_pack(st, hm, ib, ik)
subroutine, public propagation_ops_elec_restore_gauge_field(wo, namespace, space, hm, mesh, ext_partners)
subroutine, public propagation_ops_finish_unpack(st, hm, ib, ik)
subroutine calculate_mxll_dipole_field(hm, mesh, st)
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)
type(algorithmic_operation_t), parameter, public op_verlet_compute_acc
Definition: xc.F90:116
real(real64) function, dimension(3) get_field_dip(field_full)
The calculator object.
Definition: density.F90:171
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:141
int true(void)