Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Function related to MO projection in RTP calculations
10 : !> \author Guillaume Le Breton 04.2023
11 : ! **************************************************************************************************
12 : MODULE rt_projection_mo_utils
13 : USE cp_control_types, ONLY: dft_control_type,&
14 : proj_mo_type,&
15 : rtp_control_type
16 : USE cp_dbcsr_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 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Initialize the mo projection objects for time dependent run
69 : !> \param qs_env ...
70 : !> \param rtp_control ...
71 : !> \author Guillaume Le Breton (04.2023)
72 : ! **************************************************************************************************
73 4 : 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 4 : INTEGER, DIMENSION(:), POINTER :: tmp_ints
80 : TYPE(cp_logger_type), POINTER :: logger
81 4 : 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 4 : NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
86 4 : input, proj_mo_section, print_key, mos)
87 :
88 4 : CALL get_qs_env(qs_env, input=input, mos=mos)
89 :
90 4 : 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 4 : CALL section_vals_get(proj_mo_section, n_repetition=nrep)
94 64 : ALLOCATE (rtp_control%proj_mo_list(nrep))
95 :
96 56 : DO i_rep = 1, nrep
97 52 : NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
98 52 : ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
99 52 : 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 52 : i_val=reftype)
103 :
104 : CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
105 52 : 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 52 : i_val=proj_mo%ref_nlumo)
109 :
110 : ! Relevent only in EMD
111 52 : IF (.NOT. rtp_control%fixed_ions) &
112 : CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, &
113 28 : l_val=proj_mo%propagate_ref)
114 :
115 52 : IF (reftype == proj_mo_ref_scf) THEN
116 : ! If no reference .wfn is provided, using the restart SCF file:
117 46 : IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
118 38 : CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
119 38 : IF (n_rep_val > 0) THEN
120 0 : 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 38 : print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
124 38 : logger => cp_get_default_logger()
125 : proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
126 38 : 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 46 : i_vals=tmp_ints)
132 194 : 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 46 : i_val=proj_mo%ref_mo_spin)
135 :
136 : ! Read the SCF mos and store the one required
137 46 : CALL read_reference_mo_from_wfn(qs_env, proj_mo)
138 :
139 6 : ELSE IF (reftype == proj_mo_ref_xas_tdp) THEN
140 6 : 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 0 : "of the .wfn file to read the reference MO from. Please define REF_MO_FILE_NAME.")
145 : END IF
146 6 : ALLOCATE (proj_mo%ref_mo_index(1))
147 : ! XAS restart files contain only one excited state
148 6 : proj_mo%ref_mo_index(1) = 1
149 6 : proj_mo%ref_mo_spin = 1
150 : ! Read XAS TDP mos
151 6 : 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 52 : 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 52 : i_val=proj_mo%td_mo_spin)
161 52 : 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 0 : "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 52 : i_vals=tmp_ints)
169 :
170 52 : nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
171 :
172 218 : ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:))
173 52 : IF (proj_mo%td_mo_index(1) == -1) THEN
174 24 : DEALLOCATE (proj_mo%td_mo_index)
175 72 : ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
176 72 : ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
177 168 : DO j_td = 1, nbr_mo_td_max
178 144 : proj_mo%td_mo_index(j_td) = j_td
179 168 : 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 84 : ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
183 66 : proj_mo%td_mo_occ(:) = 0.0_dp
184 66 : DO j_td = 1, SIZE(proj_mo%td_mo_index)
185 38 : 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 0 : "is smaller than the MO number you have required in TD_MO_INDEX.")
189 66 : 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 108 : l_val=proj_mo%sum_on_all_td)
195 :
196 : END DO
197 :
198 8 : 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 52 : 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 52 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
219 : TYPE(dft_control_type), POINTER :: dft_control
220 52 : 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 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
224 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
225 :
226 52 : NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
227 52 : mo_ref_fmstruct, matrix_s)
228 :
229 52 : my_xasref = .FALSE.
230 6 : 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 52 : para_env=para_env)
239 :
240 52 : natom = SIZE(particle_set, 1)
241 :
242 52 : nspins = SIZE(mo_qs)
243 : ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved
244 52 : IF (my_xasref .AND. nspins < 2) THEN
245 0 : nspins = 2
246 : END IF
247 260 : ALLOCATE (mo_ref_temp(nspins))
248 :
249 156 : DO ispin = 1, nspins
250 104 : IF (my_xasref) THEN
251 12 : mo_set => mo_qs(1)
252 : ELSE
253 92 : mo_set => mo_qs(ispin)
254 : END IF
255 104 : mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
256 104 : NULLIFY (mo_ref_fmstruct)
257 : CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
258 104 : ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
259 104 : NULLIFY (mo_ref_temp(ispin)%mo_coeff)
260 104 : ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
261 104 : CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
262 104 : CALL cp_fm_struct_release(mo_ref_fmstruct)
263 :
264 104 : mo_ref_temp(ispin)%nao = mo_set%nao
265 104 : mo_ref_temp(ispin)%homo = mo_set%homo
266 104 : mo_ref_temp(ispin)%nelectron = mo_set%nelectron
267 312 : ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
268 312 : ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
269 156 : 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 52 : IF (para_env%is_source()) THEN
282 26 : INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file)
283 26 : IF (.NOT. is_file) &
284 : CALL cp_abort(__LOCATION__, &
285 0 : "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 26 : 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 52 : rst_unit=restart_unit)
297 :
298 52 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
299 :
300 52 : 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 0 : "with 2 spins separated or set REF_MO_SPIN to 1")
305 :
306 : ! Store only the mos required
307 52 : nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
308 52 : IF (proj_mo%ref_mo_index(1) == -1) THEN
309 18 : DEALLOCATE (proj_mo%ref_mo_index)
310 54 : ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
311 126 : DO i_ref = 1, nbr_mo_max
312 126 : proj_mo%ref_mo_index(i_ref) = i_ref
313 : END DO
314 : ELSE
315 78 : DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
316 44 : 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 34 : "is smaller than the MO number you have required in REF_MO_INDEX.")
320 : END DO
321 : END IF
322 52 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
323 :
324 52 : IF (nbr_ref_mo > nbr_mo_max) &
325 : CALL cp_abort(__LOCATION__, &
326 0 : "The number of reference mo is larger then the total number of available one in the .wfn file.")
327 :
328 : ! Store
329 308 : 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 52 : ncol_global=1)
334 :
335 52 : IF (dft_control%rtp_control%fixed_ions) &
336 24 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
337 :
338 204 : DO mo_index = 1, nbr_ref_mo
339 152 : real_mo_index = proj_mo%ref_mo_index(mo_index)
340 152 : IF (real_mo_index > nbr_mo_max) &
341 : CALL cp_abort(__LOCATION__, &
342 0 : "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 152 : CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
346 204 : 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 64 : target_start=1)
352 64 : 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 88 : target_start=1)
359 : END IF
360 : END DO
361 :
362 : ! Clean temporary variables
363 156 : DO ispin = 1, nspins
364 156 : CALL deallocate_mo_set(mo_ref_temp(ispin))
365 : END DO
366 52 : DEALLOCATE (mo_ref_temp)
367 :
368 52 : CALL cp_fm_struct_release(mo_ref_fmstruct)
369 52 : IF (dft_control%rtp_control%fixed_ions) &
370 24 : CALL cp_fm_release(mo_coeff_temp)
371 :
372 52 : 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 104 : 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 104 : 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 104 : 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 104 : NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
399 :
400 208 : logger => cp_get_default_logger()
401 :
402 : CALL get_qs_env(qs_env, &
403 : dft_control=dft_control, &
404 104 : input=input)
405 :
406 : ! The general section
407 104 : 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 104 : 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 104 : 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 104 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
416 : print_mo_section, ""), &
417 : cp_p_file)) &
418 : RETURN
419 :
420 102 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
421 : CALL get_qs_env(qs_env, &
422 56 : 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 56 : ncol_global=1)
427 56 : CALL cp_fm_create(S_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
428 : END IF
429 :
430 102 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
431 102 : nbr_ref_td = SIZE(proj_mo%td_mo_index)
432 306 : ALLOCATE (popu(nbr_ref_td))
433 204 : ALLOCATE (phase(nbr_ref_td))
434 :
435 102 : IF (proj_mo%sum_on_all_ref) THEN
436 48 : ALLOCATE (sum_popu_ref(nbr_ref_td))
437 168 : sum_popu_ref(:) = 0.0_dp
438 168 : DO i_ref = 1, nbr_ref_mo
439 : ! Compute SxMO_ref for the upcoming projection later on
440 144 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
441 96 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
442 96 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
443 : ELSE
444 48 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
445 : END IF
446 1032 : sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
447 : END DO
448 24 : IF (proj_mo%sum_on_all_td) THEN
449 84 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=SUM(sum_popu_ref), n_proj=n_proj)
450 : ELSE
451 12 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
452 : END IF
453 24 : DEALLOCATE (sum_popu_ref)
454 : ELSE
455 234 : DO i_ref = 1, nbr_ref_mo
456 156 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
457 80 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
458 80 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
459 : ELSE
460 76 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
461 : END IF
462 234 : IF (proj_mo%sum_on_all_td) THEN
463 588 : 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 72 : 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 102 : IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
472 56 : CALL cp_fm_struct_release(mo_ref_fmstruct)
473 56 : CALL cp_fm_release(S_mo_ref)
474 : END IF
475 102 : DEALLOCATE (popu)
476 102 : DEALLOCATE (phase)
477 :
478 208 : 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 600 : 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 300 : CALL timeset(routineN, handle)
506 :
507 300 : is_emd = .FALSE.
508 300 : IF (PRESENT(S_mo_ref)) is_emd = .TRUE.
509 :
510 300 : nbr_ref_td = SIZE(popu)
511 300 : 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 300 : ncol_global=1)
517 300 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
518 :
519 1784 : 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 1484 : target_start=1)
526 1484 : IF (is_emd) THEN
527 : ! The reference MO have to be propagated in the new basis, so the projection
528 936 : 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 548 : 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 1484 : target_start=1)
540 :
541 1484 : IF (is_emd) THEN
542 936 : CALL cp_fm_trace(mo_coeff_temp, S_mo_ref, imag_proj)
543 : ELSE
544 548 : CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
545 : END IF
546 :
547 : ! Store the result
548 1484 : phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians
549 1784 : popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
550 : END DO
551 :
552 300 : CALL cp_fm_struct_release(mo_ref_fmstruct)
553 300 : CALL cp_fm_release(mo_coeff_temp)
554 :
555 300 : CALL timestop(handle)
556 :
557 300 : 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 180 : 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 180 : NULLIFY (logger)
586 :
587 180 : logger => cp_get_default_logger()
588 180 : output_unit = cp_logger_get_default_io_unit(logger)
589 :
590 180 : IF (.NOT. (output_unit > 0)) RETURN
591 :
592 90 : IF (proj_mo%sum_on_all_ref) THEN
593 12 : 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 78 : ".dat"
600 : END IF
601 :
602 : print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
603 90 : extension=TRIM(ext))
604 :
605 90 : IF (print_unit /= output_unit) THEN
606 90 : 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 90 : "Real time propagation step:", qs_env%sim_step
612 : ELSE
613 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "PROJECTION MO"
614 : END IF
615 :
616 90 : 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 12 : TRIM(proj_mo%ref_mo_file_name)
620 12 : IF (proj_mo%sum_on_all_td) THEN
621 : WRITE (print_unit, "(T3, A, E20.12)") &
622 6 : "The sum over all the TD MOs population:", popu_tot
623 : ELSE
624 : WRITE (print_unit, "(T3,A)") &
625 6 : "For each TD MOs required is printed: Population "
626 42 : DO j_td = 1, SIZE(popu)
627 42 : 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 78 : TRIM(proj_mo%ref_mo_file_name)
636 :
637 78 : IF (proj_mo%sum_on_all_td) THEN
638 : WRITE (print_unit, "(T3, A, E20.12)") &
639 42 : "The sum over all the TD MOs population:", popu_tot
640 : ELSE
641 : WRITE (print_unit, "(T3,A)") &
642 36 : "For each TD MOs required is printed: Population & Phase [rad] "
643 94 : DO j_td = 1, SIZE(popu)
644 94 : 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 90 : 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 84 : 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 28 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: SinvB
672 : TYPE(rt_prop_type), POINTER :: rtp
673 :
674 28 : CALL get_qs_env(qs_env, rtp=rtp)
675 28 : 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 28 : ncol_global=1)
681 28 : CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
682 :
683 116 : DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
684 : ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
685 88 : CALL cp_dbcsr_sm_fm_multiply(SinvB(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
686 116 : 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 28 : CALL cp_fm_struct_release(mo_ref_fmstruct)
690 28 : CALL cp_fm_release(d_mo)
691 :
692 28 : END SUBROUTINE propagate_ref_mo
693 :
694 : END MODULE rt_projection_mo_utils
695 :
|
