Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 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_operations, ONLY: cp_dbcsr_sm_fm_multiply
17 : USE cp_files, ONLY: close_file,&
18 : open_file
19 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace
20 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
21 : cp_fm_struct_release,&
22 : cp_fm_struct_type
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_release,&
25 : cp_fm_to_fm,&
26 : cp_fm_type
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_get_default_io_unit,&
29 : cp_logger_type,&
30 : cp_to_string
31 : USE cp_output_handling, ONLY: cp_p_file,&
32 : cp_print_key_finished_output,&
33 : cp_print_key_generate_filename,&
34 : cp_print_key_should_output,&
35 : cp_print_key_unit_nr
36 : USE dbcsr_api, ONLY: dbcsr_p_type
37 : USE input_constants, ONLY: proj_mo_ref_scf,&
38 : proj_mo_ref_xas_tdp
39 : USE input_section_types, ONLY: section_vals_get,&
40 : section_vals_get_subs_vals,&
41 : section_vals_type,&
42 : section_vals_val_get
43 : USE kinds, ONLY: default_string_length,&
44 : dp
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE particle_types, ONLY: particle_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_kind_types, ONLY: qs_kind_type
50 : USE qs_mo_io, ONLY: read_mos_restart_low
51 : USE qs_mo_types, ONLY: deallocate_mo_set,&
52 : duplicate_mo_set,&
53 : mo_set_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 2 : 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, reftype
77 2 : INTEGER, DIMENSION(:), POINTER :: tmp_ints
78 : TYPE(cp_logger_type), POINTER :: logger
79 2 : 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 2 : NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
84 2 : input, proj_mo_section, print_key, mos)
85 :
86 2 : IF (.NOT. rtp_control%fixed_ions) &
87 : CALL cp_abort(__LOCATION__, &
88 0 : "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO not implemented for EMD yet.")
89 :
90 2 : CALL get_qs_env(qs_env, input=input, mos=mos)
91 :
92 2 : proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
93 :
94 : ! Read the input section and load the reference MOs
95 2 : CALL section_vals_get(proj_mo_section, n_repetition=nrep)
96 26 : ALLOCATE (rtp_control%proj_mo_list(nrep))
97 :
98 22 : DO i_rep = 1, nrep
99 20 : NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
100 20 : ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
101 20 : proj_mo => rtp_control%proj_mo_list(i_rep)%proj_mo
102 :
103 : CALL section_vals_val_get(proj_mo_section, "REFERENCE_TYPE", i_rep_section=i_rep, &
104 20 : i_val=reftype)
105 :
106 : CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
107 20 : c_val=proj_mo%ref_mo_file_name)
108 :
109 20 : IF (reftype == proj_mo_ref_scf) THEN
110 : ! If no reference .wfn is provided, using the restart SCF file:
111 18 : IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
112 16 : CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
113 16 : 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 16 : print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
118 16 : logger => cp_get_default_logger()
119 : proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
120 16 : 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 18 : i_vals=tmp_ints)
126 82 : 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 18 : i_val=proj_mo%ref_mo_spin)
129 :
130 : ! Read the SCF mos and store the one required
131 18 : CALL read_reference_mo_from_wfn(qs_env, proj_mo)
132 :
133 2 : ELSE IF (reftype == proj_mo_ref_xas_tdp) THEN
134 2 : IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
135 : CALL cp_abort(__LOCATION__, &
136 : "Input error in DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO. "// &
137 : "For REFERENCE_TYPE XAS_TDP one must define the name "// &
138 0 : "of the .wfn file to read the reference MO from. Please define REF_MO_FILE_NAME.")
139 : END IF
140 2 : ALLOCATE (proj_mo%ref_mo_index(1))
141 : ! XAS restart files contain only one excited state
142 2 : proj_mo%ref_mo_index(1) = 1
143 2 : proj_mo%ref_mo_spin = 1
144 : ! Read XAS TDP mos
145 2 : CALL read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref=.TRUE.)
146 :
147 : END IF
148 :
149 : ! Initialize the other parameters related to the TD mos.
150 : CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
151 20 : l_val=proj_mo%sum_on_all_ref)
152 :
153 : CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
154 20 : i_val=proj_mo%td_mo_spin)
155 20 : IF (proj_mo%td_mo_spin > SIZE(mos)) &
156 : CALL cp_abort(__LOCATION__, &
157 : "You asked to project the time dependent BETA spin while the "// &
158 : "real time DFT run has only one spin defined. "// &
159 0 : "Please set TD_MO_SPIN to 1 or use UKS.")
160 :
161 : CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
162 20 : i_vals=tmp_ints)
163 :
164 90 : ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:))
165 20 : nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
166 20 : IF (proj_mo%td_mo_index(1) == -1) THEN
167 8 : DEALLOCATE (proj_mo%td_mo_index)
168 24 : ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
169 56 : DO j_td = 1, nbr_mo_td_max
170 56 : proj_mo%td_mo_index(j_td) = j_td
171 : END DO
172 : ELSE
173 34 : DO j_td = 1, SIZE(proj_mo%td_mo_index)
174 22 : IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) &
175 : CALL cp_abort(__LOCATION__, &
176 : "The MO number available in the Time Dependent run "// &
177 12 : "is smaller than the MO number you have required in TD_MO_INDEX.")
178 : END DO
179 : END IF
180 :
181 : CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
182 42 : l_val=proj_mo%sum_on_all_td)
183 :
184 : END DO
185 :
186 4 : END SUBROUTINE init_mo_projection
187 :
188 : ! **************************************************************************************************
189 : !> \brief Read the MO from .wfn file and store the required mos for TD projections
190 : !> \param qs_env ...
191 : !> \param proj_mo ...
192 : !> \param xas_ref ...
193 : !> \author Guillaume Le Breton (04.2023)
194 : ! **************************************************************************************************
195 20 : SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref)
196 : TYPE(qs_environment_type), POINTER :: qs_env
197 : TYPE(proj_mo_type), POINTER :: proj_mo
198 : LOGICAL, OPTIONAL :: xas_ref
199 :
200 : INTEGER :: i_ref, ispin, mo_index, natom, &
201 : nbr_mo_max, nbr_ref_mo, nspins, &
202 : real_mo_index, restart_unit
203 : LOGICAL :: is_file, my_xasref
204 : TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
205 : TYPE(cp_fm_type) :: mo_coeff_temp
206 20 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
207 : TYPE(dft_control_type), POINTER :: dft_control
208 20 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_qs, mo_ref_temp
209 : TYPE(mp_para_env_type), POINTER :: para_env
210 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
211 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
212 :
213 20 : NULLIFY (mo_qs, mo_ref_temp, qs_kind_set, particle_set, para_env, dft_control, &
214 20 : mo_ref_fmstruct, matrix_s)
215 :
216 20 : my_xasref = .FALSE.
217 2 : IF (PRESENT(xas_ref)) my_xasref = xas_ref
218 :
219 : CALL get_qs_env(qs_env, &
220 : qs_kind_set=qs_kind_set, &
221 : particle_set=particle_set, &
222 : dft_control=dft_control, &
223 : matrix_s_kp=matrix_s, &
224 : mos=mo_qs, &
225 20 : para_env=para_env)
226 :
227 20 : natom = SIZE(particle_set, 1)
228 :
229 : ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved
230 20 : IF (my_xasref) THEN
231 2 : nspins = 2
232 6 : ALLOCATE (mo_ref_temp(nspins))
233 6 : DO ispin = 1, nspins
234 6 : CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1))
235 : END DO
236 : ELSE
237 18 : nspins = SIZE(mo_qs)
238 90 : ALLOCATE (mo_ref_temp(nspins))
239 54 : DO ispin = 1, nspins
240 54 : CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin))
241 : END DO
242 : END IF
243 :
244 20 : IF (para_env%is_source()) THEN
245 10 : INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file)
246 10 : IF (.NOT. is_file) &
247 : CALL cp_abort(__LOCATION__, &
248 0 : "Reference file not found! Name of the file CP2K looked for: "//TRIM(proj_mo%ref_mo_file_name))
249 :
250 : CALL open_file(file_name=proj_mo%ref_mo_file_name, &
251 : file_action="READ", &
252 : file_form="UNFORMATTED", &
253 : file_status="OLD", &
254 10 : unit_number=restart_unit)
255 : END IF
256 :
257 : CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
258 : particle_set=particle_set, natom=natom, &
259 20 : rst_unit=restart_unit)
260 :
261 20 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
262 :
263 20 : IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) &
264 : CALL cp_abort(__LOCATION__, &
265 : "You asked as reference spin the BETA one while the "// &
266 : "reference .wfn file has only one spin. Use a reference .wfn "// &
267 0 : "with 2 spins separated or set REF_MO_SPIN to 1")
268 :
269 : ! Store only the mos required
270 20 : nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
271 20 : IF (proj_mo%ref_mo_index(1) == -1) THEN
272 6 : DEALLOCATE (proj_mo%ref_mo_index)
273 18 : ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
274 42 : DO i_ref = 1, nbr_mo_max
275 42 : proj_mo%ref_mo_index(i_ref) = i_ref
276 : END DO
277 : ELSE
278 38 : DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
279 24 : IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) &
280 : CALL cp_abort(__LOCATION__, &
281 : "The MO number available in the Reference SCF "// &
282 14 : "is smaller than the MO number you have required in REF_MO_INDEX.")
283 : END DO
284 : END IF
285 20 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
286 :
287 20 : IF (nbr_ref_mo > nbr_mo_max) &
288 : CALL cp_abort(__LOCATION__, &
289 0 : "The number of reference mo is larger then the total number of available one in the .wfn file.")
290 :
291 : ! Store
292 120 : ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
293 : CALL cp_fm_struct_create(mo_ref_fmstruct, &
294 : context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
295 : nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
296 20 : ncol_global=1)
297 20 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
298 :
299 80 : DO mo_index = 1, nbr_ref_mo
300 60 : real_mo_index = proj_mo%ref_mo_index(mo_index)
301 60 : IF (real_mo_index > nbr_mo_max) &
302 : CALL cp_abort(__LOCATION__, &
303 0 : "One of reference mo index is larger then the total number of available mo in the .wfn file.")
304 :
305 : ! fill with the reference mo values
306 60 : CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
307 : CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
308 : ncol=1, &
309 : source_start=real_mo_index, &
310 60 : target_start=1)
311 :
312 : ! multiply with overlap matrix to save time later on
313 80 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
314 : END DO
315 :
316 : ! Clean temporary variables
317 60 : DO ispin = 1, nspins
318 60 : CALL deallocate_mo_set(mo_ref_temp(ispin))
319 : END DO
320 20 : DEALLOCATE (mo_ref_temp)
321 :
322 20 : CALL cp_fm_struct_release(mo_ref_fmstruct)
323 20 : CALL cp_fm_release(mo_coeff_temp)
324 :
325 20 : END SUBROUTINE read_reference_mo_from_wfn
326 :
327 : ! **************************************************************************************************
328 : !> \brief Compute the projection of the current MO coefficients on reference ones
329 : !> and write the results.
330 : !> \param qs_env ...
331 : !> \param mos_new ...
332 : !> \param proj_mo ...
333 : !> \param n_proj ...
334 : !> \author Guillaume Le Breton
335 : ! **************************************************************************************************
336 40 : SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
337 : TYPE(qs_environment_type), POINTER :: qs_env
338 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
339 : TYPE(proj_mo_type) :: proj_mo
340 : INTEGER :: n_proj
341 :
342 : INTEGER :: i_ref, nbr_ref_mo, nbr_ref_td
343 40 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: phase, popu, sum_popu_ref
344 : TYPE(cp_logger_type), POINTER :: logger
345 : TYPE(dft_control_type), POINTER :: dft_control
346 : TYPE(section_vals_type), POINTER :: input, print_mo_section, proj_mo_section
347 :
348 40 : NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
349 :
350 40 : logger => cp_get_default_logger()
351 :
352 : CALL get_qs_env(qs_env, &
353 : dft_control=dft_control, &
354 40 : input=input)
355 :
356 : ! The general section
357 40 : proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
358 : ! The section we are dealing in this particular subroutine call: n_proj.
359 40 : print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
360 :
361 : ! STOP if not the required time step
362 40 : IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
363 : print_mo_section, ""), &
364 : cp_p_file)) &
365 : RETURN
366 :
367 38 : IF (.NOT. dft_control%rtp_control%fixed_ions) &
368 : CALL cp_abort(__LOCATION__, &
369 0 : "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO not implemented for EMD yet.")
370 :
371 38 : nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
372 38 : nbr_ref_td = SIZE(proj_mo%td_mo_index)
373 114 : ALLOCATE (popu(nbr_ref_td))
374 114 : ALLOCATE (phase(nbr_ref_td))
375 :
376 38 : IF (proj_mo%sum_on_all_ref) THEN
377 24 : ALLOCATE (sum_popu_ref(nbr_ref_td))
378 56 : sum_popu_ref(:) = 0.0_dp
379 56 : DO i_ref = 1, nbr_ref_mo
380 48 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
381 344 : sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
382 : END DO
383 8 : IF (proj_mo%sum_on_all_td) THEN
384 28 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=SUM(sum_popu_ref), n_proj=n_proj)
385 : ELSE
386 4 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
387 : END IF
388 8 : DEALLOCATE (sum_popu_ref)
389 : ELSE
390 98 : DO i_ref = 1, nbr_ref_mo
391 68 : CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
392 :
393 98 : IF (proj_mo%sum_on_all_td) THEN
394 196 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=SUM(popu), n_proj=n_proj)
395 : ELSE
396 :
397 40 : CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
398 : END IF
399 : END DO
400 : END IF
401 :
402 38 : DEALLOCATE (popu)
403 38 : DEALLOCATE (phase)
404 :
405 40 : END SUBROUTINE compute_and_write_proj_mo
406 :
407 : ! **************************************************************************************************
408 : !> \brief Compute the projection of the current MO coefficients on reference ones
409 : !> \param popu ...
410 : !> \param phase ...
411 : !> \param mos_new ...
412 : !> \param proj_mo ...
413 : !> \param i_ref ...
414 : !> \author Guillaume Le Breton
415 : ! **************************************************************************************************
416 116 : SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
417 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: popu, phase
418 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
419 : TYPE(proj_mo_type) :: proj_mo
420 : INTEGER :: i_ref
421 :
422 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_proj_mo'
423 :
424 : INTEGER :: handle, j_td, nbr_ref_td, spin_td
425 : REAL(KIND=dp) :: imag_proj, real_proj
426 : TYPE(cp_fm_struct_type), POINTER :: mo_ref_fmstruct
427 : TYPE(cp_fm_type) :: mo_coeff_temp
428 :
429 116 : CALL timeset(routineN, handle)
430 :
431 116 : nbr_ref_td = SIZE(popu)
432 116 : spin_td = proj_mo%td_mo_spin
433 :
434 : CALL cp_fm_struct_create(mo_ref_fmstruct, &
435 : context=mos_new(1)%matrix_struct%context, &
436 : nrow_global=mos_new(1)%matrix_struct%nrow_global, &
437 116 : ncol_global=1)
438 116 : CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
439 :
440 656 : DO j_td = 1, nbr_ref_td
441 : ! Real part of the projection:
442 : real_proj = 0.0_dp
443 : CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
444 : ncol=1, &
445 : source_start=proj_mo%td_mo_index(j_td), &
446 540 : target_start=1)
447 :
448 540 : CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
449 :
450 : ! Imaginary part of the projection
451 : imag_proj = 0.0_dp
452 : CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
453 : ncol=1, &
454 : source_start=proj_mo%td_mo_index(j_td), &
455 540 : target_start=1)
456 :
457 540 : CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
458 :
459 : ! Store the result
460 540 : phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians
461 1196 : popu(j_td) = real_proj**2 + imag_proj**2
462 :
463 : END DO
464 :
465 116 : CALL cp_fm_struct_release(mo_ref_fmstruct)
466 116 : CALL cp_fm_release(mo_coeff_temp)
467 :
468 116 : CALL timestop(handle)
469 :
470 116 : END SUBROUTINE compute_proj_mo
471 :
472 : ! **************************************************************************************************
473 : !> \brief Write in one file the projection of (all) the time-dependent MO coefficients
474 : !> on one reference ones
475 : !> \param qs_env ...
476 : !> \param print_mo_section ...
477 : !> \param proj_mo ...
478 : !> \param i_ref ...
479 : !> \param popu ...
480 : !> \param phase ...
481 : !> \param popu_tot ...
482 : !> \param n_proj ...
483 : !> \author Guillaume Le Breton
484 : ! **************************************************************************************************
485 76 : SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
486 : TYPE(qs_environment_type), POINTER :: qs_env
487 : TYPE(section_vals_type), POINTER :: print_mo_section
488 : TYPE(proj_mo_type) :: proj_mo
489 : INTEGER, OPTIONAL :: i_ref
490 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: popu, phase
491 : REAL(KIND=dp), OPTIONAL :: popu_tot
492 : INTEGER, OPTIONAL :: n_proj
493 :
494 : CHARACTER(LEN=default_string_length) :: ext, filename
495 : INTEGER :: j_td, output_unit, print_unit
496 : TYPE(cp_logger_type), POINTER :: logger
497 :
498 76 : NULLIFY (logger)
499 :
500 76 : logger => cp_get_default_logger()
501 76 : output_unit = cp_logger_get_default_io_unit(logger)
502 :
503 76 : IF (.NOT. (output_unit > 0)) RETURN
504 :
505 38 : IF (proj_mo%sum_on_all_ref) THEN
506 4 : ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))//"-ALL_REF.dat"
507 : ELSE
508 : ! Filename is update wrt the reference MO number
509 : ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))// &
510 : "-REF-"// &
511 : TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
512 34 : ".dat"
513 : END IF
514 :
515 : print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
516 38 : extension=TRIM(ext))
517 :
518 38 : IF (print_unit /= output_unit) THEN
519 38 : INQUIRE (UNIT=print_unit, NAME=filename)
520 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
521 38 : "PROJECTION MO", "The projection of the TD MOs is done in the file:", &
522 76 : TRIM(filename)
523 : WRITE (UNIT=print_unit, FMT="(/,(T2,A,T40,I6))") &
524 38 : "Real time propagation step:", qs_env%sim_step
525 : ELSE
526 0 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "PROJECTION MO"
527 : END IF
528 :
529 38 : IF (proj_mo%sum_on_all_ref) THEN
530 : WRITE (print_unit, "(T3,A)") &
531 : "Projection on all the required MO number from the reference file "// &
532 4 : TRIM(proj_mo%ref_mo_file_name)
533 4 : IF (proj_mo%sum_on_all_td) THEN
534 : WRITE (print_unit, "(T3, A, E20.12)") &
535 2 : "The sum over all the TD MOs population:", popu_tot
536 : ELSE
537 : WRITE (print_unit, "(T3,A)") &
538 2 : "For each TD MOs required is printed: Population "
539 14 : DO j_td = 1, SIZE(popu)
540 14 : WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
541 : END DO
542 : END IF
543 : ELSE
544 : WRITE (print_unit, "(T3,A)") &
545 : "Projection on the MO number "// &
546 : TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
547 : " from the reference file "// &
548 34 : TRIM(proj_mo%ref_mo_file_name)
549 :
550 34 : IF (proj_mo%sum_on_all_td) THEN
551 : WRITE (print_unit, "(T3, A, E20.12)") &
552 14 : "The sum over all the TD MOs population:", popu_tot
553 : ELSE
554 : WRITE (print_unit, "(T3,A)") &
555 20 : "For each TD MOs required is printed: Population & Phase [rad] "
556 62 : DO j_td = 1, SIZE(popu)
557 62 : WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
558 : END DO
559 : END IF
560 : END IF
561 :
562 38 : CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
563 :
564 : END SUBROUTINE write_proj_mo
565 :
566 : END MODULE rt_projection_mo_utils
567 :
|
