Skip to content

Commit 0a9394b

Browse files
glb96glb96
authored
Linear density delta kick and restart
1 parent add1d78 commit 0a9394b

37 files changed

+2327
-148
lines changed

‎exts/dbcsr‎

Submodule dbcsr updated 41 files

‎src/emd/rt_delta_pulse.F‎

Lines changed: 80 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,7 @@ MODULE rt_delta_pulse
6363
USE moments_utils, ONLY: get_reference_point
6464
USE parallel_gemm_api, ONLY: parallel_gemm
6565
USE particle_types, ONLY: particle_type
66+
USE qs_dftb_matrices, ONLY: build_dftb_overlap
6667
USE qs_environment_types, ONLY: get_qs_env,&
6768
qs_environment_type
6869
USE qs_kind_types, ONLY: qs_kind_type
@@ -72,7 +73,9 @@ MODULE rt_delta_pulse
7273
build_local_magmom_matrix,&
7374
build_local_moment_matrix
7475
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
75-
USE rt_propagation_types, ONLY: rt_prop_type
76+
USE rt_propagation_types, ONLY: get_rtp,&
77+
rt_prop_create_mos,&
78+
rt_prop_type
7679
#include "../base/base_uses.f90"
7780

7881
IMPLICIT NONE
@@ -81,12 +84,79 @@ MODULE rt_delta_pulse
8184

8285
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse'
8386

84-
PUBLIC :: apply_delta_pulse_periodic, &
85-
apply_delta_pulse, &
86-
apply_delta_pulse_mag
87+
PUBLIC :: apply_delta_pulse
8788

8889
CONTAINS
8990

91+
! **************************************************************************************************
92+
!> \brief Interface to call the delta pulse depending on the type of calculation.
93+
!> \param qs_env ...
94+
!> \param rtp ...
95+
!> \param rtp_control ...
96+
!> \author Update: Guillaume Le Breton (2023.01)
97+
! **************************************************************************************************
98+
99+
SUBROUTINE apply_delta_pulse(qs_env, rtp, rtp_control)
100+
TYPE(qs_environment_type), POINTER :: qs_env
101+
TYPE(rt_prop_type), POINTER :: rtp
102+
TYPE(rtp_control_type), POINTER :: rtp_control
103+
104+
LOGICAL :: my_apply_pulse, periodic_cell
105+
TYPE(cell_type), POINTER :: cell
106+
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
107+
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
108+
TYPE(dft_control_type), POINTER :: dft_control
109+
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
110+
111+
CALL get_qs_env(qs_env, &
112+
cell=cell, &
113+
dft_control=dft_control, &
114+
matrix_s=matrix_s)
115+
periodic_cell = ANY(cell%perd > 0)
116+
my_apply_pulse = .TRUE.
117+
CALL get_qs_env(qs_env, mos=mos)
118+
IF (rtp%linear_scaling) THEN
119+
IF (.NOT. ASSOCIATED(mos)) THEN
120+
CALL cp_warn(__LOCATION__, "Delta Pulse not implemented for Linear-Scaling based ground "// &
121+
"state calculation. If you want to perform a Linear-Scaling RTP from a "// &
122+
"Linear-Scaling GS calculation you can do the following: (i) LSCF froms "// &
123+
"scratch, (ii) MO-based SCF (for 1 SCF loop for instance) with the LSCF "// &
124+
"result as a restart and (iii) linear scaling RTP + delta kick (for 1 "// &
125+
"SCF loop for instance).")
126+
my_apply_pulse = .FALSE.
127+
ELSE
128+
! create temporary mos_old and mos_new to use delta kick routine designed for MOs-based RTP
129+
CALL rt_prop_create_mos(rtp, mos, qs_env%mpools, dft_control, &
130+
init_mos_old=.TRUE., init_mos_new=.TRUE., &
131+
init_mos_next=.FALSE., init_mos_admn=.FALSE.)
132+
END IF
133+
END IF
134+
135+
IF (my_apply_pulse) THEN
136+
CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
137+
IF (rtp_control%apply_delta_pulse) THEN
138+
IF (dft_control%qs_control%dftb) &
139+
CALL build_dftb_overlap(qs_env, 1, matrix_s)
140+
IF (rtp_control%periodic) THEN
141+
CALL apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
142+
ELSE
143+
IF (periodic_cell) THEN
144+
CPWARN("This application of the delta pulse is not compatible with PBC!")
145+
END IF
146+
CALL apply_delta_pulse_electric(qs_env, mos_old, mos_new)
147+
END IF
148+
ELSE IF (rtp_control%apply_delta_pulse_mag) THEN
149+
IF (periodic_cell) THEN
150+
CPWARN("This application of the delta pulse is not compatible with PBC!")
151+
END IF
152+
CALL apply_delta_pulse_mag(qs_env, mos_old, mos_new)
153+
ELSE
154+
CPABORT("Code error: this case should not happen!")
155+
END IF
156+
END IF
157+
158+
END SUBROUTINE apply_delta_pulse
159+
90160
! **************************************************************************************************
91161
!> \brief uses perturbation theory to get the proper initial conditions
92162
!> The len_rep option is NOT compatible with periodic boundary conditions!
@@ -96,11 +166,11 @@ MODULE rt_delta_pulse
96166
!> \author Joost & Martin (2011)
97167
! **************************************************************************************************
98168

