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