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