(git:a1d4336)
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) &
107 CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, &
108 l_val=proj_mo%propagate_ref)
109
110 ! If no reference .wfn is provided, using the restart SCF file:
111 IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
112 CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
113 IF (n_rep_val > 0) THEN
114 CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", c_val=proj_mo%ref_mo_file_name)
115 ELSE
116 !try to read from the filename that is generated automatically from the printkey
117 print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
118 logger => cp_get_default_logger()
119 proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
120 extension=".wfn", my_local=.false.)
121 END IF
122 END IF
123
124 CALL section_vals_val_get(proj_mo_section, "REF_MO_INDEX", i_rep_section=i_rep, &
125 i_vals=tmp_ints)
126 ALLOCATE (proj_mo%ref_mo_index, source=tmp_ints(:))
127 CALL section_vals_val_get(proj_mo_section, "REF_MO_SPIN", i_rep_section=i_rep, &
128 i_val=proj_mo%ref_mo_spin)
129
130 ! Read the SCF mos and store the one required
131 CALL read_reference_mo_from_wfn(qs_env, proj_mo)
132
133 ! Initialize the other parameters related to the TD mos.
134 CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
135 l_val=proj_mo%sum_on_all_ref)
136
137 CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
138 i_val=proj_mo%td_mo_spin)
139 IF (proj_mo%td_mo_spin > SIZE(mos)) &
140 CALL cp_abort(__location__, &
141 "You asked to project the time dependent BETA spin while the "// &
142 "real time DFT run has only one spin defined. "// &
143 "Please set TD_MO_SPIN to 1 or use UKS.")
144
145 CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
146 i_vals=tmp_ints)
147
148 nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
149
150 ALLOCATE (proj_mo%td_mo_index, source=tmp_ints(:))
151 IF (proj_mo%td_mo_index(1) == -1) THEN
152 DEALLOCATE (proj_mo%td_mo_index)
153 ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
154 ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
155 DO j_td = 1, nbr_mo_td_max
156 proj_mo%td_mo_index(j_td) = j_td
157 proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
158 END DO
159 ELSE
160 ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
161 proj_mo%td_mo_occ(:) = 0.0_dp
162 DO j_td = 1, SIZE(proj_mo%td_mo_index)
163 IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) &
164 CALL cp_abort(__location__, &
165 "The MO number available in the Time Dependent run "// &
166 "is smaller than the MO number you have required in TD_MO_INDEX.")
167 proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
168 END DO
169 END IF
170
171 CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
172 l_val=proj_mo%sum_on_all_td)
173
174 END DO
175
176 END SUBROUTINE init_mo_projection
177
178! **************************************************************************************************
179!> \brief Read the MO from .wfn file and store the required MOs for TD projections
180!> \param qs_env ...
181!> \param proj_mo ...
182!> \author Guillaume Le Breton (04.2023)
183! **************************************************************************************************
184 SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo)
185 TYPE(qs_environment_type), POINTER :: qs_env
186 TYPE(proj_mo_type), POINTER :: proj_mo
187
188 INTEGER :: i_ref, ispin, mo_index, natom, &
189 nbr_mo_max, nbr_ref_mo, nspins, &
190 real_mo_index, restart_unit
191 LOGICAL :: is_file
192 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
193 TYPE(cp_fm_type) :: mo_coeff_temp
194 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
195 TYPE(dft_control_type), POINTER :: dft_control
196 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_qs, mo_ref_temp
197 TYPE(mo_set_type), POINTER :: mo_set
198 TYPE(mp_para_env_type), POINTER :: para_env
199 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
200 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
201
202 NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
203 mo_ref_fmstruct, matrix_s)
204
205 CALL get_qs_env(qs_env, &
206 qs_kind_set=qs_kind_set, &
207 particle_set=particle_set, &
208 dft_control=dft_control, &
209 matrix_s_kp=matrix_s, &
210 mos=mo_qs, &
211 para_env=para_env)
212
213 natom = SIZE(particle_set, 1)
214
215 nspins = SIZE(mo_qs)
216
217 ALLOCATE (mo_ref_temp(nspins))
218
219 DO ispin = 1, nspins
220 mo_set => mo_qs(ispin)
221 mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
222 NULLIFY (mo_ref_fmstruct)
223 CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
224 ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
225 NULLIFY (mo_ref_temp(ispin)%mo_coeff)
226 ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
227 CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
228 CALL cp_fm_struct_release(mo_ref_fmstruct)
229
230 mo_ref_temp(ispin)%nao = mo_set%nao
231 mo_ref_temp(ispin)%homo = mo_set%homo
232 mo_ref_temp(ispin)%nelectron = mo_set%nelectron
233 ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
234 ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
235 NULLIFY (mo_set)
236 END DO
237
238 IF (para_env%is_source()) THEN
239 INQUIRE (file=trim(proj_mo%ref_mo_file_name), exist=is_file)
240 IF (.NOT. is_file) &
241 CALL cp_abort(__location__, &
242 "Reference file not found! Name of the file CP2K looked for: "//trim(proj_mo%ref_mo_file_name))
243
244 CALL open_file(file_name=proj_mo%ref_mo_file_name, &
245 file_action="READ", &
246 file_form="UNFORMATTED", &
247 file_status="OLD", &
248 unit_number=restart_unit)
249 END IF
250
251 CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
252 particle_set=particle_set, natom=natom, &
253 rst_unit=restart_unit)
254
255 IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
256
257 IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) &
258 CALL cp_abort(__location__, &
259 "Projection on spin BETA is not possible as the reference wavefunction "// &
260 "only has one spin channel. Use a reference .wfn calculated with UKS/LSD, or set REF_MO_SPIN to 1")
261
262 ! Store only the mos required
263 nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
264 IF (proj_mo%ref_mo_index(1) == -1) THEN
265 DEALLOCATE (proj_mo%ref_mo_index)
266 ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
267 DO i_ref = 1, nbr_mo_max
268 proj_mo%ref_mo_index(i_ref) = i_ref
269 END DO
270 ELSE
271 DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
272 IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) &
273 CALL cp_abort(__location__, &
274 "The number of MOs available in the reference wavefunction "// &
275 "is smaller than the MO number you have requested in REF_MO_INDEX.")
276 END DO
277 END IF
278 nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
279
280 IF (nbr_ref_mo > nbr_mo_max) &
281 CALL cp_abort(__location__, &
282 "The total number of requested MOs is larger than what is available in the reference wavefunction. "// &
283 "If you are trying to project onto virtual states, make sure they are included in the .wfn file "// &
284 "e.g., by the ADDED_MOS keyword in the SCF section of the input when calculating your reference.")
285
286 ! Store
287 ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
288 CALL cp_fm_struct_create(mo_ref_fmstruct, &
289 context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
290 nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
291 ncol_global=1)
292
293 IF (dft_control%rtp_control%fixed_ions) &
294 CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
295
296 DO mo_index = 1, nbr_ref_mo
297 real_mo_index = proj_mo%ref_mo_index(mo_index)
298 IF (real_mo_index > nbr_mo_max) &
299 CALL cp_abort(__location__, &
300 "One of reference mo index is larger then the total number of available mo in the .wfn file.")
301
302 ! fill with the reference mo values
303 CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
304 IF (dft_control%rtp_control%fixed_ions) THEN
305 ! multiply with overlap matrix to save time later on: proj_mo%mo_ref is SxMO_ref
306 CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
307 ncol=1, &
308 source_start=real_mo_index, &
309 target_start=1)
310 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
311 ELSE
312 ! the AO will change with times: proj_mo%mo_ref are really the MOs coeffs
313 CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, proj_mo%mo_ref(mo_index), &
314 ncol=1, &
315 source_start=real_mo_index, &
316 target_start=1)
317 END IF
318 END DO
319
320 ! Clean temporary variables
321 DO ispin = 1, nspins
322 CALL deallocate_mo_set(mo_ref_temp(ispin))
323 END DO
324 DEALLOCATE (mo_ref_temp)
325
326 CALL cp_fm_struct_release(mo_ref_fmstruct)
327 IF (dft_control%rtp_control%fixed_ions) &
328 CALL cp_fm_release(mo_coeff_temp)
329
330 END SUBROUTINE read_reference_mo_from_wfn
331
332! **************************************************************************************************
333!> \brief Compute the projection of the current MO coefficients on reference ones
334!> and write the results.
335!> \param qs_env ...
336!> \param mos_new ...
337!> \param proj_mo ...
338!> \param n_proj ...
339!> \author Guillaume Le Breton
340! **************************************************************************************************
341 SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
342 TYPE(qs_environment_type), POINTER :: qs_env
343 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
344 TYPE(proj_mo_type) :: proj_mo
345 INTEGER :: n_proj
346
347 INTEGER :: i_ref, nbr_ref_mo, nbr_ref_td
348 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase, popu, sum_popu_ref
349 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
350 TYPE(cp_fm_type) :: s_mo_ref
351 TYPE(cp_logger_type), POINTER :: logger
352 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
353 TYPE(dft_control_type), POINTER :: dft_control
354 TYPE(section_vals_type), POINTER :: input, print_mo_section, proj_mo_section
355
356 NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
357
358 logger => cp_get_default_logger()
359
360 CALL get_qs_env(qs_env, &
361 dft_control=dft_control, &
362 input=input)
363
364 ! The general section
365 proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
366 ! The section we are dealing in this particular subroutine call: n_proj.
367 print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
368
369 ! Propagate the reference MO if required at each time step
370 IF (proj_mo%propagate_ref) CALL propagate_ref_mo(qs_env, proj_mo)
371
372 ! Does not compute the projection if not the required time step
373 IF (.NOT. btest(cp_print_key_should_output(logger%iter_info, &
374 print_mo_section, ""), &
375 cp_p_file)) &
376 RETURN
377
378 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
379 CALL get_qs_env(qs_env, &
380 matrix_s_kp=matrix_s)
381 CALL cp_fm_struct_create(mo_ref_fmstruct, &
382 context=proj_mo%mo_ref(1)%matrix_struct%context, &
383 nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
384 ncol_global=1)
385 CALL cp_fm_create(s_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
386 END IF
387
388 nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
389 nbr_ref_td = SIZE(proj_mo%td_mo_index)
390 ALLOCATE (popu(nbr_ref_td))
391 ALLOCATE (phase(nbr_ref_td))
392
393 IF (proj_mo%sum_on_all_ref) THEN
394 ALLOCATE (sum_popu_ref(nbr_ref_td))
395 sum_popu_ref(:) = 0.0_dp
396 DO i_ref = 1, nbr_ref_mo
397 ! Compute SxMO_ref for the upcoming projection later on
398 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
399 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
400 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
401 ELSE
402 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
403 END IF
404 sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
405 END DO
406 IF (proj_mo%sum_on_all_td) THEN
407 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=sum(sum_popu_ref), n_proj=n_proj)
408 ELSE
409 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
410 END IF
411 DEALLOCATE (sum_popu_ref)
412 ELSE
413 DO i_ref = 1, nbr_ref_mo
414 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
415 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), s_mo_ref, ncol=1)
416 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, s_mo_ref=s_mo_ref)
417 ELSE
418 CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
419 END IF
420 IF (proj_mo%sum_on_all_td) THEN
421 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=sum(popu), n_proj=n_proj)
422 ELSE
423
424 CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
425 END IF
426 END DO
427 END IF
428
429 IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
430 CALL cp_fm_struct_release(mo_ref_fmstruct)
431 CALL cp_fm_release(s_mo_ref)
432 END IF
433 DEALLOCATE (popu)
434 DEALLOCATE (phase)
435
436 END SUBROUTINE compute_and_write_proj_mo
437
438! **************************************************************************************************
439!> \brief Compute the projection of the current MO coefficients on reference ones
440!> \param popu ...
441!> \param phase ...
442!> \param mos_new ...
443!> \param proj_mo ...
444!> \param i_ref ...
445!> \param S_mo_ref ...
446!> \author Guillaume Le Breton
447! **************************************************************************************************
448 SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref)
449 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: popu, phase
450 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
451 TYPE(proj_mo_type) :: proj_mo
452 INTEGER :: i_ref
453 TYPE(cp_fm_type), OPTIONAL :: s_mo_ref
454
455 CHARACTER(len=*), PARAMETER :: routinen = 'compute_proj_mo'
456
457 INTEGER :: handle, j_td, nbr_ref_td, spin_td
458 LOGICAL :: is_emd
459 REAL(kind=dp) :: imag_proj, real_proj
460 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
461 TYPE(cp_fm_type) :: mo_coeff_temp
462
463 CALL timeset(routinen, handle)
464
465 is_emd = .false.
466 IF (PRESENT(s_mo_ref)) is_emd = .true.
467
468 nbr_ref_td = SIZE(popu)
469 spin_td = proj_mo%td_mo_spin
470
471 CALL cp_fm_struct_create(mo_ref_fmstruct, &
472 context=mos_new(1)%matrix_struct%context, &
473 nrow_global=mos_new(1)%matrix_struct%nrow_global, &
474 ncol_global=1)
475 CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
476
477 DO j_td = 1, nbr_ref_td
478 ! Real part of the projection:
479 real_proj = 0.0_dp
480 CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
481 ncol=1, &
482 source_start=proj_mo%td_mo_index(j_td), &
483 target_start=1)
484 IF (is_emd) THEN
485 ! The reference MO have to be propagated in the new basis, so the projection
486 CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, real_proj)
487 ELSE
488 ! The reference MO is time independent. proj_mo%mo_ref(i_ref) is in fact SxMO_ref already
489 CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
490 END IF
491
492 ! Imaginary part of the projection
493 imag_proj = 0.0_dp
494 CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
495 ncol=1, &
496 source_start=proj_mo%td_mo_index(j_td), &
497 target_start=1)
498
499 IF (is_emd) THEN
500 CALL cp_fm_trace(mo_coeff_temp, s_mo_ref, imag_proj)
501 ELSE
502 CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
503 END IF
504
505 ! Store the result
506 phase(j_td) = atan2(imag_proj, real_proj) ! in radians
507 popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
508 END DO
509
510 CALL cp_fm_struct_release(mo_ref_fmstruct)
511 CALL cp_fm_release(mo_coeff_temp)
512
513 CALL timestop(handle)
514
515 END SUBROUTINE compute_proj_mo
516
517! **************************************************************************************************
518!> \brief Write in one file the projection of (all) the time-dependent MO coefficients
519!> onto reference ones
520!> \param qs_env ...
521!> \param print_mo_section ...
522!> \param proj_mo ...
523!> \param i_ref ...
524!> \param popu ...
525!> \param phase ...
526!> \param popu_tot ...
527!> \param n_proj ...
528!> \author Guillaume Le Breton
529! **************************************************************************************************
530 SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
531 TYPE(qs_environment_type), POINTER :: qs_env
532 TYPE(section_vals_type), POINTER :: print_mo_section
533 TYPE(proj_mo_type) :: proj_mo
534 INTEGER, OPTIONAL :: i_ref
535 REAL(kind=dp), DIMENSION(:), OPTIONAL :: popu, phase
536 REAL(kind=dp), OPTIONAL :: popu_tot
537 INTEGER, OPTIONAL :: n_proj
538
539 CHARACTER(LEN=default_string_length) :: ext, filename
540 INTEGER :: j_td, output_unit, print_unit
541 TYPE(cp_logger_type), POINTER :: logger
542
543 NULLIFY (logger)
544
545 logger => cp_get_default_logger()
546 output_unit = cp_logger_get_default_io_unit(logger)
547
548 IF (.NOT. (output_unit > 0)) RETURN
549
550 IF (proj_mo%sum_on_all_ref) THEN
551 ext = "-"//trim(adjustl(cp_to_string(n_proj)))//"-ALL_REF.dat"
552 ELSE
553 ! Filename is updated wrt the reference MO number
554 ext = "-"//trim(adjustl(cp_to_string(n_proj)))// &
555 "-REF-"// &
556 trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
557 ".dat"
558 END IF
559
560 print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
561 extension=trim(ext))
562
563 IF (print_unit /= output_unit) THEN
564 INQUIRE (unit=print_unit, name=filename)
565 WRITE (unit=print_unit, fmt="(/,(T2,A,T40,I6))") &
566 "Real time propagation step:", qs_env%sim_step
567 ELSE
568 WRITE (unit=output_unit, fmt="(/,T2,A)") "PROJECTION MO"
569 END IF
570
571 IF (proj_mo%sum_on_all_ref) THEN
572 WRITE (print_unit, "(T3,A)") &
573 "Projection on all the required MO number from the reference file "// &
574 trim(proj_mo%ref_mo_file_name)
575 IF (proj_mo%sum_on_all_td) THEN
576 WRITE (print_unit, "(T3, A, E20.12)") &
577 "The sum over all the TD MOs population:", popu_tot
578 ELSE
579 WRITE (print_unit, "(T3,A)") &
580 "For each TD MOs required is printed: Population "
581 DO j_td = 1, SIZE(popu)
582 WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
583 END DO
584 END IF
585 ELSE
586 WRITE (print_unit, "(T3,A)") &
587 "Projection on the MO number "// &
588 trim(adjustl(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
589 " from the reference file "// &
590 trim(proj_mo%ref_mo_file_name)
591
592 IF (proj_mo%sum_on_all_td) THEN
593 WRITE (print_unit, "(T3, A, E20.12)") &
594 "The sum over all the TD MOs population:", popu_tot
595 ELSE
596 WRITE (print_unit, "(T3,A)") &
597 "For each TD MOs required is printed: Population & Phase [rad] "
598 DO j_td = 1, SIZE(popu)
599 WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
600 END DO
601 END IF
602 END IF
603
604 CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
605
606 END SUBROUTINE write_proj_mo
607
608! **************************************************************************************************
609!> \brief Propagate the reference MO in case of EMD: since the nuclei moves, the MO coeff can be
610!> propagated to represent the same MO (because the AO move with the nuclei).
611!> To do so, we use the same formula as for the electrons of the system, but without the
612!> Hamiltonian:
613!> dc^j_alpha/dt = - sum_{beta, gamma} S^{-1}_{alpha, beta} B_{beta,gamma} c^j_gamma
614!> \param qs_env ...
615!> \param proj_mo ...
616!> \author Guillaume Le Breton
617! **************************************************************************************************
618 SUBROUTINE propagate_ref_mo(qs_env, proj_mo)
619 TYPE(qs_environment_type), POINTER :: qs_env
620 TYPE(proj_mo_type) :: proj_mo
621
622 INTEGER :: i_ref
623 REAL(kind=dp) :: dt
624 TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
625 TYPE(cp_fm_type) :: d_mo
626 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sinvb
627 TYPE(rt_prop_type), POINTER :: rtp
628
629 CALL get_qs_env(qs_env, rtp=rtp)
630 CALL get_rtp(rtp=rtp, sinvb=sinvb, dt=dt)
631
632 CALL cp_fm_struct_create(mo_ref_fmstruct, &
633 context=proj_mo%mo_ref(1)%matrix_struct%context, &
634 nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
635 ncol_global=1)
636 CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
637
638 DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
639 ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
640 CALL cp_dbcsr_sm_fm_multiply(sinvb(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
641 CALL cp_fm_scale_and_add(1.0_dp, proj_mo%mo_ref(i_ref), 1.0_dp, d_mo)
642 END DO
643
644 CALL cp_fm_struct_release(mo_ref_fmstruct)
645 CALL cp_fm_release(d_mo)
646
647 END SUBROUTINE propagate_ref_mo
648
649END MODULE rt_projection_mo_utils
650
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:80
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
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: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:172
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, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
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_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, 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, 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:530
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:687
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:354
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:342
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:633
cp_control_types::proj_mo_type
Definition cp_control_types.F:69
cp_control_types::rtp_control_type
Definition cp_control_types.F:91
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:115
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:743
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76