99-
SUBROUTINE apply_delta_pulse_periodic(qs_env, mos_old, mos_new)
169+
SUBROUTINE apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
100170
TYPE(qs_environment_type), POINTER :: qs_env
101171
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new
102172

103-
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_periodic'
173+
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric_periodic'
104174

105175
INTEGER :: handle, icol, idir, irow, ispin, nao, &
106176
ncol_local, nmo, nrow_local, nvirt, &
@@ -302,7 +372,7 @@ SUBROUTINE apply_delta_pulse_periodic(qs_env, mos_old, mos_new)
302372

303373
CALL timestop(handle)
304374

305-
END SUBROUTINE apply_delta_pulse_periodic
375+
END SUBROUTINE apply_delta_pulse_electric_periodic
306376

307377
! **************************************************************************************************
308378
!> \brief applies exp(ikr) to the wavefunction.... stored in mos_old...
@@ -312,11 +382,11 @@ END SUBROUTINE apply_delta_pulse_periodic
312382
!> \author Joost & Martin (2011)
313383
! **************************************************************************************************
314384

315-
SUBROUTINE apply_delta_pulse(qs_env, mos_old, mos_new)
385+
SUBROUTINE apply_delta_pulse_electric(qs_env, mos_old, mos_new)
316386
TYPE(qs_environment_type), POINTER :: qs_env
317387
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new
318388

319-
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse'
389+
CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric'
320390

321391
INTEGER :: handle, i, nao, nmo
322392
REAL(KIND=dp), DIMENSION(3) :: kvec
@@ -382,7 +452,7 @@ SUBROUTINE apply_delta_pulse(qs_env, mos_old, mos_new)
382452

383453
CALL timestop(handle)
384454

385-
END SUBROUTINE apply_delta_pulse
455+
END SUBROUTINE apply_delta_pulse_electric
386456

387457
! **************************************************************************************************
388458
!> \brief apply magnetic delta pulse to linear order

‎src/emd/rt_propagation_utils.F‎

