(git:a145afa)
Loading...
Searching...
No Matches
rt_projection_mo_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Function related to MO projection in RTP calculations
10!> \author Guillaume Le Breton 04.2023
11! **************************************************************************************************
12MODULE rt_projection_mo_utils
13 USE cp_control_types, ONLY: dft_control_type,&
14 proj_mo_type,&
15 rtp_control_type
16 USE cp_dbcsr_api, ONLY: dbcsr_p_type
17 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
18 USE cp_files, ONLY: close_file,&
19 open_file
20 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
21 cp_fm_trace
22 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 cp_fm_struct_release,&
24 cp_fm_struct_type
25 USE cp_fm_types, ONLY: cp_fm_create,&
26 cp_fm_release,&
27 cp_fm_to_fm,&
28 cp_fm_type
29 USE cp_log_handling, ONLY: cp_get_default_logger,&
30 cp_logger_get_default_io_unit,&
31 cp_logger_type,&
32 cp_to_string
33 USE cp_output_handling, ONLY: cp_p_file,&
34 cp_print_key_finished_output,&
35 cp_print_key_generate_filename,&
36 cp_print_key_should_output,&
37 cp_print_key_unit_nr
38 USE input_section_types, ONLY: section_vals_get,&
39 section_vals_get_subs_vals,&
40 section_vals_type,&
41 section_vals_val_get
42 USE kinds, ONLY: default_string_length,&
43 dp
44 USE message_passing, ONLY: mp_para_env_type
45 USE particle_types, ONLY: particle_type
46 USE qs_environment_types, ONLY: get_qs_env,&
47 qs_environment_type
48 USE qs_kind_types, ONLY: qs_kind_type
49 USE qs_mo_io, ONLY: read_mos_restart_low
50 USE qs_mo_types, ONLY: deallocate_mo_set,&
51 mo_set_type
52 USE rt_propagation_types, ONLY: get_rtp,&
53 rt_prop_type
54#include "./../base/base_uses.f90"
55
56 IMPLICIT NONE
57 PRIVATE
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_projection_mo_utils'
60
61 PUBLIC :: init_mo_projection, compute_and_write_proj_mo
62
63CONTAINS
64
65! **************************************************************************************************
66!> \brief Initialize the mo projection objects for time dependent run
67!> \param qs_env ...
68!> \param rtp_control ...
69!> \author Guillaume Le Breton (04.2023)
70! **************************************************************************************************
71 SUBROUTINE init_mo_projection(qs_env, rtp_control)
72 TYPE(qs_environment_type), POINTER :: qs_env
73 TYPE(rtp_control_type), POINTER :: rtp_control
74
75 INTEGER :: i_rep, j_td, n_rep_val, nbr_mo_td_max, &
76 nrep
77 INTEGER, DIMENSION(:), POINTER :: tmp_ints
78 TYPE(cp_logger_type), POINTER :: logger
79 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
80 TYPE(proj_mo_type), POINTER :: proj_mo
81 TYPE(section_vals_type), POINTER :: input, print_key, proj_mo_section
82
83 NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
84 input, proj_mo_section, print_key, mos)
85
86 CALL get_qs_env(qs_env, input=input, mos=mos)
87
88 proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
89
90 ! Read the input section and load the reference MOs
91 CALL section_vals_get(proj_mo_section, n_repetition=nrep)
92 ALLOCATE (rtp_control%proj_mo_list(nrep))
93
94 DO i_rep = 1, nrep
95 NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
96 ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
97 proj_mo => rtp_control%proj_mo_list(i_rep)%proj_mo
98
99 CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
100 c_val=proj_mo%ref_mo_file_name)
101
102 CALL section_vals_val_get(proj_mo_section, "REF_ADD_LUMO", i_rep_section=i_rep, &
103 i_val=proj_mo%ref_nlumo)
104
105 ! Relevent only in EMD
106 IF (.NOT. rtp_control%fixed_ions) THEN
107 CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, &
108 l_val=proj_mo%propagate_ref)
109 END IF
110
111 ! If no reference .wfn is provided, using the restart SCF file:
112 IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
113 CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
114 IF (n_rep_val > 0) THEN
115 CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", c_val=proj_mo%ref_mo_file_name)
116 ELSE
117 !try to read from the filename that is generated automatically from the printkey
118 print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
119 logger => cp_get_default_logger()
120 proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
121 extension=".wfn", my_local=.false.)
122 END IF
123 END IF
124
125 CALL section_vals_val_get(proj_mo_section, "REF_MO_INDEX", i_rep_section=i_rep, &
126 i_vals=tmp_ints)
127 ALLOCATE (proj_mo%ref_mo_index, source=tmp_ints(:))
128 CALL section_vals_val_get(proj_mo_section, "REF_MO_SPIN", i_rep_section=i_rep, &
129 i_val=proj_mo%ref_mo_spin)
130
131 ! Read the SCF mos and store the one required
132 CALL read_reference_mo_from_wfn(qs_env, proj_mo)
133
134 ! Initialize the other parameters related to the TD mos.
135 CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
136 l_val=proj_mo%sum_on_all_ref)
137
138 CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
139 i_val=proj_mo%td_mo_spin)
140 IF (proj_mo%td_mo_spin > SIZE(mos)) THEN
141 CALL cp_abort(__location__, &
142 "You asked to project the time dependent BETA spin while the "// &
143 "real time DFT run has only one spin defined. "// &
144 "Please set TD_MO_SPIN to 1 or use UKS.")
145 END IF
146
147 CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
148 i_vals=tmp_ints)
149
150 nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
151
152 ALLOCATE (proj_mo%td_mo_index, source=tmp_ints(:))
153 IF (proj_mo%td_mo_index(1) == -1) THEN
154 DEALLOCATE (proj_mo%td_mo_index)
155 ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
156 ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
157 DO j_td = 1, nbr_mo_td_max
158 proj_mo%td_mo_index(j_td) = j_td
159 proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
160 END DO
161 ELSE
162 ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
163 proj_mo%td_mo_occ(:) = 0.0_dp
164 DO j_td = 1, SIZE(proj_mo%td_mo_index)
165 IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) THEN
166 CALL cp_abort(__location__, &
167 "The MO number available in the Time Dependent run "// &
168 "is smaller than the MO number you have required in TD_MO_INDEX.")
169 END IF
170 proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
171 END DO
172 END IF
173
174 CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
175 l_val=proj_mo%sum_on_all_td)
176
177 END DO
178
179 END SUBROUTINE init_mo_projection
180
181! **************************************************************************************************
182!> \brief Read the MO from .wfn file and store the required MOs for TD projections
183!> \param qs_env ...
184!> \param proj_mo ...
185!> \author Guillaume Le Breton (04.2023)
186! **************************************************************************************************
187 SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo)
188 TYPE(qs_environment_type), POINTER :: qs_env
189 TYPE(proj_mo_type), POINTER :: proj_mo
190
191 INTEGER :: i_ref, ispin, mo_index, natom, &
192 nbr_mo_max, nbr_ref_mo, nspins, &
193 real_mo_index, restart_unit
194 LOGICAL :: is_file
195 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
196 TYPE(cp_fm_type) :: mo_coeff_temp
197 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
198 TYPE(dft_control_type), POINTER :: dft_control
199 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_qs, mo_ref_temp
200 TYPE(mo_set_type), POINTER :: mo_set
201 TYPE(mp_para_env_type), POINTER :: para_env
202 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
203 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
204
205 NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
206 mo_ref_fmstruct, matrix_s)
207
208 CALL get_qs_env(qs_env, &
209 qs_kind_set=qs_kind_set, &
210 particle_set=particle_set, &
211 dft_control=dft_control, &
212 matrix_s_kp=matrix_s, &
213 mos=mo_qs, &
214 para_env=para_env)
215
216 natom = SIZE(particle_set, 1)
217
218 nspins = SIZE(mo_qs)
219
220 ALLOCATE (mo_ref_temp(nspins))
221
222 DO ispin = 1, nspins
223 mo_set => mo_qs(ispin)
224 mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
225 NULLIFY (mo_ref_fmstruct)
226 CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
227 ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
228 NULLIFY (mo_ref_temp(ispin)%mo_coeff)
229 ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
230 CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
231 CALL cp_fm_struct_release(mo_ref_fmstruct)
232
233 mo_ref_temp(ispin)%nao = mo_set%nao
234 mo_ref_temp(ispin)%homo = mo_set%homo
235 mo_ref_temp(ispin)%nelectron = mo_set%nelectron
236 ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
237 ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
238 NULLIFY (mo_set)
239 END DO
240
241 IF (para_env%is_source()) THEN
242 INQUIRE (file=trim(proj_mo%ref_mo_file_name), exist=is_file)
243 IF (.NOT. is_file) THEN
244 CALL cp_abort(__location__, &
245 "Reference file not found! Name of the file CP2K looked for: "//trim(proj_mo%ref_mo_file_name))
246 END IF
247
248 CALL open_file(file_name=proj_mo%ref_mo_file_name, &
249 file_action="READ", &
250 file_form="UNFORMATTED", &
251 file_status="OLD", &
252 unit_number=restart_unit)
253 END IF
254
255 CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
256 particle_set=particle_set, natom=natom, &
257 rst_unit=restart_unit)
258
259 IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
260
261 IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) THEN
262 CALL cp_abort(__location__, &
263 "Projection on spin BETA is not possible as the reference wavefunction "// &
264 "only has one spin channel. Use a reference .wfn calculated with UKS/LSD, or set REF_MO_SPIN to 1")
265 END IF
266
267 ! Store only the mos required
268 nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
269 IF (proj_mo%ref_mo_index(1) == -1) THEN
270 DEALLOCATE (proj_mo%ref_mo_index)
271 ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
272 DO i_ref = 1, nbr_mo_max
273 proj_mo%ref_mo_index(i_ref) = i_ref
274 END DO
275 ELSE
276 DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
277 IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) THEN
278 CALL cp_abort(__location__, &
279 "The number of MOs available in the reference wavefunction "// &
280 "is smaller than the MO number you have requested in REF_MO_INDEX.")
281 END IF
282 END DO
283 END IF
284 nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
285
286 IF (nbr_ref_mo > nbr_mo_max) THEN
287 CALL cp_abort(__location__, &
288 "The total number of requested MOs is larger than what is available in the reference wavefunction. "// &
289 "If you are trying to project onto virtual states, make sure they are included in the .wfn file "// &
290 "e.g., by the ADDED_MOS keyword in the SCF section of the input when calculating your reference.")
291 END IF
292
293 ! Store
294 ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
295 CALL cp_fm_struct_create(mo_ref_fmstruct, &
296 context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
297 nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
298 ncol_global=1)
299
300 IF (dft_control%rtp_control%fixed_ions) THEN
301 CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
302 END IF
303
304 DO mo_index = 1, nbr_ref_mo
305 real_mo_index = proj_mo%ref_mo_index(mo_index)
306 IF (real_mo_index > nbr_mo_max) THEN
307 CALL cp_abort(__location__, &
308 "One of reference mo index is larger then the total number of available mo in the .wfn file.")
309 END IF
310
311 ! fill with the reference mo values
312 CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
313 IF (dft_control%rtp_control%fixed_ions) THEN
314 ! multiply with overlap matrix to save time later on: proj_mo%mo_ref is SxMO_ref
315 CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
316 ncol=1, &
317 source_start=real_mo_index, &
318 target_start=1)
319 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
320 ELSE
321 ! the AO will change with times: proj_mo%mo_ref are really the MOs coeffs
322 CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, proj_mo%mo_ref(mo_index), &
323 ncol=1, &
324 source_start=real_mo_index, &
325 target_start=1)
326 END IF
327 END DO
328
329 ! Clean temporary variables
330 DO ispin = 1, nspins
331 CALL deallocate_mo_set(mo_ref_temp(ispin))
332 END DO
333 DEALLOCATE (mo_ref_temp)
334
335 CALL cp_fm_struct_release(mo_ref_fmstruct)
336 IF (dft_control%rtp_control%fixed_ions) THEN
337 CALL cp_fm_release(mo_coeff_temp)
338 END IF
339
340 END SUBROUTINE read_reference_mo_from_wfn
341
342! **************************************************************************************************
343!> \brief Compute the projection of the current MO coefficients on reference ones
344!> and write the results.
345!> \param qs_env ...
346!> \param mos_new ...
347!> \param proj_mo ...
348!> \param n_proj ...
349!> \author Guillaume Le Breton
350! **************************************************************************************************
351 SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
352 TYPE(qs_environment_type), POINTER :: qs_env
353 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
354 TYPE(proj_mo_type) :: proj_mo
355 INTEGER :: n_proj
356
357 INTEGER :: i_ref, nbr_ref_mo, nbr_ref_td
358 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase, popu, sum_popu_ref
359 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
360 TYPE(cp_fm_type) :: s_mo_ref
361 TYPE(cp_logger_type), POINTER :: logger
362 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
363 TYPE(dft_control_type), POINTER :: dft_control
364 TYPE(section_vals_type), POINTER :: input, print_mo_section, proj_mo_section
365
366 NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
367
368 logger => cp_get_default_logger()
369
370 CALL get_qs_env(qs_env, &
371 dft_control=dft_control, &
372 input=input)
373
374 ! The general section
375 proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
376 ! The section we are dealing in this particular subroutine call: n_proj.
377 print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
378
379 ! Propagate the reference MO if required at each time step
380 IF (proj_mo%propagate_ref) CALL propagate_ref_mo(qs_env, proj_mo)
381
382 ! Does not compute the projection if not the required time step
383 IF (.NOT. btest(cp_print_key_should_output(logger%iter_info, &
384 print_mo_section, ""), &
385 cp_p_file)) THEN
386 RETURN
387 END IF
388
389 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
390 CALL get_qs_env(qs_env, &
391 matrix_s_kp=matrix_s)
392 CALL cp_fm_struct_create(mo_ref_fmstruct, &
393 context=proj_mo%mo_ref(1)%matrix_struct%context, &
394 nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
395 ncol_global=1)
396 CALL cp_fm_create(s_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
397 END IF
398
399 nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
400 nbr_ref_td = SIZE(proj_mo%td_mo_index)
401 ALLOCATE (popu(nbr_ref_td))
402 ALLOCATE (phase(nbr_ref_td))
403
404 IF (proj_mo%sum_on_all_ref) THEN
405 ALLOCATE (sum_popu_ref(nbr_ref_td))
406 sum_popu_ref(:) = 0.0_dp
407 DO i_ref = 1, nbr_ref_mo
408 ! Compute SxMO_ref for the upcoming projection later on
409 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
410 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
411 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
412 ELSE
413 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
414 END IF
415 sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
416 END DO
417 IF (proj_mo%sum_on_all_td) THEN
418 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=sum(sum_popu_ref), n_proj=n_proj)
419 ELSE
420 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
421 END IF
422 DEALLOCATE (sum_popu_ref)
423 ELSE
424 DO i_ref = 1, nbr_ref_mo
425 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
426 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
427 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
428 ELSE
429 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
430 END IF
431 IF (proj_mo%sum_on_all_td) THEN
432 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=sum(popu), n_proj=n_proj)
433 ELSE
434
435 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
436 END IF
437 END DO
438 END IF
439
440 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
441 CALL cp_fm_struct_release(mo_ref_fmstruct)
442 CALL cp_fm_release(s_mo_ref)
443 END IF
444 DEALLOCATE (popu)
445 DEALLOCATE (phase)
446
447 END SUBROUTINE compute_and_write_proj_mo
448
449! **************************************************************************************************
450!> \brief Compute the projection of the current MO coefficients on reference ones
451!> \param popu ...
452!> \param phase ...
453!> \param mos_new ...
454!> \param proj_mo ...
455!> \param i_ref ...
456!> \param S_mo_ref ...
457!> \author Guillaume Le Breton
458! **************************************************************************************************
459 SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref)
460 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: popu, phase
461 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
462 TYPE(proj_mo_type) :: proj_mo
463 INTEGER :: i_ref
464 TYPE(cp_fm_type), OPTIONAL :: s_mo_ref
465
466 CHARACTER(len=*), PARAMETER :: routinen = 'compute_proj_mo'
467
468 INTEGER :: handle, j_td, nbr_ref_td, spin_td
469 LOGICAL :: is_emd
470 REAL(kind=dp) :: imag_proj, real_proj
471 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
472 TYPE(cp_fm_type) :: mo_coeff_temp
473
474 CALL timeset(routinen, handle)
475
476 is_emd = .false.
477 IF (PRESENT(s_mo_ref)) is_emd = .true.
478
479 nbr_ref_td = SIZE(popu)
480 spin_td = proj_mo%td_mo_spin
481
482 CALL cp_fm_struct_create(mo_ref_fmstruct, &
483 context=mos_new(1)%matrix_struct%context, &
484 nrow_global=mos_new(1)%matrix_struct%nrow_global, &
485 ncol_global=1)
486 CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
487
488 DO j_td = 1, nbr_ref_td
489 ! Real part of the projection:
490 real_proj = 0.0_dp
491 CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
492 ncol=1, &
493 source_start=proj_mo%td_mo_index(j_td), &
494 target_start=1)
495 IF (is_emd) THEN
496 ! The reference MO have to be propagated in the new basis, so the projection
497 CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, real_proj)
498 ELSE
499 ! The reference MO is time independent. proj_mo%mo_ref(i_ref) is in fact SxMO_ref already
500 CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
501 END IF
502
503 ! Imaginary part of the projection
504 imag_proj = 0.0_dp
505 CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
506 ncol=1, &
507 source_start=proj_mo%td_mo_index(j_td), &
508 target_start=1)
509
510 IF (is_emd) THEN
511 CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, imag_proj)
512 ELSE
513 CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
514 END IF
515
516 ! Store the result
517 phase(j_td) = atan2(imag_proj, real_proj) ! in radians
518 popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
519 END DO
520
521 CALL cp_fm_struct_release(mo_ref_fmstruct)
522 CALL cp_fm_release(mo_coeff_temp)
523
524 CALL timestop(handle)
525
526 END SUBROUTINE compute_proj_mo
527
528! **************************************************************************************************
529!> \brief Write in one file the projection of (all) the time-dependent MO coefficients
530!> onto reference ones
531!> \param qs_env ...
532!> \param print_mo_section ...
533!> \param proj_mo ...
534!> \param i_ref ...
535!> \param popu ...
536!> \param phase ...
537!> \param popu_tot ...
538!> \param n_proj ...
539!> \author Guillaume Le Breton
540! **************************************************************************************************
541 SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
542 TYPE(qs_environment_type), POINTER :: qs_env
543 TYPE(section_vals_type), POINTER :: print_mo_section
544 TYPE(proj_mo_type) :: proj_mo
545 INTEGER, OPTIONAL :: i_ref
546 REAL(kind=dp), DIMENSION(:), OPTIONAL :: popu, phase
547 REAL(kind=dp), OPTIONAL :: popu_tot
548 INTEGER, OPTIONAL :: n_proj
549
550 CHARACTER(LEN=default_string_length) :: ext, filename
551 INTEGER :: j_td, output_unit, print_unit
552 TYPE(cp_logger_type), POINTER :: logger
553
554 NULLIFY (logger)
555
556 logger => cp_get_default_logger()
557 output_unit = cp_logger_get_default_io_unit(logger)
558
559 IF (.NOT. (output_unit > 0)) RETURN
560
561 IF (proj_mo%sum_on_all_ref) THEN
562 ext = "-"//trim(adjustl(cp_to_string(n_proj)))//"-ALL_REF.dat"
563 ELSE
564 ! Filename is updated wrt the reference MO number
565 ext = "-"//trim(adjustl(cp_to_string(n_proj)))// &
566 "-REF-"// &
567 trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
568 ".dat"
569 END IF
570
571 print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
572 extension=trim(ext))
573
574 IF (print_unit /= output_unit) THEN
575 INQUIRE (unit=print_unit, name=filename)
576 WRITE (unit=print_unit, fmt="(/,(T2,A,T40,I6))") &
577 "Real time propagation step:", qs_env%sim_step
578 ELSE
579 WRITE (unit=output_unit, fmt="(/,T2,A)") "PROJECTION MO"
580 END IF
581
582 IF (proj_mo%sum_on_all_ref) THEN
583 WRITE (print_unit, "(T3,A)") &
584 "Projection on all the required MO number from the reference file "// &
585 trim(proj_mo%ref_mo_file_name)
586 IF (proj_mo%sum_on_all_td) THEN
587 WRITE (print_unit, "(T3, A, E20.12)") &
588 "The sum over all the TD MOs population:", popu_tot
589 ELSE
590 WRITE (print_unit, "(T3,A)") &
591 "For each TD MOs required is printed: Population "
592 DO j_td = 1, SIZE(popu)
593 WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
594 END DO
595 END IF
596 ELSE
597 WRITE (print_unit, "(T3,A)") &
598 "Projection on the MO number "// &
599 trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
600 " from the reference file "// &
601 trim(proj_mo%ref_mo_file_name)
602
603 IF (proj_mo%sum_on_all_td) THEN
604 WRITE (print_unit, "(T3, A, E20.12)") &
605 "The sum over all the TD MOs population:", popu_tot
606 ELSE
607 WRITE (print_unit, "(T3,A)") &
608 "For each TD MOs required is printed: Population & Phase [rad] "
609 DO j_td = 1, SIZE(popu)
610 WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
611 END DO
612 END IF
613 END IF
614
615 CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
616
617 END SUBROUTINE write_proj_mo
618
619! **************************************************************************************************
620!> \brief Propagate the reference MO in case of EMD: since the nuclei moves, the MO coeff can be
621!> propagated to represent the same MO (because the AO move with the nuclei).
622!> To do so, we use the same formula as for the electrons of the system, but without the
623!> Hamiltonian:
624!> dc^j_alpha/dt = - sum_{beta, gamma} S^{-1}_{alpha, beta} B_{beta,gamma} c^j_gamma
625!> \param qs_env ...
626!> \param proj_mo ...
627!> \author Guillaume Le Breton
628! **************************************************************************************************
629 SUBROUTINE propagate_ref_mo(qs_env, proj_mo)
630 TYPE(qs_environment_type), POINTER :: qs_env
631 TYPE(proj_mo_type) :: proj_mo
632
633 INTEGER :: i_ref
634 REAL(kind=dp) :: dt
635 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
636 TYPE(cp_fm_type) :: d_mo
637 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sinvb
638 TYPE(rt_prop_type), POINTER :: rtp
639
640 CALL get_qs_env(qs_env, rtp=rtp)
641 CALL get_rtp(rtp=rtp, sinvb=sinvb, dt=dt)
642
643 CALL cp_fm_struct_create(mo_ref_fmstruct, &
644 context=proj_mo%mo_ref(1)%matrix_struct%context, &
645 nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
646 ncol_global=1)
647 CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
648
649 DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
650 ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
651 CALL cp_dbcsr_sm_fm_multiply(sinvb(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
652 CALL cp_fm_scale_and_add(1.0_dp, proj_mo%mo_ref(i_ref), 1.0_dp, d_mo)
653 END DO
654
655 CALL cp_fm_struct_release(mo_ref_fmstruct)
656 CALL cp_fm_release(d_mo)
657
658 END SUBROUTINE propagate_ref_mo
659
660END MODULE rt_projection_mo_utils
661
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:78
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::read_mos_restart_low
subroutine, public read_mos_restart_low(mos, para_env, qs_kind_set, particle_set, natom, rst_unit, multiplicity, rt_mos, natom_mismatch)
Reading the mos from apreviously defined restart file.
Definition qs_mo_io.F:695
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:356
rt_projection_mo_utils
Function related to MO projection in RTP calculations.
Definition rt_projection_mo_utils.F:12
rt_projection_mo_utils::init_mo_projection
subroutine, public init_mo_projection(qs_env, rtp_control)
Initialize the mo projection objects for time dependent run.
Definition rt_projection_mo_utils.F:72
rt_projection_mo_utils::compute_and_write_proj_mo
subroutine, public compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
Compute the projection of the current MO coefficients on reference ones and write the results.
Definition rt_projection_mo_utils.F:352
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_control_types::proj_mo_type
Definition cp_control_types.F:72
cp_control_types::rtp_control_type
Definition cp_control_types.F:94
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:185
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76