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 : 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 8260 : 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 8260 : IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "RESTART"), cp_p_file)) THEN
92 1468 : CALL timeset(routineN, handle)
93 :
94 1468 : nspins = SIZE(evects, 1)
95 1468 : 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 1468 : CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
104 3112 : DO ispin = 1, nspins
105 3112 : 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 1468 : do_backup=.TRUE., file_form="UNFORMATTED")
111 :
112 1468 : IF (ounit > 0) THEN
113 734 : WRITE (ounit) nstates, nspins, nao
114 734 : WRITE (ounit) nmo_active(1:nspins)
115 734 : WRITE (ounit) evals
116 : END IF
117 :
118 5222 : DO istate = 1, nstates
119 9600 : 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 4378 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
130 :
131 4378 : CALL cp_fm_write_unformatted(evects(ispin, istate), ounit)
132 :
133 8132 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
134 : END DO
135 : END DO
136 :
137 1468 : CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "RESTART")
138 :
139 1468 : CALL timestop(handle)
140 : END IF
141 :
142 8260 : 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()) THEN
298 1 : CALL close_file(unit_number=iunit)
299 : END IF
300 :
301 2 : CALL timestop(handle)
302 :
303 12 : END FUNCTION tddfpt_read_restart
304 : ! **************************************************************************************************
305 : !> \brief Write Ritz vectors to a binary restart file.
306 : !> \param evects vectors to store
307 : !> \param evals TDDFPT eigenvalues
308 : !> \param gs_mos structure that holds ground state occupied and virtual
309 : !> molecular orbitals
310 : !> \param logger a logger object
311 : !> \param tddfpt_print_section TDDFPT%PRINT input section
312 : !> \param matrix_s ...
313 : !> \param S_evects ...
314 : !> \param sub_env ...
315 : ! **************************************************************************************************
316 2 : SUBROUTINE tddfpt_write_newtonx_output(evects, evals, gs_mos, logger, tddfpt_print_section, &
317 2 : matrix_s, S_evects, sub_env)
318 :
319 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
320 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
321 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
322 : INTENT(in) :: gs_mos
323 : TYPE(cp_logger_type), INTENT(in), POINTER :: logger
324 : TYPE(section_vals_type), INTENT(in), POINTER :: tddfpt_print_section
325 : TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
326 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: S_evects
327 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
328 :
329 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_write_newtonx_output'
330 :
331 : INTEGER :: handle, iocc, ispin, istate, ivirt, nao, &
332 : nspins, nstates, ounit
333 : INTEGER, DIMENSION(maxspins) :: nmo_active, nmo_occ, nmo_virt
334 : LOGICAL :: print_phases, print_virtuals, &
335 : scale_with_phases
336 2 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: phase_evects
337 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
338 2 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: evects_mo
339 :
340 2 : IF (BTEST(cp_print_key_should_output(logger%iter_info, tddfpt_print_section, "NAMD_PRINT"), cp_p_file)) THEN
341 2 : CALL timeset(routineN, handle)
342 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_VIRTUALS", l_val=print_virtuals)
343 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%PRINT_PHASES", l_val=print_phases)
344 2 : CALL section_vals_val_get(tddfpt_print_section, "NAMD_PRINT%SCALE_WITH_PHASES", l_val=scale_with_phases)
345 :
346 2 : nspins = SIZE(evects, 1)
347 2 : nstates = SIZE(evects, 2)
348 :
349 : IF (debug_this_module) THEN
350 : CPASSERT(SIZE(evals) == nstates)
351 : CPASSERT(nspins > 0)
352 : CPASSERT(nstates > 0)
353 : END IF
354 :
355 2 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
356 :
357 2 : IF (sub_env%is_split) THEN
358 : CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when states"// &
359 0 : " are distributed to different CPU pools.")
360 : END IF
361 :
362 : ! test for reduced active orbitals
363 4 : DO ispin = 1, nspins
364 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
365 2 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nmo_active(ispin))
366 4 : IF (nmo_occ(ispin) /= nmo_active(ispin)) THEN
367 : CALL cp_abort(__LOCATION__, "NEWTONX interface print not possible when using"// &
368 0 : " a reduced set of active occupied orbitals.")
369 : END IF
370 : END DO
371 :
372 : ounit = cp_print_key_unit_nr(logger, tddfpt_print_section, "NAMD_PRINT", &
373 2 : extension=".inp", file_form="FORMATTED", file_action="WRITE", file_status="REPLACE")
374 : IF (debug_this_module) CALL tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
375 :
376 : ! print eigenvectors
377 2 : IF (print_virtuals) THEN
378 12 : ALLOCATE (evects_mo(nspins, nstates))
379 4 : DO istate = 1, nstates
380 6 : DO ispin = 1, nspins
381 :
382 : ! transform eigenvectors
383 2 : NULLIFY (fmstruct)
384 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
385 2 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
386 : CALL cp_fm_struct_create(fmstruct, para_env=sub_env%para_env, &
387 : context=sub_env%blacs_env, &
388 2 : nrow_global=nmo_virt(ispin), ncol_global=nmo_occ(ispin))
389 2 : CALL cp_fm_create(evects_mo(ispin, istate), fmstruct)
390 2 : CALL cp_fm_struct_release(fmstruct)
391 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, istate), S_evects(ispin, istate), &
392 4 : ncol=nmo_occ(ispin), alpha=1.0_dp, beta=0.0_dp)
393 : END DO
394 : END DO
395 4 : DO istate = 1, nstates
396 6 : DO ispin = 1, nspins
397 : CALL parallel_gemm("T", "N", &
398 : nmo_virt(ispin), &
399 : nmo_occ(ispin), &
400 : nao, &
401 : 1.0_dp, &
402 : gs_mos(ispin)%mos_virt, &
403 : S_evects(ispin, istate), & !this also needs to be orthogonalized
404 : 0.0_dp, &
405 4 : evects_mo(ispin, istate))
406 : END DO
407 : END DO
408 : END IF
409 :
410 4 : DO istate = 1, nstates
411 6 : DO ispin = 1, nspins
412 :
413 2 : IF (.NOT. print_virtuals) THEN
414 0 : CALL cp_fm_column_scale(evects(ispin, istate), gs_mos(ispin)%phases_occ)
415 0 : IF (ounit > 0) THEN
416 0 : WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
417 0 : CALL cp_fm_write_info(evects(ispin, istate), ounit)
418 : END IF
419 0 : CALL cp_fm_write_formatted(evects(ispin, istate), ounit, "ES EIGENVECTORS")
420 : ELSE
421 2 : CALL cp_fm_column_scale(evects_mo(ispin, istate), gs_mos(ispin)%phases_occ)
422 2 : IF (ounit > 0) THEN
423 1 : WRITE (ounit, "(/,A)") "ES EIGENVECTORS SIZE"
424 1 : CALL cp_fm_write_info(evects_mo(ispin, istate), ounit)
425 : END IF
426 2 : CALL cp_fm_write_formatted(evects_mo(ispin, istate), ounit, "ES EIGENVECTORS")
427 : END IF
428 :
429 : ! compute and print phase of eigenvectors
430 2 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
431 6 : ALLOCATE (phase_evects(nmo_occ(ispin)))
432 2 : IF (print_virtuals) THEN
433 2 : CALL compute_phase_eigenvectors(evects_mo(ispin, istate), phase_evects, sub_env)
434 : ELSE
435 0 : CALL compute_phase_eigenvectors(evects(ispin, istate), phase_evects, sub_env)
436 : END IF
437 2 : IF (ounit > 0) THEN
438 1 : WRITE (ounit, "(/,A,/)") "PHASES ES EIGENVECTORS"
439 5 : DO iocc = 1, nmo_occ(ispin)
440 5 : WRITE (ounit, "(F20.14)") phase_evects(iocc)
441 : END DO
442 : END IF
443 4 : DEALLOCATE (phase_evects)
444 :
445 : END DO
446 : END DO
447 :
448 2 : IF (print_virtuals) THEN
449 2 : CALL cp_fm_release(evects_mo)
450 : END IF
451 :
452 4 : DO ispin = 1, nspins
453 2 : IF (ounit > 0) THEN
454 1 : WRITE (ounit, "(/,A)") "OCCUPIED MOS SIZE"
455 1 : CALL cp_fm_write_info(gs_mos(ispin)%mos_occ, ounit)
456 : END IF
457 4 : CALL cp_fm_write_formatted(gs_mos(ispin)%mos_occ, ounit, "OCCUPIED MO COEFFICIENTS")
458 : END DO
459 :
460 2 : IF (ounit > 0) THEN
461 1 : WRITE (ounit, "(A)") "OCCUPIED MO EIGENVALUES"
462 2 : DO ispin = 1, nspins
463 1 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
464 6 : DO iocc = 1, nmo_occ(ispin)
465 5 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_occ(iocc)
466 : END DO
467 : END DO
468 : END IF
469 : !
470 2 : IF (print_virtuals) THEN
471 4 : DO ispin = 1, nspins
472 2 : IF (ounit > 0) THEN
473 1 : WRITE (ounit, "(/,A)") "VIRTUAL MOS SIZE"
474 1 : CALL cp_fm_write_info(gs_mos(ispin)%mos_virt, ounit)
475 : END IF
476 4 : CALL cp_fm_write_formatted(gs_mos(ispin)%mos_virt, ounit, "VIRTUAL MO COEFFICIENTS")
477 : END DO
478 :
479 2 : IF (ounit > 0) THEN
480 1 : WRITE (ounit, "(A)") "VIRTUAL MO EIGENVALUES"
481 2 : DO ispin = 1, nspins
482 1 : nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
483 21 : DO ivirt = 1, nmo_virt(ispin)
484 20 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%evals_virt(ivirt)
485 : END DO
486 : END DO
487 : END IF
488 : END IF
489 :
490 : ! print phases of molecular orbitals
491 :
492 2 : IF (print_phases) THEN
493 0 : IF (ounit > 0) THEN
494 0 : WRITE (ounit, "(A)") "PHASES OCCUPIED ORBITALS"
495 0 : DO ispin = 1, nspins
496 0 : DO iocc = 1, nmo_occ(ispin)
497 0 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_occ(iocc)
498 : END DO
499 : END DO
500 0 : IF (print_virtuals) THEN
501 0 : WRITE (ounit, "(A)") "PHASES VIRTUAL ORBITALS"
502 0 : DO ispin = 1, nspins
503 0 : DO ivirt = 1, nmo_virt(ispin)
504 0 : WRITE (ounit, "(F20.14)") gs_mos(ispin)%phases_virt(ivirt)
505 : END DO
506 : END DO
507 : END IF
508 : END IF
509 : END IF
510 :
511 2 : CALL cp_print_key_finished_output(ounit, logger, tddfpt_print_section, "NAMD_PRINT")
512 :
513 2 : CALL timestop(handle)
514 : END IF
515 :
516 2 : END SUBROUTINE tddfpt_write_newtonx_output
517 : ! **************************************************************************************************
518 : !> \brief ...
519 : !> \param evects ...
520 : !> \param ounit ...
521 : !> \param S_evects ...
522 : !> \param matrix_s ...
523 : ! **************************************************************************************************
524 0 : SUBROUTINE tddfpt_check_orthonormality(evects, ounit, S_evects, matrix_s)
525 :
526 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
527 : INTEGER, INTENT(in) :: ounit
528 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: S_evects
529 : TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
530 :
531 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_check_orthonormality'
532 :
533 : INTEGER :: handle, ispin, ivect, jvect, nspins, &
534 : nvects_total
535 : INTEGER, DIMENSION(maxspins) :: nactive
536 : REAL(kind=dp) :: norm
537 : REAL(kind=dp), DIMENSION(maxspins) :: weights
538 :
539 0 : CALL timeset(routineN, handle)
540 :
541 0 : nspins = SIZE(evects, 1)
542 0 : nvects_total = SIZE(evects, 2)
543 :
544 : IF (debug_this_module) THEN
545 : CPASSERT(SIZE(S_evects, 1) == nspins)
546 : CPASSERT(SIZE(S_evects, 2) == nvects_total)
547 : END IF
548 :
549 0 : DO ispin = 1, nspins
550 0 : CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
551 : END DO
552 :
553 0 : DO jvect = 1, nvects_total
554 : ! <psi1_i | psi1_j>
555 0 : DO ivect = 1, jvect - 1
556 0 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
557 0 : norm = SUM(weights(1:nspins))
558 :
559 0 : DO ispin = 1, nspins
560 0 : CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect), -norm, evects(ispin, ivect))
561 : END DO
562 : END DO
563 :
564 : ! <psi1_j | psi1_j>
565 0 : DO ispin = 1, nspins
566 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
567 0 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
568 : END DO
569 :
570 0 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
571 :
572 0 : norm = SUM(weights(1:nspins))
573 : norm = 1.0_dp/SQRT(norm)
574 :
575 0 : IF ((ounit > 0) .AND. debug_this_module) WRITE (ounit, '(A,F10.8)') "norm", norm
576 :
577 : END DO
578 :
579 0 : CALL timestop(handle)
580 :
581 0 : END SUBROUTINE tddfpt_check_orthonormality
582 : ! **************************************************************************************************
583 : !> \brief ...
584 : !> \param evects ...
585 : !> \param phase_evects ...
586 : !> \param sub_env ...
587 : ! **************************************************************************************************
588 2 : SUBROUTINE compute_phase_eigenvectors(evects, phase_evects, sub_env)
589 :
590 : ! copied from parts of tddgpt_init_ground_state_mos by S. Chulkov
591 :
592 : TYPE(cp_fm_type), INTENT(in) :: evects
593 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: phase_evects
594 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
595 :
596 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_phase_eigenvectors'
597 : REAL(kind=dp), PARAMETER :: eps_dp = EPSILON(0.0_dp)
598 :
599 : INTEGER :: handle, icol_global, icol_local, irow_global, irow_local, ncol_global, &
600 : ncol_local, nrow_global, nrow_local, sign_int
601 : INTEGER, ALLOCATABLE, DIMENSION(:) :: minrow_neg_array, minrow_pos_array, &
602 : sum_sign_array
603 2 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
604 : REAL(kind=dp) :: element
605 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
606 2 : POINTER :: my_block
607 :
608 2 : CALL timeset(routineN, handle)
609 :
610 : ! compute and print the phase of excited-state eigenvectors:
611 : CALL cp_fm_get_info(evects, nrow_global=nrow_global, ncol_global=ncol_global, &
612 : nrow_local=nrow_local, ncol_local=ncol_local, local_data=my_block, &
613 2 : row_indices=row_indices, col_indices=col_indices) ! nrow_global either nao or nocc
614 :
615 10 : ALLOCATE (minrow_neg_array(ncol_global), minrow_pos_array(ncol_global), sum_sign_array(ncol_global))
616 10 : minrow_neg_array(:) = nrow_global
617 10 : minrow_pos_array(:) = nrow_global
618 2 : sum_sign_array(:) = 0
619 :
620 10 : DO icol_local = 1, ncol_local
621 8 : icol_global = col_indices(icol_local)
622 :
623 86 : DO irow_local = 1, nrow_local
624 76 : irow_global = row_indices(irow_local)
625 :
626 76 : element = my_block(irow_local, icol_local)
627 :
628 76 : sign_int = 0
629 76 : IF (element >= eps_dp) THEN
630 : sign_int = 1
631 36 : ELSE IF (element <= -eps_dp) THEN
632 36 : sign_int = -1
633 : END IF
634 :
635 76 : sum_sign_array(icol_global) = sum_sign_array(icol_global) + sign_int
636 :
637 84 : IF (sign_int > 0) THEN
638 40 : IF (minrow_pos_array(icol_global) > irow_global) THEN
639 8 : minrow_pos_array(icol_global) = irow_global
640 : END IF
641 36 : ELSE IF (sign_int < 0) THEN
642 36 : IF (minrow_neg_array(icol_global) > irow_global) THEN
643 8 : minrow_neg_array(icol_global) = irow_global
644 : END IF
645 : END IF
646 :
647 : END DO
648 : END DO
649 :
650 2 : CALL sub_env%para_env%sum(sum_sign_array)
651 2 : CALL sub_env%para_env%min(minrow_neg_array)
652 2 : CALL sub_env%para_env%min(minrow_pos_array)
653 :
654 10 : DO icol_global = 1, ncol_global
655 :
656 10 : IF (sum_sign_array(icol_global) > 0) THEN
657 : ! most of the expansion coefficients are positive => MO's phase = +1
658 6 : phase_evects(icol_global) = 1.0_dp
659 2 : ELSE IF (sum_sign_array(icol_global) < 0) THEN
660 : ! most of the expansion coefficients are negative => MO's phase = -1
661 2 : phase_evects(icol_global) = -1.0_dp
662 : ELSE
663 : ! equal number of positive and negative expansion coefficients
664 0 : IF (minrow_pos_array(icol_global) <= minrow_neg_array(icol_global)) THEN
665 : ! the first positive expansion coefficient has a lower index then
666 : ! the first negative expansion coefficient; MO's phase = +1
667 0 : phase_evects(icol_global) = 1.0_dp
668 : ELSE
669 : ! MO's phase = -1
670 0 : phase_evects(icol_global) = -1.0_dp
671 : END IF
672 : END IF
673 :
674 : END DO
675 :
676 2 : DEALLOCATE (minrow_neg_array, minrow_pos_array, sum_sign_array)
677 :
678 2 : CALL timestop(handle)
679 :
680 2 : END SUBROUTINE compute_phase_eigenvectors
681 :
682 : END MODULE qs_tddfpt2_restart
|