Lines changed: 24 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -266,6 +266,7 @@ SUBROUTINE get_restart_wfn(qs_env)
266266
unit_nr
267267
REAL(KIND=dp) :: alpha, cs_pos
268268
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
269+
TYPE(cp_fm_type) :: mos_occ
269270
TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
270271
TYPE(cp_logger_type), POINTER :: logger
271272
TYPE(cp_para_env_type), POINTER :: para_env
@@ -315,15 +316,29 @@ SUBROUTINE get_restart_wfn(qs_env)
315316
re = 2*ispin - 1
316317
im = 2*ispin
317318
CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, ncol_global=ncol)
318-
alpha = 1.0_dp
319-
IF (SIZE(mo_array) == 1) alpha = 2*alpha
320-
CALL cp_dbcsr_plus_fm_fm_t( &
321-
sparse_matrix=rho_old(re)%matrix, &
322-
matrix_v=mo_array(ispin)%mo_coeff, matrix_g=mo_array(ispin)%mo_coeff, ncol=ncol, &
323-
keep_sparsity=.FALSE., alpha=alpha)
324-
END DO
325-
DO i = 1, nspin
326-
CALL dbcsr_copy(rho_new(i)%matrix, rho_old(i)%matrix)
319+
CALL cp_fm_create(mos_occ, &
320+
matrix_struct=mo_array(ispin)%mo_coeff%matrix_struct, &
321+
name="mos_occ")
322+
CALL cp_fm_to_fm(mo_array(ispin)%mo_coeff, mos_occ)
323+
IF (mo_array(ispin)%uniform_occupation) THEN
324+
alpha = 3.0_dp - REAL(nspin, dp)
325+
CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
326+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
327+
matrix_v=mos_occ, &
328+
ncol=ncol, &
329+
alpha=alpha, keep_sparsity=.FALSE.)
330+
ELSE
331+
alpha = 1.0_dp
332+
CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
333+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
334+
matrix_v=mo_array(ispin)%mo_coeff, &
335+
matrix_g=mos_occ, &
336+
ncol=ncol, &
337+
alpha=alpha, keep_sparsity=.FALSE.)
338+
END IF
339+
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
340+
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
341+
CALL cp_fm_release(mos_occ)
327342
END DO
328343
CALL calc_update_rho_sparse(qs_env)
329344
ELSE

‎src/emd/rt_propagator_init.F‎

Lines changed: 108 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,7 @@ MODULE rt_propagator_init
5454
put_data_to_history,&
5555
s_matrices_create
5656
USE rt_propagation_types, ONLY: get_rtp,&
57+
rt_prop_release_mos,&
5758
rt_prop_type
5859
#include "../base/base_uses.f90"
5960

@@ -406,54 +407,131 @@ END SUBROUTINE backtransform_matrix
406407

407408
! **************************************************************************************************
408409
!> \brief Computes the density matrix from the mos
410+
!> Update: Initialized the density matrix from complex mos (for
411+
!> instance after delta kick)
409412
!> \param rtp ...
410413
!> \param mos ...
414+
!> \param mos_old ...
411415
!> \author Samuel Andermatt (08.15)
416+
!> \author Guillaume Le Breton (01.23)
412417
! **************************************************************************************************
413418

414-
SUBROUTINE rt_initialize_rho_from_mos(rtp, mos)
419+
SUBROUTINE rt_initialize_rho_from_mos(rtp, mos, mos_old)
415420

416421
TYPE(rt_prop_type), POINTER :: rtp
417422
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
423+
TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mos_old
418424

419-
INTEGER :: ispin, ncol, re
425+
INTEGER :: im, ispin, ncol, re
420426
REAL(KIND=dp) :: alpha
421-
TYPE(cp_fm_type) :: mos_occ
427+
TYPE(cp_fm_type) :: mos_occ, mos_occ_im
422428
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_new, rho_old
423429

424430
CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
425431

