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