(git:ed6f26b)
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-2025 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_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
65CONTAINS
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
694END MODULE rt_projection_mo_utils
695
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
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:552
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: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:351
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:164
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: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:1063
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:950
input_constants::proj_mo_ref_scf
integer, parameter, public proj_mo_ref_scf
Definition input_constants.F:950
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:731
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:704
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:1047
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_pp, 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, 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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_control_types::proj_mo_type
Definition cp_control_types.F:51
cp_control_types::rtp_control_type
Definition cp_control_types.F:73
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
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:113
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:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76