426-
DO ispin = 1, SIZE(mos)
427-
re = 2*ispin - 1
428-
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
429-
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
430-
431-
CALL cp_fm_create(mos_occ, &
432-
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
433-
name="mos_occ")
434-
CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
435-
IF (mos(ispin)%uniform_occupation) THEN
436-
alpha = 3.0_dp - REAL(SIZE(mos), dp)
437-
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
438-
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
439-
matrix_v=mos_occ, &
440-
ncol=ncol, &
441-
alpha=alpha, keep_sparsity=.FALSE.)
442-
ELSE
443-
alpha = 1.0_dp
432+
IF (PRESENT(mos_old)) THEN
433+
! Used the mos from delta kick. Initialize both real and im part
434+
DO ispin = 1, SIZE(mos_old)/2
435+
re = 2*ispin - 1; im = 2*ispin
436+
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
437+
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
438+
CALL cp_fm_create(mos_occ, &
439+
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
440+
name="mos_occ")
441+
!Real part of rho
442+
CALL cp_fm_to_fm(mos_old(re), mos_occ)
443+
IF (mos(ispin)%uniform_occupation) THEN
444+
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
445+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
446+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
447+
matrix_v=mos_occ, &
448+
ncol=ncol, &
449+
alpha=alpha, keep_sparsity=.FALSE.)
450+
ELSE
451+
alpha = 1.0_dp
452+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
453+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
454+
matrix_v=mos_old(re), &
455+
matrix_g=mos_occ, &
456+
ncol=ncol, &
457+
alpha=alpha, keep_sparsity=.FALSE.)
458+
END IF
459+
460+
! Add complex part of MOs, i*i=-1
461+
CALL cp_fm_to_fm(mos_old(im), mos_occ)
462+
IF (mos(ispin)%uniform_occupation) THEN
463+
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
464+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
465+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
466+
matrix_v=mos_occ, &
467+
ncol=ncol, &
468+
alpha=alpha, keep_sparsity=.FALSE.)
469+
ELSE
470+
alpha = 1.0_dp
471+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
472+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
473+
matrix_v=mos_old(im), &
474+
matrix_g=mos_occ, &
475+
ncol=ncol, &
476+
alpha=alpha, keep_sparsity=.FALSE.)
477+
END IF
478+
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
479+
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
480+
481+
! Imaginary part of rho
482+
CALL cp_fm_create(mos_occ_im, &
483+
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
484+
name="mos_occ_im")
485+
alpha = 3.0_dp - REAL(SIZE(mos_old)/2, dp)
486+
CALL cp_fm_to_fm(mos_old(re), mos_occ)
444487
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
445-
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
446-
matrix_v=mos(ispin)%mo_coeff, &
488+
CALL cp_fm_to_fm(mos_old(im), mos_occ_im)
489+
CALL cp_fm_column_scale(mos_occ_im, mos(ispin)%occupation_numbers/alpha)
490+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(im)%matrix, &
491+
matrix_v=mos_occ_im, &
447492
matrix_g=mos_occ, &
448493
ncol=ncol, &
449-
alpha=alpha, keep_sparsity=.FALSE.)
450-
END IF
494+
alpha=2.0_dp*alpha, &
495+
symmetry_mode=-1, keep_sparsity=.FALSE.)
451496

452-
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
453-
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
454-
CALL cp_fm_release(mos_occ)
455-
456-
END DO
497+
CALL dbcsr_filter(rho_old(im)%matrix, rtp%filter_eps)
498+
CALL dbcsr_copy(rho_new(im)%matrix, rho_old(im)%matrix)
499+
CALL cp_fm_release(mos_occ_im)
500+
CALL cp_fm_release(mos_occ)
501+
END DO
502+
! Release the mos used to apply the delta kick, no longer required
503+
CALL rt_prop_release_mos(rtp)
504+
ELSE
505+
DO ispin = 1, SIZE(mos)
506+
re = 2*ispin - 1
507+
CALL dbcsr_set(rho_old(re)%matrix, 0.0_dp)
508+
CALL cp_fm_get_info(mos(ispin)%mo_coeff, ncol_global=ncol)
509+
510+
CALL cp_fm_create(mos_occ, &
511+
matrix_struct=mos(ispin)%mo_coeff%matrix_struct, &
512+
name="mos_occ")
513+
CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos_occ)
514+
IF (mos(ispin)%uniform_occupation) THEN
515+
alpha = 3.0_dp - REAL(SIZE(mos), dp)
516+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
517+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
518+
matrix_v=mos_occ, &
519+
ncol=ncol, &
520+
alpha=alpha, keep_sparsity=.FALSE.)
521+
ELSE
522+
alpha = 1.0_dp
523+
CALL cp_fm_column_scale(mos_occ, mos(ispin)%occupation_numbers/alpha)
524+
CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
525+
matrix_v=mos(ispin)%mo_coeff, &
526+
matrix_g=mos_occ, &
527+
ncol=ncol, &
528+
alpha=alpha, keep_sparsity=.FALSE.)
529+
END IF
530+
CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
531+
CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
532+
CALL cp_fm_release(mos_occ)
533+
END DO
534+
END IF
457535

458536
END SUBROUTINE rt_initialize_rho_from_mos
459537

0 commit comments

Comments
 (0)