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) &
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 :
110 : ! If no reference .wfn is provided, using the restart SCF file:
111 46 : IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
112 38 : CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
113 38 : IF (n_rep_val > 0) THEN
114 0 : 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 38 : print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
118 38 : logger => cp_get_default_logger()
119 : proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
120 38 : 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 46 : i_vals=tmp_ints)
126 194 : 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 46 : i_val=proj_mo%ref_mo_spin)
129 :
130 : ! Read the SCF mos and store the one required
131 46 : 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 46 : 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 46 : i_val=proj_mo%td_mo_spin)
139 46 : 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 0 : "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 46 : i_vals=tmp_ints)
147 :
148 46 : nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
149 :
150 194 : ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:))
151 46 : IF (proj_mo%td_mo_index(1) == -1) THEN
152 18 : DEALLOCATE (proj_mo%td_mo_index)
153 54 : ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
154 54 : ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
155 126 : DO j_td = 1, nbr_mo_td_max
156 108 : proj_mo%td_mo_index(j_td) = j_td
157 126 : 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 84 : ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
161 66 : proj_mo%td_mo_occ(:) = 0.0_dp
162 66 : DO j_td = 1, SIZE(proj_mo%td_mo_index)
163 38 : 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 0 : "is smaller than the MO number you have required in TD_MO_INDEX.")
167 66 : 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 50 : l_val=proj_mo%sum_on_all_td)
173 :
174 : END DO
175 :
176 8 : 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 46 : 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 46 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
195 : TYPE(dft_control_type), POINTER :: dft_control
196 46 : 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 46 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
200 46 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
201 :
202 46 : NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
203 46 : 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 46 : para_env=para_env)
212 :
213 46 : natom = SIZE(particle_set, 1)
214 :
215 46 : nspins = SIZE(mo_qs)
216 :
217 230 : ALLOCATE (mo_ref_temp(nspins))
218 :
219 138 : DO ispin = 1, nspins
220 92 : mo_set => mo_qs(ispin)
221 92 : mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
222 92 : NULLIFY (mo_ref_fmstruct)
223 : CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
224 92 : ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
225 92 : NULLIFY (mo_ref_temp(ispin)%mo_coeff)
226 92 : ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
227 92 : CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
228 92 : CALL cp_fm_struct_release(mo_ref_fmstruct)
229 :
230 92 : mo_ref_temp(ispin)%nao = mo_set%nao
231 92 : mo_ref_temp(ispin)%homo = mo_set%homo
232 92 : mo_ref_temp(ispin)%nelectron = mo_set%nelectron
233 276 : ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
234 276 : ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
235 138 : NULLIFY (mo_set)
236 : END DO
237 :
238 46 : IF (para_env%is_source()) THEN
239 23 : INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file)
240 23 : IF (.NOT. is_file) &
241 : CALL cp_abort(__LOCATION__, &
242 0 : "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 23 : 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 46 : rst_unit=restart_unit)
254 :
255 46 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
256 :
257 46 : 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 0 : "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 46 : nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
264 46 : IF (proj_mo%ref_mo_index(1) == -1) THEN
265 18 : DEALLOCATE (proj_mo%ref_mo_index)
266 54 : ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
267 126 : DO i_ref = 1, nbr_mo_max
268 126 : proj_mo%ref_mo_index(i_ref) = i_ref
269 : END DO
270 : ELSE
271 66 : DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
272 38 : 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 28 : "is smaller than the MO number you have requested in REF_MO_INDEX.")
276 : END DO
277 : END IF
278 46 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
279 :
280 46 : 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 0 : "e.g., by the ADDED_MOS keyword in the SCF section of the input when calculating your reference.")
285 :
286 : ! Store
287 284 : 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 46 : ncol_global=1)
292 :
293 46 : IF (dft_control%rtp_control%fixed_ions) &
294 22 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
295 :
296 192 : DO mo_index = 1, nbr_ref_mo
297 146 : real_mo_index = proj_mo%ref_mo_index(mo_index)
298 146 : IF (real_mo_index > nbr_mo_max) &
299 : CALL cp_abort(__LOCATION__, &
300 0 : "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 146 : CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
304 192 : 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 62 : target_start=1)
310 62 : 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 84 : target_start=1)
317 : END IF
318 : END DO
319 :
320 : ! Clean temporary variables
321 138 : DO ispin = 1, nspins
322 138 : CALL deallocate_mo_set(mo_ref_temp(ispin))
323 : END DO
324 46 : DEALLOCATE (mo_ref_temp)
325 :
326 46 : CALL cp_fm_struct_release(mo_ref_fmstruct)
327 46 : IF (dft_control%rtp_control%fixed_ions) &
328 22 : CALL cp_fm_release(mo_coeff_temp)
329 :
330 46 : 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 92 : 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 92 : 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 92 : 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 92 : NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
357 :
358 184 : logger => cp_get_default_logger()
359 :
360 : CALL get_qs_env(qs_env, &
361 : dft_control=dft_control, &
362 92 : input=input)
363 :
364 : ! The general section
365 92 : 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 92 : 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 92 : 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 92 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
374 : print_mo_section, ""), &
375 : cp_p_file)) &
376 : RETURN
377 :
378 90 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
379 : CALL get_qs_env(qs_env, &
380 48 : 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 48 : ncol_global=1)
385 48 : CALL cp_fm_create(S_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
386 : END IF
387 :
388 90 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
389 90 : nbr_ref_td = SIZE(proj_mo%td_mo_index)
390 270 : ALLOCATE (popu(nbr_ref_td))
391 180 : ALLOCATE (phase(nbr_ref_td))
392 :
393 90 : IF (proj_mo%sum_on_all_ref) THEN
394 48 : ALLOCATE (sum_popu_ref(nbr_ref_td))
395 168 : sum_popu_ref(:) = 0.0_dp
396 168 : DO i_ref = 1, nbr_ref_mo
397 : ! Compute SxMO_ref for the upcoming projection later on
398 144 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
399 96 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
400 96 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
401 : ELSE
402 48 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
403 : END IF
404 1032 : sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
405 : END DO
406 24 : IF (proj_mo%sum_on_all_td) THEN
407 84 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=SUM(sum_popu_ref), n_proj=n_proj)
408 : ELSE
409 12 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
410 : END IF
411 24 : DEALLOCATE (sum_popu_ref)
412 : ELSE
413 210 : DO i_ref = 1, nbr_ref_mo
414 144 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
415 72 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
416 72 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
417 : ELSE
418 72 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
419 : END IF
420 210 : IF (proj_mo%sum_on_all_td) THEN
421 504 : 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 72 : 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 90 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
430 48 : CALL cp_fm_struct_release(mo_ref_fmstruct)
431 48 : CALL cp_fm_release(S_mo_ref)
432 : END IF
433 90 : DEALLOCATE (popu)
434 90 : DEALLOCATE (phase)
435 :
436 184 : 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 576 : 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 288 : CALL timeset(routineN, handle)
464 :
465 288 : is_emd = .FALSE.
466 288 : IF (PRESENT(S_mo_ref)) is_emd = .TRUE.
467 :
468 288 : nbr_ref_td = SIZE(popu)
469 288 : 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 288 : ncol_global=1)
475 288 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
476 :
477 1700 : 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 1412 : target_start=1)
484 1412 : IF (is_emd) THEN
485 : ! The reference MO have to be propagated in the new basis, so the projection
486 888 : 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 524 : 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 1412 : target_start=1)
498 :
499 1412 : IF (is_emd) THEN
500 888 : CALL cp_fm_trace(mo_coeff_temp, S_mo_ref, imag_proj)
501 : ELSE
502 524 : CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
503 : END IF
504 :
505 : ! Store the result
506 1412 : phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians
507 1700 : popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
508 : END DO
509 :
510 288 : CALL cp_fm_struct_release(mo_ref_fmstruct)
511 288 : CALL cp_fm_release(mo_coeff_temp)
512 :
513 288 : CALL timestop(handle)
514 :
515 288 : 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 168 : 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 168 : NULLIFY (logger)
544 :
545 168 : logger => cp_get_default_logger()
546 168 : output_unit = cp_logger_get_default_io_unit(logger)
547 :
548 168 : IF (.NOT. (output_unit > 0)) RETURN
549 :
550 84 : IF (proj_mo%sum_on_all_ref) THEN
551 12 : 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 72 : ".dat"
558 : END IF
559 :
560 : print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
561 84 : extension=TRIM(ext))
562 :
563 84 : IF (print_unit /= output_unit) THEN
564 84 : INQUIRE (UNIT=print_unit, NAME=filename)
565 : WRITE (UNIT=print_unit, FMT="(/,(T2,A,T40,I6))") &
566 84 : "Real time propagation step:", qs_env%sim_step
567 : ELSE
568 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "PROJECTION MO"
569 : END IF
570 :
571 84 : 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 12 : TRIM(proj_mo%ref_mo_file_name)
575 12 : IF (proj_mo%sum_on_all_td) THEN
576 : WRITE (print_unit, "(T3, A, E20.12)") &
577 6 : "The sum over all the TD MOs population:", popu_tot
578 : ELSE
579 : WRITE (print_unit, "(T3,A)") &
580 6 : "For each TD MOs required is printed: Population "
581 42 : DO j_td = 1, SIZE(popu)
582 42 : 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 72 : TRIM(proj_mo%ref_mo_file_name)
591 :
592 72 : IF (proj_mo%sum_on_all_td) THEN
593 : WRITE (print_unit, "(T3, A, E20.12)") &
594 36 : "The sum over all the TD MOs population:", popu_tot
595 : ELSE
596 : WRITE (print_unit, "(T3,A)") &
597 36 : "For each TD MOs required is printed: Population & Phase [rad] "
598 94 : DO j_td = 1, SIZE(popu)
599 94 : 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 84 : 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 72 : 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 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: SinvB
627 : TYPE(rt_prop_type), POINTER :: rtp
628 :
629 24 : CALL get_qs_env(qs_env, rtp=rtp)
630 24 : 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 24 : ncol_global=1)
636 24 : CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
637 :
638 108 : DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
639 : ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
640 84 : CALL cp_dbcsr_sm_fm_multiply(SinvB(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
641 108 : 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 24 : CALL cp_fm_struct_release(mo_ref_fmstruct)
645 24 : CALL cp_fm_release(d_mo)
646 :
647 24 : END SUBROUTINE propagate_ref_mo
648 :
649 : END MODULE rt_projection_mo_utils
650 :
|
