Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE qs_tddfpt2_restart
9 : USE cp_blacs_env, ONLY: cp_blacs_env_type
10 : USE cp_dbcsr_api, ONLY: dbcsr_type
11 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
12 : USE cp_files, ONLY: close_file,&
13 : open_file
14 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
15 : cp_fm_scale_and_add,&
16 : cp_fm_trace
17 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
18 : fm_pool_create_fm,&
19 : fm_pool_give_back_fm
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_get_info,&
25 : cp_fm_read_unformatted,&
26 : cp_fm_release,&
27 : cp_fm_type,&
28 : cp_fm_write_formatted,&
29 : cp_fm_write_info,&
30 : cp_fm_write_unformatted
31 : USE cp_log_handling, ONLY: cp_logger_type
32 : USE cp_output_handling, ONLY: cp_p_file,&
33 : cp_print_key_finished_output,&
34 : cp_print_key_generate_filename,&
35 : cp_print_key_should_output,&
36 : cp_print_key_unit_nr
37 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
38 : section_vals_type,&
39 : section_vals_val_get
40 : USE kinds, ONLY: default_path_length,&
41 : dp
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE parallel_gemm_api, ONLY: parallel_gemm
44 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
45 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
46 : USE string_utilities, ONLY: integer_to_string
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 :
51 : PRIVATE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_restart'
54 :
55 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
56 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
57 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
58 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
59 :
60 : PUBLIC :: tddfpt_write_restart, tddfpt_read_restart, tddfpt_write_newtonx_output, tddfpt_check_orthonormality
61 :
62 : ! **************************************************************************************************
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief Write Ritz vectors to a binary restart file.
68 : !> \param evects vectors to store
69 : !> \param evals TDDFPT eigenvalues
70 : !> \param gs_mos structure that holds ground state occupied and virtual
71 : !> molecular orbitals
72 : !> \param logger a logger object
73 : !> \param tddfpt_print_section TDDFPT%PRINT input section
74 : !> \par History
75 : !> * 08.2016 created [Sergey Chulkov]
76 : ! **************************************************************************************************
77 6556 : SUBROUTINE tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
78 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
79 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
80 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
81 : INTENT(in) :: gs_mos
82 : TYPE(cp_logger_type), POINTER :: logger
83 : TYPE(section_vals_type), POINTER :: tddfpt_print_section
84 :
85 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_restart'
86 :
87 : INTEGER :: handle, ispin, istate, nao, nspins, &
88 : nstates, ounit
89 : INTEGER, DIMENSION(maxspins) :: nmo_active
90 :
91 6556 : IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "RESTART"), cp_p_file)) THEN
92 1160 : CALL timeset(routineN, handle)
93 :
94 1160 : nspins = SIZE(evects, 1)
95 1160 : nstates = SIZE(evects, 2)
96 :
97 : IF (debug_this_module) THEN
98 : CPASSERT(SIZE(evals) == nstates)
99 : CPASSERT(nspins > 0)
100 : CPASSERT(nstates > 0)
101 : END IF
102 :
103 1160 : CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
104 2450 : DO ispin = 1, nspins
105 2450 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nmo_active(ispin))
106 : END DO
107 :
108 : ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "RESTART", &
109 : extension=".tdwfn", file_status="REPLACE", file_action="WRITE", &
110 1160 : do_backup=.TRUE., file_form="UNFORMATTED")
111 :
112 1160 : IF (ounit > 0) THEN
113 580 : WRITE (ounit) nstates, nspins, nao
114 580 : WRITE (ounit) nmo_active(1:nspins)
115 580 : WRITE (ounit) evals
116 : END IF
117 :
118 4294 : DO istate = 1, nstates
119 7832 : DO ispin = 1, nspins
120 : ! TDDFPT wave function is actually stored as a linear combination of virtual MOs
121 : ! that replaces the corresponding deoccupied MO. Unfortunately, the phase
122 : ! of the occupied MOs varies depending on the eigensolver used as well as
123 : ! how eigenvectors are distributed across computational cores. The phase is important
124 : ! because TDDFPT wave functions are used to compute a response electron density
125 : ! \rho^{-} = 1/2 * [C_{0} * evect^T + evect * C_{0}^{-}], where C_{0} is the expansion
126 : ! coefficients of the reference ground-state wave function. To make the restart file
127 : ! transferable, TDDFPT wave functions are stored in assumption that all ground state
128 : ! MOs have a positive phase.
129 3538 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
130 :
131 3538 : CALL cp_fm_write_unformatted(evects(ispin, istate), ounit)
132 :
133 6672 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
134 : END DO
135 : END DO
136 :
137 1160 : CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "RESTART")
138 :
139 1160 : CALL timestop(handle)
140 : END IF
141 :
142 6556 : END SUBROUTINE tddfpt_write_restart
143 :
144 : ! **************************************************************************************************
145 : !> \brief Initialise initial guess vectors by reading (un-normalised) Ritz vectors
146 : !> from a binary restart file.
147 : !> \param evects vectors to initialise (initialised on exit)
148 : !> \param evals TDDFPT eigenvalues (initialised on exit)
149 : !> \param gs_mos structure that holds ground state occupied and virtual
150 : !> molecular orbitals
151 : !> \param logger a logger object
152 : !> \param tddfpt_section TDDFPT input section
153 : !> \param tddfpt_print_section TDDFPT%PRINT input section
154 : !> \param fm_pool_ao_mo_active pools of dense matrices with shape [nao x nmo_active(spin)]
155 : !> \param blacs_env_global BLACS parallel environment involving all the processor
156 : !> \return the number of excited states found in the restart file
157 : !> \par History
158 : !> * 08.2016 created [Sergey Chulkov]
159 : ! **************************************************************************************************
160 12 : FUNCTION tddfpt_read_restart(evects, evals, gs_mos, logger, tddfpt_section, tddfpt_print_section, &
161 6 : fm_pool_ao_mo_active, blacs_env_global) RESULT(nstates_read)
162 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
163 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals
164 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
165 : INTENT(in) :: gs_mos
166 : TYPE(cp_logger_type), POINTER :: logger
167 : TYPE(section_vals_type), POINTER :: tddfpt_section, tddfpt_print_section
168 : TYPE(cp_fm_pool_p_type), DIMENSION(:), INTENT(in) :: fm_pool_ao_mo_active
169 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
170 : INTEGER :: nstates_read
171 :
172 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_read_restart'
173 :
174 : CHARACTER(len=20) :: read_str, ref_str
175 : CHARACTER(LEN=default_path_length) :: filename
176 : INTEGER :: handle, ispin, istate, iunit, n_rep_val, &
177 : nao, nao_read, nspins, nspins_read, &
178 : nstates
179 : INTEGER, DIMENSION(maxspins) :: nmo_active, nmo_active_read
180 : LOGICAL :: file_exists
181 6 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_read
182 : TYPE(cp_fm_type) :: evtest
183 : TYPE(mp_para_env_type), POINTER :: para_env_global
184 : TYPE(section_vals_type), POINTER :: print_key
185 :
186 6 : CALL timeset(routineN, handle)
187 :
188 6 : CPASSERT(ASSOCIATED(tddfpt_section))
189 :
190 : ! generate restart file name
191 6 : CALL section_vals_val_get(tddfpt_section, "WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
192 6 : IF (n_rep_val > 0) THEN
193 0 : CALL section_vals_val_get(tddfpt_section, "WFN_RESTART_FILE_NAME", c_val=filename)
194 : ELSE
195 6 : print_key => section_vals_get_subs_vals(tddfpt_print_section, "RESTART")
196 : filename = cp_print_key_generate_filename(logger, print_key, &
197 6 : extension=".tdwfn", my_local=.FALSE.)
198 : END IF
199 :
200 6 : CALL blacs_env_global%get(para_env=para_env_global)
201 :
202 6 : IF (para_env_global%is_source()) THEN
203 3 : INQUIRE (FILE=filename, exist=file_exists)
204 :
205 3 : IF (.NOT. file_exists) THEN
206 2 : nstates_read = 0
207 2 : CALL para_env_global%bcast(nstates_read)
208 :
209 : CALL cp_warn(__LOCATION__, &
210 : "User requested to restart the TDDFPT wave functions from the file '"//TRIM(filename)// &
211 2 : "' which does not exist. Guess wave functions will be constructed using Kohn-Sham orbitals.")
212 2 : CALL timestop(handle)
213 2 : RETURN
214 : END IF
215 :
216 : CALL open_file(file_name=filename, file_action="READ", file_form="UNFORMATTED", &
217 1 : file_status="OLD", unit_number=iunit)
218 : END IF
219 :
220 4 : nspins = SIZE(evects, 1)
221 4 : nstates = SIZE(evects, 2)
222 :
223 12 : DO ispin = 1, nspins
224 8 : CALL fm_pool_create_fm(fm_pool_ao_mo_active(ispin)%pool, evtest)
225 8 : CALL cp_fm_get_info(evtest, nrow_global=nao, ncol_global=nmo_active(ispin))
226 12 : CALL fm_pool_give_back_fm(fm_pool_ao_mo_active(ispin)%pool, evtest)
227 : END DO
228 :
229 4 : IF (para_env_global%is_source()) THEN
230 1 : READ (iunit) nstates_read, nspins_read, nao_read
231 :
232 1 : IF (nspins_read /= nspins) THEN
233 0 : CALL integer_to_string(nspins, ref_str)
234 0 : CALL integer_to_string(nspins_read, read_str)
235 : CALL cp_abort(__LOCATION__, &
236 : "Restarted TDDFPT wave function contains incompatible number of spin components ("// &
237 0 : TRIM(read_str)//" instead of "//TRIM(ref_str)//").")
238 : END IF
239 :
240 1 : IF (nao_read /= nao) THEN
241 0 : CALL integer_to_string(nao, ref_str)
242 0 : CALL integer_to_string(nao_read, read_str)
243 : CALL cp_abort(__LOCATION__, &
244 0 : "Incompatible number of atomic orbitals ("//TRIM(read_str)//" instead of "//TRIM(ref_str)//").")
245 : END IF
246 :
247 1 : READ (iunit) nmo_active_read(1:nspins)
248 :
249 3 : DO ispin = 1, nspins
250 3 : IF (nmo_active_read(ispin) /= nmo_active(ispin)) THEN
251 : CALL cp_abort(__LOCATION__, &
252 0 : "Incompatible number of electrons and/or multiplicity.")
253 : END IF
254 : END DO
255 :
256 1 : IF (nstates_read /= nstates) THEN
257 0 : CALL integer_to_string(nstates, ref_str)
258 0 : CALL integer_to_string(nstates_read, read_str)
259 : CALL cp_warn(__LOCATION__, &
260 : "TDDFPT restart file contains "//TRIM(read_str)// &
261 : " wave function(s) however "//TRIM(ref_str)// &
262 0 : " excited states were requested.")
263 : END IF
264 : END IF
265 4 : CALL para_env_global%bcast(nstates_read)
266 :
267 : ! exit if restart file does not exist
268 4 : IF (nstates_read <= 0) THEN
269 2 : CALL timestop(handle)
270 2 : RETURN
271 : END IF
272 :
273 2 : IF (para_env_global%is_source()) THEN
274 3 : ALLOCATE (evals_read(nstates_read))
275 1 : READ (iunit) evals_read
276 1 : IF (nstates_read <= nstates) THEN
277 4 : evals(1:nstates_read) = evals_read(1:nstates_read)
278 : ELSE
279 0 : evals(1:nstates) = evals_read(1:nstates)
280 : END IF
281 1 : DEALLOCATE (evals_read)
282 : END IF
283 14 : CALL para_env_global%bcast(evals)
284 :
285 8 : DO istate = 1, nstates_read
286 20 : DO ispin = 1, nspins
287 18 : IF (istate <= nstates) THEN
288 12 : CALL fm_pool_create_fm(fm_pool_ao_mo_active(ispin)%pool, evects(ispin, istate))
289 :
290 12 : CALL cp_fm_read_unformatted(evects(ispin, istate), iunit)
291 :
292 12 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
293 : END IF
294 : END DO
295 : END DO
296 :
297 2 : IF (para_env_global%is_source()) &
298 1 : CALL close_file(unit_number=iunit)
299 :
300 2 : CALL timestop(handle)
301 :
302 12 : END FUNCTION tddfpt_read_restart
303 : ! **************************************************************************************************
304 : !> \brief Write Ritz vectors to a binary restart file.
305 : !> \param evects vectors to store
306 : !> \param evals TDDFPT eigenvalues
307 : !> \param gs_mos structure that holds ground state occupied and virtual
308 : !> molecular orbitals
309 : !> \param logger a logger object
310 : !> \param tddfpt_print_section TDDFPT%PRINT input section
311 : !> \param matrix_s ...
312 : !> \param S_evects ...
313 : !> \param sub_env ...
314 : ! **************************************************************************************************
315 2 : SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, &
316 2 : matrix_s, S_evects, sub_env)
317 :
318 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
319 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
320 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
321 : INTENT(in) :: gs_mos
322 : TYPE(cp_logger_type), INTENT(in), POINTER :: logger
323 : TYPE(section_vals_type), INTENT(in), POINTER :: tddfpt_print_section
324 : TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
325 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: S_evects
326 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
327 :
328 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_newtonx_output'
329 :
330 : INTEGER :: handle, iocc, ispin, istate, ivirt, nao, &
331 : nspins, nstates, ounit
332 : INTEGER, DIMENSION(maxspins) :: nmo_active, nmo_occ, nmo_virt
333 : LOGICAL :: print_phases, print_virtuals, &
334 : scale_with_phases
335 2 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase_evects
336 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
337 2 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects_mo
338 :
339 2 : IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "NAMD_PRINT"), cp_p_file)) THEN
340 2 : CALL timeset(routineN, handle)
341 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals)
342 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
343 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%SCALE_WITH_PHASES", l_val=scale_with_phases)
344 :
345 2 : nspins = SIZE(evects, 1)
346 2 : nstates = SIZE(evects, 2)
347 :
348 : IF (debug_this_module) THEN
349 : CPASSERT(SIZE(evals) == nstates)
350 : CPASSERT(nspins > 0)
351 : CPASSERT(nstates > 0)
352 : END IF
353 :
354 2 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
355 :
356 2 : IF (sub_env%is_split) THEN
357 : CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when states"// &
358 0 : " are distributed to different CPU pools.")
359 : END IF
360 :
361 : ! test for reduced active orbitals
362 4 : DO ispin = 1, nspins
363 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
364 2 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nmo_active(ispin))
365 4 : IF (nmo_occ(ispin) /= nmo_active(ispin)) THEN
366 : CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when using"// &
367 0 : " a reduced set of active occupied orbitals.")
368 : END IF
369 : END DO
370 :
371 : ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "NAMD_PRINT", &
372 2 : extension=".inp", file_form="FORMATTED", file_action="WRITE", file_status="REPLACE")
373 : IF (debug_this_module) CALL tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
374 :
375 : ! print eigenvectors
376 2 : IF (print_virtuals) THEN
377 12 : ALLOCATE (evects_mo(nspins, nstates))
378 4 : DO istate = 1, nstates
379 6 : DO ispin = 1, nspins
380 :
381 : ! transform eigenvectors
382 2 : NULLIFY (fmstruct)
383 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
384 2 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
385 : CALL cp_fm_struct_create(fmstruct, para_env=sub_env%para_env, &
386 : context=sub_env%blacs_env, &
387 2 : nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin))
388 2 : CALL cp_fm_create(evects_mo(ispin, istate), fmstruct)
389 2 : CALL cp_fm_struct_release(fmstruct)
390 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, istate), S_evects(ispin, istate), &
391 4 : ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
392 : END DO
393 : END DO
394 4 : DO istate = 1, nstates
395 6 : DO ispin = 1, nspins
396 : CALL parallel_gemm("T", "N", &
397 : nmo_virt(ispin), &
398 : nmo_occ(ispin), &
399 : nao, &
400 : 1.0_dp, &
401 : gs_mos(ispin)%mos_virt, &
402 : S_evects(ispin, istate), & !this also needs to be orthogonalized
403 : 0.0_dp, &
404 4 : evects_mo(ispin, istate))
405 : END DO
406 : END DO
407 : END IF
408 :
409 4 : DO istate = 1, nstates
410 6 : DO ispin = 1, nspins
411 :
412 2 : IF (.NOT. print_virtuals) THEN
413 0 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
414 0 : IF (ounit > 0) THEN
415 0 : WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
416 0 : CALL cp_fm_write_info(evects(ispin, istate), ounit)
417 : END IF
418 0 : CALL cp_fm_write_formatted(evects(ispin, istate), ounit, "ES EIGENVECTORS")
419 : ELSE
420 2 : CALL cp_fm_column_scale(evects_mo(ispin, istate), gs_mos(ispin)%phases_occ)
421 2 : IF (ounit > 0) THEN
422 1 : WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
423 1 : CALL cp_fm_write_info(evects_mo(ispin, istate), ounit)
424 : END IF
425 2 : CALL cp_fm_write_formatted(evects_mo(ispin, istate), ounit, "ES EIGENVECTORS")
426 : END IF
427 :
428 : ! compute and print phase of eigenvectors
429 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
430 6 : ALLOCATE (phase_evects(nmo_occ(ispin)))
431 2 : IF (print_virtuals) THEN
432 2 : CALL compute_phase_eigenvectors(evects_mo(ispin, istate), phase_evects, sub_env)
433 : ELSE
434 0 : CALL compute_phase_eigenvectors(evects(ispin, istate), phase_evects, sub_env)
435 : END IF
436 2 : IF (ounit > 0) THEN
437 1 : WRITE (ounit, "(/,A,/)") "PHASES ES EIGENVECTORS"
438 5 : DO iocc = 1, nmo_occ(ispin)
439 5 : WRITE (ounit, "(F20.14)") phase_evects(iocc)
440 : END DO
441 : END IF
442 4 : DEALLOCATE (phase_evects)
443 :
444 : END DO
445 : END DO
446 :
447 2 : IF (print_virtuals) THEN
448 2 : CALL cp_fm_release(evects_mo)
449 : END IF
450 :
451 4 : DO ispin = 1, nspins
452 2 : IF (ounit > 0) THEN
453 1 : WRITE (ounit, "(/,A)") "OCCUPIED MOS SIZE"
454 1 : CALL cp_fm_write_info(gs_mos(ispin)%mos_occ, ounit)
455 : END IF
456 4 : CALL cp_fm_write_formatted(gs_mos(ispin)%mos_occ, ounit, "OCCUPIED MO COEFFICIENTS")
457 : END DO
458 :
459 2 : IF (ounit > 0) THEN
460 1 : WRITE (ounit, "(A)") "OCCUPIED MO EIGENVALUES"
461 2 : DO ispin = 1, nspins
462 1 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
463 6 : DO iocc = 1, nmo_occ(ispin)
464 5 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_occ(iocc)
465 : END DO
466 : END DO
467 : END IF
468 : !
469 2 : IF (print_virtuals) THEN
470 4 : DO ispin = 1, nspins
471 2 : IF (ounit > 0) THEN
472 1 : WRITE (ounit, "(/,A)") "VIRTUAL MOS SIZE"
473 1 : CALL cp_fm_write_info(gs_mos(ispin)%mos_virt, ounit)
474 : END IF
475 4 : CALL cp_fm_write_formatted(gs_mos(ispin)%mos_virt, ounit, "VIRTUAL MO COEFFICIENTS")
476 : END DO
477 :
478 2 : IF (ounit > 0) THEN
479 1 : WRITE (ounit, "(A)") "VIRTUAL MO EIGENVALUES"
480 2 : DO ispin = 1, nspins
481 1 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
482 21 : DO ivirt = 1, nmo_virt(ispin)
483 20 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_virt(ivirt)
484 : END DO
485 : END DO
486 : END IF
487 : END IF
488 :
489 : ! print phases of molecular orbitals
490 :
491 2 : IF (print_phases) THEN
492 0 : IF (ounit > 0) THEN
493 0 : WRITE (ounit, "(A)") "PHASES OCCUPIED ORBITALS"
494 0 : DO ispin = 1, nspins
495 0 : DO iocc = 1, nmo_occ(ispin)
496 0 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_occ(iocc)
497 : END DO
498 : END DO
499 0 : IF (print_virtuals) THEN
500 0 : WRITE (ounit, "(A)") "PHASES VIRTUAL ORBITALS"
501 0 : DO ispin = 1, nspins
502 0 : DO ivirt = 1, nmo_virt(ispin)
503 0 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_virt(ivirt)
504 : END DO
505 : END DO
506 : END IF
507 : END IF
508 : END IF
509 :
510 2 : CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "NAMD_PRINT")
511 :
512 2 : CALL timestop(handle)
513 : END IF
514 :
515 2 : END SUBROUTINE tddfpt_write_newtonx_output
516 : ! **************************************************************************************************
517 : !> \brief ...
518 : !> \param evects ...
519 : !> \param ounit ...
520 : !> \param S_evects ...
521 : !> \param matrix_s ...
522 : ! **************************************************************************************************
523 0 : SUBROUTINE tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
524 :
525 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
526 : INTEGER, INTENT(in) :: ounit
527 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: S_evects
528 : TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
529 :
530 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_check_orthonormality'
531 :
532 : INTEGER :: handle, ispin, ivect, jvect, nspins, &
533 : nvects_total
534 : INTEGER, DIMENSION(maxspins) :: nactive
535 : REAL(kind=dp) :: norm
536 : REAL(kind=dp), DIMENSION(maxspins) :: weights
537 :
538 0 : CALL timeset(routineN, handle)
539 :
540 0 : nspins = SIZE(evects, 1)
541 0 : nvects_total = SIZE(evects, 2)
542 :
543 : IF (debug_this_module) THEN
544 : CPASSERT(SIZE(S_evects, 1) == nspins)
545 : CPASSERT(SIZE(S_evects, 2) == nvects_total)
546 : END IF
547 :
548 0 : DO ispin = 1, nspins
549 0 : CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
550 : END DO
551 :
552 0 : DO jvect = 1, nvects_total
553 : ! <psi1_i | psi1_j>
554 0 : DO ivect = 1, jvect - 1
555 0 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
556 0 : norm = SUM(weights(1:nspins))
557 :
558 0 : DO ispin = 1, nspins
559 0 : CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect), -norm, evects(ispin, ivect))
560 : END DO
561 : END DO
562 :
563 : ! <psi1_j | psi1_j>
564 0 : DO ispin = 1, nspins
565 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
566 0 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
567 : END DO
568 :
569 0 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
570 :
571 0 : norm = SUM(weights(1:nspins))
572 : norm = 1.0_dp/SQRT(norm)
573 :
574 0 : IF ((ounit > 0) .AND. debug_this_module) WRITE (ounit, '(A,F10.8)') "norm", norm
575 :
576 : END DO
577 :
578 0 : CALL timestop(handle)
579 :
580 0 : END SUBROUTINE tddfpt_check_orthonormality
581 : ! **************************************************************************************************
582 : !> \brief ...
583 : !> \param evects ...
584 : !> \param phase_evects ...
585 : !> \param sub_env ...
586 : ! **************************************************************************************************
587 2 : SUBROUTINE compute_phase_eigenvectors(evects, phase_evects, sub_env)
588 :
589 : ! copied from parts of tddgpt_init_ground_state_mos by S. Chulkov
590 :
591 : TYPE(cp_fm_type), INTENT(in) :: evects
592 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: phase_evects
593 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
594 :
595 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_phase_eigenvectors'
596 : REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp)
597 :
598 : INTEGER :: handle, icol_global, icol_local, irow_global, irow_local, ncol_global, &
599 : ncol_local, nrow_global, nrow_local, sign_int
600 : INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
601 : sum_sign_array
602 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
603 : REAL(kind=dp) :: element
604 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
605 2 : POINTER :: my_block
606 :
607 2 : CALL timeset(routineN, handle)
608 :
609 : ! compute and print the phase of excited-state eigenvectors:
610 : CALL cp_fm_get_info(evects, nrow_global=nrow_global, ncol_global=ncol_global, &
611 : nrow_local=nrow_local, ncol_local=ncol_local, local_data=my_block, &
612 2 : row_indices=row_indices, col_indices=col_indices) ! nrow_global either nao or nocc
613 :
614 14 : ALLOCATE (minrow_neg_array(ncol_global), minrow_pos_array(ncol_global), sum_sign_array(ncol_global))
615 10 : minrow_neg_array(:) = nrow_global
616 10 : minrow_pos_array(:) = nrow_global
617 10 : sum_sign_array(:) = 0
618 :
619 10 : DO icol_local = 1, ncol_local
620 8 : icol_global = col_indices(icol_local)
621 :
622 86 : DO irow_local = 1, nrow_local
623 76 : irow_global = row_indices(irow_local)
624 :
625 76 : element = my_block(irow_local, icol_local)
626 :
627 76 : sign_int = 0
628 76 : IF (element >= eps_dp) THEN
629 : sign_int = 1
630 36 : ELSE IF (element <= -eps_dp) THEN
631 36 : sign_int = -1
632 : END IF
633 :
634 76 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
635 :
636 84 : IF (sign_int > 0) THEN
637 40 : IF (minrow_pos_array(icol_global) > irow_global) &
638 8 : minrow_pos_array(icol_global) = irow_global
639 36 : ELSE IF (sign_int < 0) THEN
640 36 : IF (minrow_neg_array(icol_global) > irow_global) &
641 8 : minrow_neg_array(icol_global) = irow_global
642 : END IF
643 :
644 : END DO
645 : END DO
646 :
647 2 : CALL sub_env%para_env%sum(sum_sign_array)
648 2 : CALL sub_env%para_env%min(minrow_neg_array)
649 2 : CALL sub_env%para_env%min(minrow_pos_array)
650 :
651 10 : DO icol_global = 1, ncol_global
652 :
653 10 : IF (sum_sign_array(icol_global) > 0) THEN
654 : ! most of the expansion coefficients are positive => MO's phase = +1
655 6 : phase_evects(icol_global) = 1.0_dp
656 2 : ELSE IF (sum_sign_array(icol_global) < 0) THEN
657 : ! most of the expansion coefficients are negative => MO's phase = -1
658 2 : phase_evects(icol_global) = -1.0_dp
659 : ELSE
660 : ! equal number of positive and negative expansion coefficients
661 0 : IF (minrow_pos_array(icol_global) <= minrow_neg_array(icol_global)) THEN
662 : ! the first positive expansion coefficient has a lower index then
663 : ! the first negative expansion coefficient; MO's phase = +1
664 0 : phase_evects(icol_global) = 1.0_dp
665 : ELSE
666 : ! MO's phase = -1
667 0 : phase_evects(icol_global) = -1.0_dp
668 : END IF
669 : END IF
670 :
671 : END DO
672 :
673 2 : DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
674 :
675 2 : CALL timestop(handle)
676 :
677 2 : END SUBROUTINE compute_phase_eigenvectors
678 :
679 : END MODULE qs_tddfpt2_restart
|
