(git:6a2e663)
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-2024 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 ! **************************************************************************************************
12 MODULE 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_operations, ONLY: cp_dbcsr_sm_fm_multiply
17  USE cp_files, ONLY: close_file,&
18  open_file
19  USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
20  cp_fm_trace
21  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22  cp_fm_struct_release,&
23  cp_fm_struct_type
24  USE cp_fm_types, ONLY: cp_fm_create,&
25  cp_fm_release,&
26  cp_fm_to_fm,&
27  cp_fm_type
28  USE cp_log_handling, ONLY: cp_get_default_logger,&
29  cp_logger_get_default_io_unit,&
30  cp_logger_type,&
31  cp_to_string
32  USE cp_output_handling, ONLY: cp_p_file,&
33  cp_print_key_finished_output,&
34  cp_print_key_generate_filename,&
35  cp_print_key_should_output,&
36  cp_print_key_unit_nr
37  USE dbcsr_api, ONLY: dbcsr_p_type
38  USE input_constants, ONLY: proj_mo_ref_scf,&
39  proj_mo_ref_xas_tdp
40  USE input_section_types, ONLY: section_vals_get,&
41  section_vals_get_subs_vals,&
42  section_vals_type,&
43  section_vals_val_get
44  USE kinds, ONLY: default_string_length,&
45  dp
46  USE message_passing, ONLY: mp_para_env_type
47  USE particle_types, ONLY: particle_type
48  USE qs_environment_types, ONLY: get_qs_env,&
49  qs_environment_type
50  USE qs_kind_types, ONLY: qs_kind_type
51  USE qs_mo_io, ONLY: read_mos_restart_low
52  USE qs_mo_types, ONLY: deallocate_mo_set,&
53  mo_set_type
54  USE rt_propagation_types, ONLY: get_rtp,&
55  rt_prop_type
56 #include "./../base/base_uses.f90"
57 
58  IMPLICIT NONE
59  PRIVATE
60 
61  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_projection_mo_utils'
62 
63  PUBLIC :: init_mo_projection, compute_and_write_proj_mo
64 
65 CONTAINS
66 
67 ! **************************************************************************************************
68 !> \brief Initialize the mo projection objects for time dependent run
69 !> \param qs_env ...
70 !> \param rtp_control ...
71 !> \author Guillaume Le Breton (04.2023)
72 ! **************************************************************************************************
73  SUBROUTINE init_mo_projection(qs_env, rtp_control)
74  TYPE(qs_environment_type), POINTER :: qs_env
75  TYPE(rtp_control_type), POINTER :: rtp_control
76 
77  INTEGER :: i_rep, j_td, n_rep_val, nbr_mo_td_max, &
78  nrep, reftype
79  INTEGER, DIMENSION(:), POINTER :: tmp_ints
80  TYPE(cp_logger_type), POINTER :: logger
81  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
82  TYPE(proj_mo_type), POINTER :: proj_mo
83  TYPE(section_vals_type), POINTER :: input, print_key, proj_mo_section
84 
85  NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
86  input, proj_mo_section, print_key, mos)
87 
88  CALL get_qs_env(qs_env, input=input, mos=mos)
89 
90  proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
91 
92  ! Read the input section and load the reference MOs
93  CALL section_vals_get(proj_mo_section, n_repetition=nrep)
94  ALLOCATE (rtp_control%proj_mo_list(nrep))
95 
96  DO i_rep = 1, nrep
97  NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
98  ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
99  proj_mo => rtp_control%proj_mo_list(i_rep)%proj_mo
100 
101  CALL section_vals_val_get(proj_mo_section, "REFERENCE_TYPE", i_rep_section=i_rep, &
102  i_val=reftype)
103 
104  CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
105  c_val=proj_mo%ref_mo_file_name)
106 
107  CALL section_vals_val_get(proj_mo_section, "REF_ADD_LUMO", i_rep_section=i_rep, &
108  i_val=proj_mo%ref_nlumo)
109 
110  ! Relevent only in EMD
111  IF (.NOT. rtp_control%fixed_ions) &
112  CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, &
113  l_val=proj_mo%propagate_ref)
114 
115  IF (reftype == proj_mo_ref_scf) THEN
116  ! If no reference .wfn is provided, using the restart SCF file:
117  IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
118  CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
119  IF (n_rep_val > 0) THEN
120  CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", c_val=proj_mo%ref_mo_file_name)
121  ELSE
122  !try to read from the filename that is generated automatically from the printkey
123  print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
124  logger => cp_get_default_logger()
125  proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
126  extension=".wfn", my_local=.false.)
127  END IF
128  END IF
129 
130  CALL section_vals_val_get(proj_mo_section, "REF_MO_INDEX", i_rep_section=i_rep, &
131  i_vals=tmp_ints)
132  ALLOCATE (proj_mo%ref_mo_index, source=tmp_ints(:))
133  CALL section_vals_val_get(proj_mo_section, "REF_MO_SPIN", i_rep_section=i_rep, &
134  i_val=proj_mo%ref_mo_spin)
135 
136  ! Read the SCF mos and store the one required
137  CALL read_reference_mo_from_wfn(qs_env, proj_mo)
138 
139  ELSE IF (reftype == proj_mo_ref_xas_tdp) THEN
140  IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
141  CALL cp_abort(__location__, &
142  "Input error in DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO. "// &
143  "For REFERENCE_TYPE XAS_TDP one must define the name "// &
144  "of the .wfn file to read the reference MO from. Please define REF_MO_FILE_NAME.")
145  END IF
146  ALLOCATE (proj_mo%ref_mo_index(1))
147  ! XAS restart files contain only one excited state
148  proj_mo%ref_mo_index(1) = 1
149  proj_mo%ref_mo_spin = 1
150  ! Read XAS TDP mos
151  CALL read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref=.true.)
152 
153  END IF
154 
155  ! Initialize the other parameters related to the TD mos.
156  CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
157  l_val=proj_mo%sum_on_all_ref)
158 
159  CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
160  i_val=proj_mo%td_mo_spin)
161  IF (proj_mo%td_mo_spin > SIZE(mos)) &
162  CALL cp_abort(__location__, &
163  "You asked to project the time dependent BETA spin while the "// &
164  "real time DFT run has only one spin defined. "// &
165  "Please set TD_MO_SPIN to 1 or use UKS.")
166 
167  CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
168  i_vals=tmp_ints)
169 
170  nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
171 
172  ALLOCATE (proj_mo%td_mo_index, source=tmp_ints(:))
173  IF (proj_mo%td_mo_index(1) == -1) THEN
174  DEALLOCATE (proj_mo%td_mo_index)
175  ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
176  ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
177  DO j_td = 1, nbr_mo_td_max
178  proj_mo%td_mo_index(j_td) = j_td
179  proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
180  END DO
181  ELSE
182  ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
183  proj_mo%td_mo_occ(:) = 0.0_dp
184  DO j_td = 1, SIZE(proj_mo%td_mo_index)
185  IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) &
186  CALL cp_abort(__location__, &
187  "The MO number available in the Time Dependent run "// &
188  "is smaller than the MO number you have required in TD_MO_INDEX.")
189  proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
190  END DO
191  END IF
192 
193  CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
194  l_val=proj_mo%sum_on_all_td)
195 
196  END DO
197 
198  END SUBROUTINE init_mo_projection
199 
200 ! **************************************************************************************************
201 !> \brief Read the MO from .wfn file and store the required mos for TD projections
202 !> \param qs_env ...
203 !> \param proj_mo ...
204 !> \param xas_ref ...
205 !> \author Guillaume Le Breton (04.2023)
206 ! **************************************************************************************************
207  SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref)
208  TYPE(qs_environment_type), POINTER :: qs_env
209  TYPE(proj_mo_type), POINTER :: proj_mo
210  LOGICAL, OPTIONAL :: xas_ref
211 
212  INTEGER :: i_ref, ispin, mo_index, natom, &
213  nbr_mo_max, nbr_ref_mo, nspins, &
214  real_mo_index, restart_unit
215  LOGICAL :: is_file, my_xasref
216  TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
217  TYPE(cp_fm_type) :: mo_coeff_temp
218  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
219  TYPE(dft_control_type), POINTER :: dft_control
220  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_qs, mo_ref_temp
221  TYPE(mo_set_type), POINTER :: mo_set
222  TYPE(mp_para_env_type), POINTER :: para_env
223  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
224  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
225 
226  NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
227  mo_ref_fmstruct, matrix_s)
228 
229  my_xasref = .false.
230  IF (PRESENT(xas_ref)) my_xasref = xas_ref
231 
232  CALL get_qs_env(qs_env, &
233  qs_kind_set=qs_kind_set, &
234  particle_set=particle_set, &
235  dft_control=dft_control, &
236  matrix_s_kp=matrix_s, &
237  mos=mo_qs, &
238  para_env=para_env)
239 
240  natom = SIZE(particle_set, 1)
241 
242  nspins = SIZE(mo_qs)
243  ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved
244  IF (my_xasref .AND. nspins < 2) THEN
245  nspins = 2
246  END IF
247  ALLOCATE (mo_ref_temp(nspins))
248 
249  DO ispin = 1, nspins
250  IF (my_xasref) THEN
251  mo_set => mo_qs(1)
252  ELSE
253  mo_set => mo_qs(ispin)
254  END IF
255  mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
256  NULLIFY (mo_ref_fmstruct)
257  CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
258  ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
259  NULLIFY (mo_ref_temp(ispin)%mo_coeff)
260  ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
261  CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
262  CALL cp_fm_struct_release(mo_ref_fmstruct)
263 
264  mo_ref_temp(ispin)%nao = mo_set%nao
265  mo_ref_temp(ispin)%homo = mo_set%homo
266  mo_ref_temp(ispin)%nelectron = mo_set%nelectron
267  ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
268  ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
269  NULLIFY (mo_set)
270  END DO
271 
272 ! DO ispin = 1, nspins
273 ! CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1))
274 ! END DO
275 ! ELSE
276 ! DO ispin = 1, nspins
277 ! CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin))
278 ! END DO
279 ! END IF
280 
281  IF (para_env%is_source()) THEN
282  INQUIRE (file=trim(proj_mo%ref_mo_file_name), exist=is_file)
283  IF (.NOT. is_file) &
284  CALL cp_abort(__location__, &
285  "Reference file not found! Name of the file CP2K looked for: "//trim(proj_mo%ref_mo_file_name))
286 
287  CALL open_file(file_name=proj_mo%ref_mo_file_name, &
288  file_action="READ", &
289  file_form="UNFORMATTED", &
290  file_status="OLD", &
291  unit_number=restart_unit)
292  END IF
293 
294  CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
295  particle_set=particle_set, natom=natom, &
296  rst_unit=restart_unit)
297 
298  IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
299 
300  IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) &
301  CALL cp_abort(__location__, &
302  "You asked as reference spin the BETA one while the "// &
303  "reference .wfn file has only one spin. Use a reference .wfn "// &
304  "with 2 spins separated or set REF_MO_SPIN to 1")
305 
306  ! Store only the mos required
307  nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
308  IF (proj_mo%ref_mo_index(1) == -1) THEN
309  DEALLOCATE (proj_mo%ref_mo_index)
310  ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
311  DO i_ref = 1, nbr_mo_max
312  proj_mo%ref_mo_index(i_ref) = i_ref
313  END DO
314  ELSE
315  DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
316  IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) &
317  CALL cp_abort(__location__, &
318  "The MO number available in the Reference SCF "// &
319  "is smaller than the MO number you have required in REF_MO_INDEX.")
320  END DO
321  END IF
322  nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
323 
324  IF (nbr_ref_mo > nbr_mo_max) &
325  CALL cp_abort(__location__, &
326  "The number of reference mo is larger then the total number of available one in the .wfn file.")
327 
328  ! Store
329  ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
330  CALL cp_fm_struct_create(mo_ref_fmstruct, &
331  context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
332  nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
333  ncol_global=1)
334 
335  IF (dft_control%rtp_control%fixed_ions) &
336  CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
337 
338  DO mo_index = 1, nbr_ref_mo
339  real_mo_index = proj_mo%ref_mo_index(mo_index)
340  IF (real_mo_index > nbr_mo_max) &
341  CALL cp_abort(__location__, &
342  "One of reference mo index is larger then the total number of available mo in the .wfn file.")
343 
344  ! fill with the reference mo values
345  CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
346  IF (dft_control%rtp_control%fixed_ions) THEN
347  ! multiply with overlap matrix to save time later on: proj_mo%mo_ref is SxMO_ref
348  CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
349  ncol=1, &
350  source_start=real_mo_index, &
351  target_start=1)
352  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
353  ELSE
354  ! the AO will change with times: proj_mo%mo_ref are really the MOs coeffs
355  CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, proj_mo%mo_ref(mo_index), &
356  ncol=1, &
357  source_start=real_mo_index, &
358  target_start=1)
359  END IF
360  END DO
361 
362  ! Clean temporary variables
363  DO ispin = 1, nspins
364  CALL deallocate_mo_set(mo_ref_temp(ispin))
365  END DO
366  DEALLOCATE (mo_ref_temp)
367 
368  CALL cp_fm_struct_release(mo_ref_fmstruct)
369  IF (dft_control%rtp_control%fixed_ions) &
370  CALL cp_fm_release(mo_coeff_temp)
371 
372  END SUBROUTINE read_reference_mo_from_wfn
373 
374 ! **************************************************************************************************
375 !> \brief Compute the projection of the current MO coefficients on reference ones
376 !> and write the results.
377 !> \param qs_env ...
378 !> \param mos_new ...
379 !> \param proj_mo ...
380 !> \param n_proj ...
381 !> \author Guillaume Le Breton
382 ! **************************************************************************************************
383  SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
384  TYPE(qs_environment_type), POINTER :: qs_env
385  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
386  TYPE(proj_mo_type) :: proj_mo
387  INTEGER :: n_proj
388 
389  INTEGER :: i_ref, nbr_ref_mo, nbr_ref_td
390  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase, popu, sum_popu_ref
391  TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
392  TYPE(cp_fm_type) :: s_mo_ref
393  TYPE(cp_logger_type), POINTER :: logger
394  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
395  TYPE(dft_control_type), POINTER :: dft_control
396  TYPE(section_vals_type), POINTER :: input, print_mo_section, proj_mo_section
397 
398  NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
399 
400  logger => cp_get_default_logger()
401 
402  CALL get_qs_env(qs_env, &
403  dft_control=dft_control, &
404  input=input)
405 
406  ! The general section
407  proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
408  ! The section we are dealing in this particular subroutine call: n_proj.
409  print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
410 
411  ! Propagate the reference MO if required at each time step
412  IF (proj_mo%propagate_ref) CALL propagate_ref_mo(qs_env, proj_mo)
413 
414  ! Does not compute the projection if not the required time step
415  IF (.NOT. btest(cp_print_key_should_output(logger%iter_info, &
416  print_mo_section, ""), &
417  cp_p_file)) &
418  RETURN
419 
420  IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
421  CALL get_qs_env(qs_env, &
422  matrix_s_kp=matrix_s)
423  CALL cp_fm_struct_create(mo_ref_fmstruct, &
424  context=proj_mo%mo_ref(1)%matrix_struct%context, &
425  nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
426  ncol_global=1)
427  CALL cp_fm_create(s_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
428  END IF
429 
430  nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
431  nbr_ref_td = SIZE(proj_mo%td_mo_index)
432  ALLOCATE (popu(nbr_ref_td))
433  ALLOCATE (phase(nbr_ref_td))
434 
435  IF (proj_mo%sum_on_all_ref) THEN
436  ALLOCATE (sum_popu_ref(nbr_ref_td))
437  sum_popu_ref(:) = 0.0_dp
438  DO i_ref = 1, nbr_ref_mo
439  ! Compute SxMO_ref for the upcoming projection later on
440  IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
441  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
442  CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
443  ELSE
444  CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
445  END IF
446  sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
447  END DO
448  IF (proj_mo%sum_on_all_td) THEN
449  CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=sum(sum_popu_ref), n_proj=n_proj)
450  ELSE
451  CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
452  END IF
453  DEALLOCATE (sum_popu_ref)
454  ELSE
455  DO i_ref = 1, nbr_ref_mo
456  IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
457  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
458  CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
459  ELSE
460  CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
461  END IF
462  IF (proj_mo%sum_on_all_td) THEN
463  CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=sum(popu), n_proj=n_proj)
464  ELSE
465 
466  CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
467  END IF
468  END DO
469  END IF
470 
471  IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
472  CALL cp_fm_struct_release(mo_ref_fmstruct)
473  CALL cp_fm_release(s_mo_ref)
474  END IF
475  DEALLOCATE (popu)
476  DEALLOCATE (phase)
477 
478  END SUBROUTINE compute_and_write_proj_mo
479 
480 ! **************************************************************************************************
481 !> \brief Compute the projection of the current MO coefficients on reference ones
482 !> \param popu ...
483 !> \param phase ...
484 !> \param mos_new ...
485 !> \param proj_mo ...
486 !> \param i_ref ...
487 !> \param S_mo_ref ...
488 !> \author Guillaume Le Breton
489 ! **************************************************************************************************
490  SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref)
491  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: popu, phase
492  TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
493  TYPE(proj_mo_type) :: proj_mo
494  INTEGER :: i_ref
495  TYPE(cp_fm_type), OPTIONAL :: s_mo_ref
496 
497  CHARACTER(len=*), PARAMETER :: routinen = 'compute_proj_mo'
498 
499  INTEGER :: handle, j_td, nbr_ref_td, spin_td
500  LOGICAL :: is_emd
501  REAL(kind=dp) :: imag_proj, real_proj
502  TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
503  TYPE(cp_fm_type) :: mo_coeff_temp
504 
505  CALL timeset(routinen, handle)
506 
507  is_emd = .false.
508  IF (PRESENT(s_mo_ref)) is_emd = .true.
509 
510  nbr_ref_td = SIZE(popu)
511  spin_td = proj_mo%td_mo_spin
512 
513  CALL cp_fm_struct_create(mo_ref_fmstruct, &
514  context=mos_new(1)%matrix_struct%context, &
515  nrow_global=mos_new(1)%matrix_struct%nrow_global, &
516  ncol_global=1)
517  CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
518 
519  DO j_td = 1, nbr_ref_td
520  ! Real part of the projection:
521  real_proj = 0.0_dp
522  CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
523  ncol=1, &
524  source_start=proj_mo%td_mo_index(j_td), &
525  target_start=1)
526  IF (is_emd) THEN
527  ! The reference MO have to be propagated in the new basis, so the projection
528  CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, real_proj)
529  ELSE
530  ! The reference MO is time independent. proj_mo%mo_ref(i_ref) is in fact SxMO_ref already
531  CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
532  END IF
533 
534  ! Imaginary part of the projection
535  imag_proj = 0.0_dp
536  CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
537  ncol=1, &
538  source_start=proj_mo%td_mo_index(j_td), &
539  target_start=1)
540 
541  IF (is_emd) THEN
542  CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, imag_proj)
543  ELSE
544  CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
545  END IF
546 
547  ! Store the result
548  phase(j_td) = atan2(imag_proj, real_proj) ! in radians
549  popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
550  END DO
551 
552  CALL cp_fm_struct_release(mo_ref_fmstruct)
553  CALL cp_fm_release(mo_coeff_temp)
554 
555  CALL timestop(handle)
556 
557  END SUBROUTINE compute_proj_mo
558 
559 ! **************************************************************************************************
560 !> \brief Write in one file the projection of (all) the time-dependent MO coefficients
561 !> on one reference ones
562 !> \param qs_env ...
563 !> \param print_mo_section ...
564 !> \param proj_mo ...
565 !> \param i_ref ...
566 !> \param popu ...
567 !> \param phase ...
568 !> \param popu_tot ...
569 !> \param n_proj ...
570 !> \author Guillaume Le Breton
571 ! **************************************************************************************************
572  SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
573  TYPE(qs_environment_type), POINTER :: qs_env
574  TYPE(section_vals_type), POINTER :: print_mo_section
575  TYPE(proj_mo_type) :: proj_mo
576  INTEGER, OPTIONAL :: i_ref
577  REAL(kind=dp), DIMENSION(:), OPTIONAL :: popu, phase
578  REAL(kind=dp), OPTIONAL :: popu_tot
579  INTEGER, OPTIONAL :: n_proj
580 
581  CHARACTER(LEN=default_string_length) :: ext, filename
582  INTEGER :: j_td, output_unit, print_unit
583  TYPE(cp_logger_type), POINTER :: logger
584 
585  NULLIFY (logger)
586 
587  logger => cp_get_default_logger()
588  output_unit = cp_logger_get_default_io_unit(logger)
589 
590  IF (.NOT. (output_unit > 0)) RETURN
591 
592  IF (proj_mo%sum_on_all_ref) THEN
593  ext = "-"//trim(adjustl(cp_to_string(n_proj)))//"-ALL_REF.dat"
594  ELSE
595  ! Filename is update wrt the reference MO number
596  ext = "-"//trim(adjustl(cp_to_string(n_proj)))// &
597  "-REF-"// &
598  trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
599  ".dat"
600  END IF
601 
602  print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
603  extension=trim(ext))
604 
605  IF (print_unit /= output_unit) THEN
606  INQUIRE (unit=print_unit, name=filename)
607 ! WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
608 ! "PROJECTION MO", "The projection of the TD MOs is done in the file:", &
609 ! TRIM(filename)
610  WRITE (unit=print_unit, fmt="(/,(T2,A,T40,I6))") &
611  "Real time propagation step:", qs_env%sim_step
612  ELSE
613  WRITE (unit=output_unit, fmt="(/,T2,A)") "PROJECTION MO"
614  END IF
615 
616  IF (proj_mo%sum_on_all_ref) THEN
617  WRITE (print_unit, "(T3,A)") &
618  "Projection on all the required MO number from the reference file "// &
619  trim(proj_mo%ref_mo_file_name)
620  IF (proj_mo%sum_on_all_td) THEN
621  WRITE (print_unit, "(T3, A, E20.12)") &
622  "The sum over all the TD MOs population:", popu_tot
623  ELSE
624  WRITE (print_unit, "(T3,A)") &
625  "For each TD MOs required is printed: Population "
626  DO j_td = 1, SIZE(popu)
627  WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
628  END DO
629  END IF
630  ELSE
631  WRITE (print_unit, "(T3,A)") &
632  "Projection on the MO number "// &
633  trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
634  " from the reference file "// &
635  trim(proj_mo%ref_mo_file_name)
636 
637  IF (proj_mo%sum_on_all_td) THEN
638  WRITE (print_unit, "(T3, A, E20.12)") &
639  "The sum over all the TD MOs population:", popu_tot
640  ELSE
641  WRITE (print_unit, "(T3,A)") &
642  "For each TD MOs required is printed: Population & Phase [rad] "
643  DO j_td = 1, SIZE(popu)
644  WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
645  END DO
646  END IF
647  END IF
648 
649  CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
650 
651  END SUBROUTINE write_proj_mo
652 
653 ! **************************************************************************************************
654 !> \brief Propagate the reference MO in case of EMD: since the nuclei moves, the MO coeff can be
655 !> propagate to represent the same MO (because the AO move with the nuclei).
656 !> To do so, we use the same formula as for the electrons of the system, but without the
657 !> Hamiltonian:
658 !> dc^j_alpha/dt = - sum_{beta, gamma} S^{-1}_{alpha, beta} B_{beta,gamma} c^j_gamma
659 !> \param qs_env ...
660 !> \param proj_mo ...
661 !> \author Guillaume Le Breton
662 ! **************************************************************************************************
663  SUBROUTINE propagate_ref_mo(qs_env, proj_mo)
664  TYPE(qs_environment_type), POINTER :: qs_env
665  TYPE(proj_mo_type) :: proj_mo
666 
667  INTEGER :: i_ref
668  REAL(kind=dp) :: dt
669  TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
670  TYPE(cp_fm_type) :: d_mo
671  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sinvb
672  TYPE(rt_prop_type), POINTER :: rtp
673 
674  CALL get_qs_env(qs_env, rtp=rtp)
675  CALL get_rtp(rtp=rtp, sinvb=sinvb, dt=dt)
676 
677  CALL cp_fm_struct_create(mo_ref_fmstruct, &
678  context=proj_mo%mo_ref(1)%matrix_struct%context, &
679  nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
680  ncol_global=1)
681  CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
682 
683  DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
684  ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
685  CALL cp_dbcsr_sm_fm_multiply(sinvb(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
686  CALL cp_fm_scale_and_add(1.0_dp, proj_mo%mo_ref(i_ref), 1.0_dp, d_mo)
687  END DO
688 
689  CALL cp_fm_struct_release(mo_ref_fmstruct)
690  CALL cp_fm_release(d_mo)
691 
692  END SUBROUTINE propagate_ref_mo
693 
694 END MODULE rt_projection_mo_utils
695 
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_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:744
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:308
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:119
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:167
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
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:515
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:803
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:704
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:1013
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_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::proj_mo_ref_xas_tdp
integer, parameter, public proj_mo_ref_xas_tdp
Definition: input_constants.F:937
input_constants::proj_mo_ref_scf
integer, parameter, public proj_mo_ref_scf
Definition: input_constants.F:937
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:724
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:697
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:1040
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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:670
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:352
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:74
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:384
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:408