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 : ! **************************************************************************************************
9 : !> \brief Definition and initialisation of the mo data type.
10 : !> \par History
11 : !> - adapted to the new QS environment data structure (02.04.2002,MK)
12 : !> - set_mo_occupation added (17.04.02,MK)
13 : !> - correct_mo_eigenvalues added (18.04.02,MK)
14 : !> - calculate_density_matrix moved from qs_scf to here (22.04.02,MK)
15 : !> - mo_set_p_type added (23.04.02,MK)
16 : !> - PRIVATE attribute set for TYPE mo_set_type (23.04.02,MK)
17 : !> - started conversion to LSD (1.2003, Joost VandeVondele)
18 : !> - Split of from qs_mo_types (07.2014, JGH)
19 : !> \author Matthias Krack (09.05.2001,MK)
20 : ! **************************************************************************************************
21 : MODULE qs_mo_io
22 :
23 : USE atomic_kind_types, ONLY: get_atomic_kind
24 : USE basis_set_types, ONLY: get_gto_basis_set,&
25 : gto_basis_set_type
26 : USE cp_dbcsr_api, ONLY: dbcsr_binary_write,&
27 : dbcsr_create,&
28 : dbcsr_p_type,&
29 : dbcsr_release,&
30 : dbcsr_type
31 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
32 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
33 : copy_fm_to_dbcsr
34 : USE cp_files, ONLY: close_file,&
35 : open_file
36 : USE cp_fm_types, ONLY: cp_fm_get_info,&
37 : cp_fm_get_submatrix,&
38 : cp_fm_set_all,&
39 : cp_fm_set_submatrix,&
40 : cp_fm_to_fm,&
41 : cp_fm_type,&
42 : cp_fm_write_unformatted
43 : USE cp_log_handling, ONLY: cp_get_default_logger,&
44 : cp_logger_get_default_unit_nr,&
45 : cp_logger_type,&
46 : cp_to_string
47 : USE cp_output_handling, ONLY: cp_p_file,&
48 : cp_print_key_finished_output,&
49 : cp_print_key_generate_filename,&
50 : cp_print_key_should_output,&
51 : cp_print_key_unit_nr
52 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
53 : section_vals_type,&
54 : section_vals_val_get
55 : USE kahan_sum, ONLY: accurate_sum
56 : USE kinds, ONLY: default_path_length,&
57 : default_string_length,&
58 : dp
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE orbital_pointers, ONLY: indco,&
61 : nco,&
62 : nso
63 : USE orbital_symbols, ONLY: cgf_symbol,&
64 : sgf_symbol
65 : USE orbital_transformation_matrices, ONLY: orbtramat
66 : USE particle_types, ONLY: particle_type
67 : USE physcon, ONLY: evolt
68 : USE qs_density_matrices, ONLY: calculate_density_matrix
69 : USE qs_dftb_types, ONLY: qs_dftb_atom_type
70 : USE qs_dftb_utils, ONLY: get_dftb_atom_param
71 : USE qs_kind_types, ONLY: get_qs_kind,&
72 : get_qs_kind_set,&
73 : qs_kind_type
74 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
75 : USE qs_mo_occupation, ONLY: set_mo_occupation
76 : USE qs_mo_types, ONLY: get_mo_set,&
77 : mo_set_type
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 :
82 : PRIVATE
83 :
84 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_io'
85 :
86 : PUBLIC :: wfn_restart_file_name, &
87 : write_rt_mos_to_restart, &
88 : read_rt_mos_from_restart, &
89 : write_dm_binary_restart, &
90 : write_mo_set_to_output_unit, &
91 : write_mo_set_to_restart, &
92 : read_mo_set_from_restart, &
93 : read_mos_restart_low, &
94 : write_mo_set_low
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief ...
100 : !> \param mo_array ...
101 : !> \param particle_set ...
102 : !> \param dft_section ...
103 : !> \param qs_kind_set ...
104 : !> \param matrix_ks ...
105 : ! **************************************************************************************************
106 155421 : SUBROUTINE write_mo_set_to_restart(mo_array, particle_set, dft_section, qs_kind_set, matrix_ks)
107 :
108 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
109 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
110 : TYPE(section_vals_type), POINTER :: dft_section
111 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
112 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
113 : POINTER :: matrix_ks
114 :
115 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_to_restart'
116 : CHARACTER(LEN=30), DIMENSION(2), PARAMETER :: &
117 : keys = ["SCF%PRINT%RESTART_HISTORY", "SCF%PRINT%RESTART "]
118 :
119 : INTEGER :: handle, ikey, ires, ispin
120 : TYPE(cp_logger_type), POINTER :: logger
121 :
122 155421 : CALL timeset(routineN, handle)
123 :
124 155421 : logger => cp_get_default_logger()
125 :
126 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
127 155421 : dft_section, keys(1)), cp_p_file) .OR. &
128 : BTEST(cp_print_key_should_output(logger%iter_info, &
129 : dft_section, keys(2)), cp_p_file)) THEN
130 :
131 16719 : IF (mo_array(1)%use_mo_coeff_b) THEN
132 : ! we are using the dbcsr mo_coeff
133 : ! we copy it to the fm for anycase
134 12310 : DO ispin = 1, SIZE(mo_array)
135 6675 : IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) THEN
136 0 : CPASSERT(.FALSE.)
137 : END IF
138 : CALL copy_dbcsr_to_fm(mo_array(ispin)%mo_coeff_b, &
139 12310 : mo_array(ispin)%mo_coeff) !fm->dbcsr
140 : END DO
141 : END IF
142 :
143 50157 : DO ikey = 1, SIZE(keys)
144 33438 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
145 155421 : dft_section, keys(ikey)), cp_p_file)) THEN
146 : ires = cp_print_key_unit_nr(logger, dft_section, keys(ikey), &
147 : extension=".wfn", file_status="REPLACE", file_action="WRITE", &
148 16731 : do_backup=.TRUE., file_form="UNFORMATTED")
149 16731 : IF (PRESENT(matrix_ks)) THEN
150 : CALL write_mo_set_low(mo_array, particle_set=particle_set, qs_kind_set=qs_kind_set, &
151 5621 : ires=ires, matrix_ks=matrix_ks)
152 : ELSE
153 : CALL write_mo_set_low(mo_array, particle_set=particle_set, qs_kind_set=qs_kind_set, &
154 11110 : ires=ires)
155 : END IF
156 16731 : CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
157 : END IF
158 : END DO
159 : END IF
160 :
161 155421 : CALL timestop(handle)
162 :
163 155421 : END SUBROUTINE write_mo_set_to_restart
164 :
165 : ! **************************************************************************************************
166 : !> \brief calculates density matrix from mo set and writes the density matrix
167 : !> into a binary restart file
168 : !> \param mo_array mos
169 : !> \param dft_section dft input section
170 : !> \param tmpl_matrix template dbcsr matrix
171 : !> \author Mohammad Hossein Bani-Hashemian
172 : ! **************************************************************************************************
173 9831 : SUBROUTINE write_dm_binary_restart(mo_array, dft_section, tmpl_matrix)
174 :
175 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
176 : TYPE(section_vals_type), POINTER :: dft_section
177 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: tmpl_matrix
178 :
179 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_dm_binary_restart'
180 :
181 : CHARACTER(LEN=default_path_length) :: file_name, project_name
182 : INTEGER :: handle, ispin, unit_nr
183 : LOGICAL :: do_dm_restart
184 : REAL(KIND=dp) :: cs_pos
185 : TYPE(cp_logger_type), POINTER :: logger
186 : TYPE(dbcsr_type), POINTER :: matrix_p_tmp
187 :
188 9831 : CALL timeset(routineN, handle)
189 9831 : logger => cp_get_default_logger()
190 9831 : IF (logger%para_env%is_source()) THEN
191 5043 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
192 : ELSE
193 : unit_nr = -1
194 : END IF
195 :
196 9831 : project_name = logger%iter_info%project_name
197 9831 : CALL section_vals_val_get(dft_section, "SCF%PRINT%DM_RESTART_WRITE", l_val=do_dm_restart)
198 9831 : NULLIFY (matrix_p_tmp)
199 :
200 9831 : IF (do_dm_restart) THEN
201 0 : ALLOCATE (matrix_p_tmp)
202 0 : DO ispin = 1, SIZE(mo_array)
203 0 : CALL dbcsr_create(matrix_p_tmp, template=tmpl_matrix(ispin)%matrix, name="DM RESTART")
204 :
205 0 : IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) CPABORT("mo_coeff_b NOT ASSOCIATED")
206 :
207 0 : CALL copy_fm_to_dbcsr(mo_array(ispin)%mo_coeff, mo_array(ispin)%mo_coeff_b)
208 : CALL calculate_density_matrix(mo_array(ispin), matrix_p_tmp, &
209 0 : use_dbcsr=.TRUE., retain_sparsity=.FALSE.)
210 :
211 0 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_SCF_DM_SPIN_", ispin, "_RESTART.dm"
212 0 : cs_pos = dbcsr_checksum(matrix_p_tmp, pos=.TRUE.)
213 0 : IF (unit_nr > 0) THEN
214 0 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
215 : END IF
216 0 : CALL dbcsr_binary_write(matrix_p_tmp, file_name)
217 :
218 0 : CALL dbcsr_release(matrix_p_tmp)
219 : END DO
220 0 : DEALLOCATE (matrix_p_tmp)
221 : END IF
222 :
223 9831 : CALL timestop(handle)
224 :
225 9831 : END SUBROUTINE write_dm_binary_restart
226 :
227 : ! **************************************************************************************************
228 : !> \brief ...
229 : !> \param mo_array ...
230 : !> \param rt_mos ...
231 : !> \param particle_set ...
232 : !> \param dft_section ...
233 : !> \param qs_kind_set ...
234 : ! **************************************************************************************************
235 432 : SUBROUTINE write_rt_mos_to_restart(mo_array, rt_mos, particle_set, dft_section, qs_kind_set)
236 :
237 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
238 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: rt_mos
239 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
240 : TYPE(section_vals_type), POINTER :: dft_section
241 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
242 :
243 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_mos_to_restart'
244 : CHARACTER(LEN=43), DIMENSION(2), PARAMETER :: keys = [ &
245 : "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY", &
246 : "REAL_TIME_PROPAGATION%PRINT%RESTART "]
247 :
248 : INTEGER :: handle, ikey, ires
249 : TYPE(cp_logger_type), POINTER :: logger
250 :
251 432 : CALL timeset(routineN, handle)
252 432 : logger => cp_get_default_logger()
253 :
254 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
255 432 : dft_section, keys(1)), cp_p_file) .OR. &
256 : BTEST(cp_print_key_should_output(logger%iter_info, &
257 : dft_section, keys(2)), cp_p_file)) THEN
258 :
259 324 : DO ikey = 1, SIZE(keys)
260 :
261 216 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
262 432 : dft_section, keys(ikey)), cp_p_file)) THEN
263 : ires = cp_print_key_unit_nr(logger, dft_section, keys(ikey), &
264 : extension=".rtpwfn", file_status="REPLACE", file_action="WRITE", &
265 108 : do_backup=.TRUE., file_form="UNFORMATTED")
266 : CALL write_mo_set_low(mo_array, qs_kind_set=qs_kind_set, particle_set=particle_set, &
267 108 : ires=ires, rt_mos=rt_mos)
268 108 : CALL cp_print_key_finished_output(ires, logger, dft_section, TRIM(keys(ikey)))
269 : END IF
270 : END DO
271 : END IF
272 :
273 432 : CALL timestop(handle)
274 :
275 432 : END SUBROUTINE write_rt_mos_to_restart
276 :
277 : ! **************************************************************************************************
278 : !> \brief ...
279 : !> \param mo_array ...
280 : !> \param qs_kind_set ...
281 : !> \param particle_set ...
282 : !> \param ires ...
283 : !> \param rt_mos ...
284 : !> \param matrix_ks ...
285 : ! **************************************************************************************************
286 16847 : SUBROUTINE write_mo_set_low(mo_array, qs_kind_set, particle_set, ires, rt_mos, matrix_ks)
287 :
288 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mo_array
289 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
290 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
291 : INTEGER :: ires
292 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
293 : OPTIONAL :: rt_mos
294 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
295 : POINTER :: matrix_ks
296 :
297 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_low'
298 :
299 : INTEGER :: handle, iatom, ikind, imat, iset, &
300 : ishell, ispin, lmax, lshell, &
301 : max_block, nao, natom, nmo, nset, &
302 : nset_max, nshell_max, nspin
303 16847 : INTEGER, DIMENSION(:), POINTER :: nset_info, nshell
304 16847 : INTEGER, DIMENSION(:, :), POINTER :: l, nshell_info
305 16847 : INTEGER, DIMENSION(:, :, :), POINTER :: nso_info
306 16847 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupation_numbers
307 : TYPE(cp_fm_type), POINTER :: mo_coeff
308 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
309 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
310 :
311 16847 : CALL timeset(routineN, handle)
312 :
313 16847 : NULLIFY (mo_coeff)
314 16847 : NULLIFY (mo_eigenvalues)
315 16847 : NULLIFY (mo_occupation_numbers)
316 :
317 16847 : nspin = SIZE(mo_array)
318 16847 : nao = mo_array(1)%nao
319 :
320 16847 : IF (ires > 0) THEN
321 : ! Create some info about the basis set first
322 8595 : natom = SIZE(particle_set, 1)
323 8595 : nset_max = 0
324 8595 : nshell_max = 0
325 :
326 54872 : DO iatom = 1, natom
327 46277 : NULLIFY (orb_basis_set, dftb_parameter)
328 46277 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
329 : CALL get_qs_kind(qs_kind_set(ikind), &
330 : basis_set=orb_basis_set, &
331 46277 : dftb_parameter=dftb_parameter)
332 101149 : IF (ASSOCIATED(orb_basis_set)) THEN
333 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
334 : nset=nset, &
335 : nshell=nshell, &
336 38739 : l=l)
337 38739 : nset_max = MAX(nset_max, nset)
338 109695 : DO iset = 1, nset
339 109695 : nshell_max = MAX(nshell_max, nshell(iset))
340 : END DO
341 7538 : ELSEIF (ASSOCIATED(dftb_parameter)) THEN
342 7537 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
343 7537 : nset_max = MAX(nset_max, 1)
344 7537 : nshell_max = MAX(nshell_max, lmax + 1)
345 : ELSE
346 : ! We assume here an atom without a basis set
347 : ! CPABORT("Unknown basis type. ")
348 : END IF
349 : END DO
350 :
351 42975 : ALLOCATE (nso_info(nshell_max, nset_max, natom))
352 318415 : nso_info(:, :, :) = 0
353 :
354 34380 : ALLOCATE (nshell_info(nset_max, natom))
355 142786 : nshell_info(:, :) = 0
356 :
357 25785 : ALLOCATE (nset_info(natom))
358 54872 : nset_info(:) = 0
359 :
360 54872 : DO iatom = 1, natom
361 46277 : NULLIFY (orb_basis_set, dftb_parameter)
362 46277 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
363 : CALL get_qs_kind(qs_kind_set(ikind), &
364 46277 : basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
365 101149 : IF (ASSOCIATED(orb_basis_set)) THEN
366 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
367 : nset=nset, &
368 : nshell=nshell, &
369 38739 : l=l)
370 38739 : nset_info(iatom) = nset
371 109695 : DO iset = 1, nset
372 70956 : nshell_info(iset, iatom) = nshell(iset)
373 210607 : DO ishell = 1, nshell(iset)
374 100912 : lshell = l(ishell, iset)
375 171868 : nso_info(ishell, iset, iatom) = nso(lshell)
376 : END DO
377 : END DO
378 7538 : ELSEIF (ASSOCIATED(dftb_parameter)) THEN
379 7537 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
380 7537 : nset_info(iatom) = 1
381 7537 : nshell_info(1, iatom) = lmax + 1
382 17656 : DO ishell = 1, lmax + 1
383 10119 : lshell = ishell - 1
384 17656 : nso_info(ishell, 1, iatom) = nso(lshell)
385 : END DO
386 : ELSE
387 : ! We assume here an atom without a basis set
388 : ! CPABORT("Unknown basis type. ")
389 : END IF
390 : END DO
391 :
392 8595 : WRITE (ires) natom, nspin, nao, nset_max, nshell_max
393 54872 : WRITE (ires) nset_info
394 142786 : WRITE (ires) nshell_info
395 318415 : WRITE (ires) nso_info
396 :
397 8595 : DEALLOCATE (nset_info)
398 :
399 8595 : DEALLOCATE (nshell_info)
400 :
401 8595 : DEALLOCATE (nso_info)
402 : END IF
403 :
404 : ! Use the ScaLAPACK block size as a default for buffering columns
405 16847 : CALL cp_fm_get_info(mo_array(1)%mo_coeff, ncol_block=max_block)
406 36612 : DO ispin = 1, nspin
407 19765 : mo_coeff => mo_array(ispin)%mo_coeff
408 19765 : nmo = mo_array(ispin)%nmo
409 19765 : IF (nmo > 0) THEN
410 19535 : mo_eigenvalues => mo_array(ispin)%eigenvalues
411 19535 : mo_occupation_numbers => mo_array(ispin)%occupation_numbers
412 19535 : IF (PRESENT(matrix_ks)) THEN
413 : ! With OT: use the Kohn-Sham matrix for the update of the MO eigenvalues
414 : CALL calculate_subspace_eigenvalues(orbitals=mo_coeff, &
415 : ks_matrix=matrix_ks(ispin)%matrix, &
416 6619 : evals_arg=mo_eigenvalues)
417 : END IF
418 19535 : IF (ires > 0) THEN
419 9948 : WRITE (ires) nmo, &
420 9948 : mo_array(ispin)%homo, &
421 9948 : mo_array(ispin)%lfomo, &
422 19896 : mo_array(ispin)%nelectron
423 9948 : WRITE (ires) mo_eigenvalues(1:nmo), mo_occupation_numbers(1:nmo)
424 : END IF
425 : END IF
426 36612 : IF (PRESENT(rt_mos)) THEN
427 414 : DO imat = 2*ispin - 1, 2*ispin
428 414 : CALL cp_fm_write_unformatted(rt_mos(imat), ires)
429 : END DO
430 : ELSE
431 19627 : CALL cp_fm_write_unformatted(mo_coeff, ires)
432 : END IF
433 : END DO
434 :
435 16847 : CALL timestop(handle)
436 :
437 16847 : END SUBROUTINE write_mo_set_low
438 :
439 : ! **************************************************************************************************
440 : !> \brief ...
441 : !> \param filename ...
442 : !> \param exist ...
443 : !> \param section ...
444 : !> \param logger ...
445 : !> \param kp ...
446 : !> \param xas ...
447 : !> \param rtp ...
448 : ! **************************************************************************************************
449 1214 : SUBROUTINE wfn_restart_file_name(filename, exist, section, logger, kp, xas, rtp)
450 : CHARACTER(LEN=default_path_length), INTENT(OUT) :: filename
451 : LOGICAL, INTENT(OUT) :: exist
452 : TYPE(section_vals_type), POINTER :: section
453 : TYPE(cp_logger_type), POINTER :: logger
454 : LOGICAL, INTENT(IN), OPTIONAL :: kp, xas, rtp
455 :
456 : INTEGER :: n_rep_val
457 : LOGICAL :: my_kp, my_rtp, my_xas
458 : TYPE(section_vals_type), POINTER :: print_key
459 :
460 607 : my_kp = .FALSE.
461 607 : my_xas = .FALSE.
462 607 : my_rtp = .FALSE.
463 607 : IF (PRESENT(kp)) my_kp = kp
464 607 : IF (PRESENT(xas)) my_xas = xas
465 607 : IF (PRESENT(rtp)) my_rtp = rtp
466 :
467 607 : exist = .FALSE.
468 607 : CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
469 607 : IF (n_rep_val > 0) THEN
470 462 : CALL section_vals_val_get(section, "WFN_RESTART_FILE_NAME", c_val=filename)
471 : ELSE
472 145 : IF (my_xas) THEN
473 : ! try to read from the filename that is generated automatically from the printkey
474 4 : print_key => section_vals_get_subs_vals(section, "PRINT%RESTART")
475 : filename = cp_print_key_generate_filename(logger, print_key, &
476 4 : extension="", my_local=.FALSE.)
477 141 : ELSE IF (my_rtp) THEN
478 : ! try to read from the filename that is generated automatically from the printkey
479 3 : print_key => section_vals_get_subs_vals(section, "REAL_TIME_PROPAGATION%PRINT%RESTART")
480 : filename = cp_print_key_generate_filename(logger, print_key, &
481 3 : extension=".rtpwfn", my_local=.FALSE.)
482 138 : ELSE IF (my_kp) THEN
483 : ! try to read from the filename that is generated automatically from the printkey
484 14 : print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
485 : filename = cp_print_key_generate_filename(logger, print_key, &
486 14 : extension=".kp", my_local=.FALSE.)
487 : ELSE
488 : ! try to read from the filename that is generated automatically from the printkey
489 124 : print_key => section_vals_get_subs_vals(section, "SCF%PRINT%RESTART")
490 : filename = cp_print_key_generate_filename(logger, print_key, &
491 124 : extension=".wfn", my_local=.FALSE.)
492 : END IF
493 : END IF
494 607 : IF (.NOT. my_xas) THEN
495 601 : INQUIRE (FILE=filename, exist=exist)
496 : END IF
497 :
498 607 : END SUBROUTINE wfn_restart_file_name
499 :
500 : ! **************************************************************************************************
501 : !> \brief ...
502 : !> \param mo_array ...
503 : !> \param qs_kind_set ...
504 : !> \param particle_set ...
505 : !> \param para_env ...
506 : !> \param id_nr ...
507 : !> \param multiplicity ...
508 : !> \param dft_section ...
509 : !> \param natom_mismatch ...
510 : !> \param cdft ...
511 : !> \param out_unit ...
512 : ! **************************************************************************************************
513 521 : SUBROUTINE read_mo_set_from_restart(mo_array, qs_kind_set, particle_set, &
514 : para_env, id_nr, multiplicity, dft_section, natom_mismatch, &
515 : cdft, out_unit)
516 :
517 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array
518 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
519 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
520 : TYPE(mp_para_env_type), POINTER :: para_env
521 : INTEGER, INTENT(IN) :: id_nr, multiplicity
522 : TYPE(section_vals_type), POINTER :: dft_section
523 : LOGICAL, INTENT(OUT), OPTIONAL :: natom_mismatch
524 : LOGICAL, INTENT(IN), OPTIONAL :: cdft
525 : INTEGER, INTENT(IN), OPTIONAL :: out_unit
526 :
527 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_mo_set_from_restart'
528 :
529 : CHARACTER(LEN=default_path_length) :: file_name
530 : INTEGER :: handle, ispin, my_out_unit, natom, &
531 : nspin, restart_unit
532 : LOGICAL :: exist, my_cdft
533 : TYPE(cp_logger_type), POINTER :: logger
534 :
535 521 : CALL timeset(routineN, handle)
536 521 : logger => cp_get_default_logger()
537 521 : my_cdft = .FALSE.
538 521 : IF (PRESENT(cdft)) my_cdft = cdft
539 521 : my_out_unit = -1
540 521 : IF (PRESENT(out_unit)) my_out_unit = out_unit
541 :
542 521 : nspin = SIZE(mo_array)
543 521 : restart_unit = -1
544 :
545 521 : IF (para_env%is_source()) THEN
546 :
547 279 : natom = SIZE(particle_set, 1)
548 279 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger)
549 279 : IF (id_nr /= 0) THEN
550 : ! Is it one of the backup files?
551 1 : file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr))
552 : END IF
553 :
554 : CALL open_file(file_name=file_name, &
555 : file_action="READ", &
556 : file_form="UNFORMATTED", &
557 : file_status="OLD", &
558 279 : unit_number=restart_unit)
559 :
560 : END IF
561 :
562 : CALL read_mos_restart_low(mo_array, para_env=para_env, qs_kind_set=qs_kind_set, &
563 : particle_set=particle_set, natom=natom, &
564 521 : rst_unit=restart_unit, multiplicity=multiplicity, natom_mismatch=natom_mismatch)
565 :
566 521 : IF (PRESENT(natom_mismatch)) THEN
567 : ! read_mos_restart_low only the io_node returns natom_mismatch, must broadcast it
568 491 : CALL para_env%bcast(natom_mismatch)
569 491 : IF (natom_mismatch) THEN
570 0 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
571 0 : CALL timestop(handle)
572 0 : RETURN
573 : END IF
574 : END IF
575 :
576 : ! Close restart file
577 521 : IF (para_env%is_source()) THEN
578 279 : IF (my_out_unit > 0) THEN
579 : WRITE (UNIT=my_out_unit, FMT="(T2,A)") &
580 6 : "WFN_RESTART| Restart file "//TRIM(file_name)//" read"
581 : END IF
582 279 : CALL close_file(unit_number=restart_unit)
583 : END IF
584 :
585 : ! CDFT has no real dft_section and does not need to print
586 521 : IF (.NOT. my_cdft) THEN
587 1396 : DO ispin = 1, nspin
588 : CALL write_mo_set_to_output_unit(mo_array(ispin), qs_kind_set, particle_set, &
589 1396 : dft_section, 4, 0, final_mos=.FALSE.)
590 : END DO
591 : END IF
592 :
593 521 : CALL timestop(handle)
594 :
595 : END SUBROUTINE read_mo_set_from_restart
596 :
597 : ! **************************************************************************************************
598 : !> \brief ...
599 : !> \param mo_array ...
600 : !> \param rt_mos ...
601 : !> \param qs_kind_set ...
602 : !> \param particle_set ...
603 : !> \param para_env ...
604 : !> \param id_nr ...
605 : !> \param multiplicity ...
606 : !> \param dft_section ...
607 : ! **************************************************************************************************
608 8 : SUBROUTINE read_rt_mos_from_restart(mo_array, rt_mos, qs_kind_set, particle_set, para_env, id_nr, multiplicity, dft_section)
609 :
610 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mo_array
611 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: rt_mos
612 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
613 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
614 : TYPE(mp_para_env_type), POINTER :: para_env
615 : INTEGER, INTENT(IN) :: id_nr, multiplicity
616 : TYPE(section_vals_type), POINTER :: dft_section
617 :
618 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_rt_mos_from_restart'
619 :
620 : CHARACTER(LEN=default_path_length) :: file_name
621 : INTEGER :: handle, ispin, natom, nspin, &
622 : restart_unit, unit_nr
623 : LOGICAL :: exist
624 : TYPE(cp_logger_type), POINTER :: logger
625 :
626 8 : CALL timeset(routineN, handle)
627 8 : logger => cp_get_default_logger()
628 :
629 8 : nspin = SIZE(mo_array)
630 8 : restart_unit = -1
631 :
632 8 : IF (para_env%is_source()) THEN
633 :
634 4 : natom = SIZE(particle_set, 1)
635 4 : CALL wfn_restart_file_name(file_name, exist, dft_section, logger, rtp=.TRUE.)
636 4 : IF (id_nr /= 0) THEN
637 : ! Is it one of the backup files?
638 0 : file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr))
639 : END IF
640 :
641 4 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
642 4 : IF (unit_nr > 0) THEN
643 4 : WRITE (unit_nr, '(T2,A)') "Read RTP restart from the file: "//TRIM(file_name)
644 : END IF
645 :
646 : CALL open_file(file_name=file_name, &
647 : file_action="READ", &
648 : file_form="UNFORMATTED", &
649 : file_status="OLD", &
650 4 : unit_number=restart_unit)
651 :
652 : END IF
653 :
654 : CALL read_mos_restart_low(mo_array, rt_mos=rt_mos, para_env=para_env, &
655 : particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom, &
656 8 : rst_unit=restart_unit, multiplicity=multiplicity)
657 :
658 : ! Close restart file
659 8 : IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
660 :
661 16 : DO ispin = 1, nspin
662 : CALL write_mo_set_to_output_unit(mo_array(ispin), qs_kind_set, particle_set, &
663 16 : dft_section, 4, 0, final_mos=.FALSE.)
664 : END DO
665 :
666 8 : CALL timestop(handle)
667 :
668 8 : END SUBROUTINE read_rt_mos_from_restart
669 :
670 : ! **************************************************************************************************
671 : !> \brief Reading the mos from apreviously defined restart file
672 : !> \param mos ...
673 : !> \param para_env ...
674 : !> \param qs_kind_set ...
675 : !> \param particle_set ...
676 : !> \param natom ...
677 : !> \param rst_unit ...
678 : !> \param multiplicity ...
679 : !> \param rt_mos ...
680 : !> \param natom_mismatch ...
681 : !> \par History
682 : !> 12.2007 created [MI]
683 : !> \author MI
684 : ! **************************************************************************************************
685 585 : SUBROUTINE read_mos_restart_low(mos, para_env, qs_kind_set, particle_set, natom, rst_unit, &
686 : multiplicity, rt_mos, natom_mismatch)
687 :
688 : TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
689 : TYPE(mp_para_env_type), POINTER :: para_env
690 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
691 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
692 : INTEGER, INTENT(IN) :: natom, rst_unit
693 : INTEGER, INTENT(in), OPTIONAL :: multiplicity
694 : TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: rt_mos
695 : LOGICAL, INTENT(OUT), OPTIONAL :: natom_mismatch
696 :
697 : INTEGER :: homo, homo_read, i, iatom, ikind, imat, irow, iset, iset_read, ishell, &
698 : ishell_read, iso, ispin, lfomo_read, lmax, lshell, my_mult, nao, nao_read, natom_read, &
699 : nelectron, nelectron_read, nmo, nmo_read, nnshell, nset, nset_max, nshell_max, nspin, &
700 : nspin_read, offset_read
701 585 : INTEGER, DIMENSION(:), POINTER :: nset_info, nshell
702 585 : INTEGER, DIMENSION(:, :), POINTER :: l, nshell_info
703 585 : INTEGER, DIMENSION(:, :, :), POINTER :: nso_info, offset_info
704 : LOGICAL :: minbas, natom_match, use_this
705 585 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read
706 585 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer_read
707 : TYPE(cp_logger_type), POINTER :: logger
708 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
709 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
710 :
711 1170 : logger => cp_get_default_logger()
712 :
713 585 : nspin = SIZE(mos)
714 585 : nao = mos(1)%nao
715 585 : my_mult = 0
716 585 : IF (PRESENT(multiplicity)) my_mult = multiplicity
717 :
718 585 : IF (para_env%is_source()) THEN
719 311 : READ (rst_unit) natom_read, nspin_read, nao_read, nset_max, nshell_max
720 311 : IF (PRESENT(rt_mos)) THEN
721 4 : IF (nspin_read /= nspin) THEN
722 0 : CPABORT("To change nspin is not possible. ")
723 : END IF
724 : ELSE
725 : ! we should allow for restarting with different spin settings
726 307 : IF (nspin_read /= nspin) THEN
727 : WRITE (cp_logger_get_default_unit_nr(logger), *) &
728 0 : "READ RESTART : WARNING : nspin is not equal "
729 : END IF
730 : ! this case needs fixing of homo/lfomo/nelec/occupations ...
731 307 : IF (nspin_read > nspin) THEN
732 0 : CPABORT("Reducing nspin is not possible. ")
733 : END IF
734 : END IF
735 :
736 311 : natom_match = (natom_read == natom)
737 :
738 311 : IF (natom_match) THEN ! actually do the read read
739 :
740 : ! Let's make it possible to change the basis set
741 1555 : ALLOCATE (nso_info(nshell_max, nset_max, natom_read))
742 1244 : ALLOCATE (nshell_info(nset_max, natom_read))
743 933 : ALLOCATE (nset_info(natom_read))
744 1244 : ALLOCATE (offset_info(nshell_max, nset_max, natom_read))
745 :
746 311 : IF (nao_read /= nao) THEN
747 : WRITE (cp_logger_get_default_unit_nr(logger), *) &
748 1 : " READ RESTART : WARNING : DIFFERENT # AOs ", nao, nao_read
749 1 : IF (PRESENT(rt_mos)) &
750 0 : CPABORT("To change basis is not possible. ")
751 : END IF
752 :
753 1140 : READ (rst_unit) nset_info
754 2865 : READ (rst_unit) nshell_info
755 6943 : READ (rst_unit) nso_info
756 :
757 311 : i = 1
758 1140 : DO iatom = 1, natom
759 2681 : DO iset = 1, nset_info(iatom)
760 4900 : DO ishell = 1, nshell_info(iset, iatom)
761 2530 : offset_info(ishell, iset, iatom) = i
762 4071 : i = i + nso_info(ishell, iset, iatom)
763 : END DO
764 : END DO
765 : END DO
766 :
767 933 : ALLOCATE (vecbuffer_read(1, nao_read))
768 :
769 : END IF ! natom_match
770 : END IF ! ionode
771 :
772 : ! make natom_match and natom_mismatch uniform across all nodes
773 585 : CALL para_env%bcast(natom_match)
774 585 : IF (PRESENT(natom_mismatch)) natom_mismatch = .NOT. natom_match
775 : ! handle natom_match false
776 585 : IF (.NOT. natom_match) THEN
777 0 : IF (PRESENT(natom_mismatch)) THEN
778 : WRITE (cp_logger_get_default_unit_nr(logger), *) &
779 0 : " READ RESTART : WARNING : DIFFERENT natom, returning ", natom, natom_read
780 : RETURN
781 : ELSE
782 0 : CPABORT("Incorrect number of atoms in restart file. ")
783 : END IF
784 : END IF
785 :
786 585 : CALL para_env%bcast(nspin_read)
787 :
788 1755 : ALLOCATE (vecbuffer(1, nao))
789 :
790 1602 : DO ispin = 1, nspin
791 :
792 1017 : nmo = mos(ispin)%nmo
793 1017 : homo = mos(ispin)%homo
794 5340 : mos(ispin)%eigenvalues(:) = 0.0_dp
795 5340 : mos(ispin)%occupation_numbers(:) = 0.0_dp
796 1017 : CALL cp_fm_set_all(mos(ispin)%mo_coeff, 0.0_dp)
797 :
798 1017 : IF (para_env%is_source() .AND. (nmo > 0)) THEN
799 536 : READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read
800 2144 : ALLOCATE (eig_read(nmo_read), occ_read(nmo_read))
801 2749 : eig_read = 0.0_dp
802 2749 : occ_read = 0.0_dp
803 :
804 536 : nmo = MIN(nmo, nmo_read)
805 : IF (nmo_read < nmo) &
806 : CALL cp_warn(__LOCATION__, &
807 : "The number of MOs on the restart unit is smaller than the number of "// &
808 : "the allocated MOs. The MO set will be padded with zeros!")
809 536 : IF (nmo_read > nmo) &
810 : CALL cp_warn(__LOCATION__, &
811 : "The number of MOs on the restart unit is greater than the number of "// &
812 7 : "the allocated MOs. The read MO set will be truncated!")
813 :
814 536 : READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read)
815 2725 : mos(ispin)%eigenvalues(1:nmo) = eig_read(1:nmo)
816 2725 : mos(ispin)%occupation_numbers(1:nmo) = occ_read(1:nmo)
817 536 : DEALLOCATE (eig_read, occ_read)
818 :
819 536 : mos(ispin)%homo = homo_read
820 536 : mos(ispin)%lfomo = lfomo_read
821 536 : IF (homo_read > nmo) THEN
822 0 : IF (nelectron_read == mos(ispin)%nelectron) THEN
823 : CALL cp_warn(__LOCATION__, &
824 : "The number of occupied MOs on the restart unit is larger than "// &
825 0 : "the allocated MOs. The read MO set will be truncated and the occupation numbers recalculated!")
826 0 : CALL set_mo_occupation(mo_set=mos(ispin))
827 : ELSE
828 : ! can not make this a warning i.e. homo must be smaller than nmo
829 : ! otherwise e.g. set_mo_occupation will go out of bounds
830 0 : CPABORT("Number of occupied MOs on restart unit larger than allocated MOs. ")
831 : END IF
832 : END IF
833 : END IF
834 :
835 1017 : CALL para_env%bcast(nmo)
836 1017 : CALL para_env%bcast(mos(ispin)%homo)
837 1017 : CALL para_env%bcast(mos(ispin)%lfomo)
838 1017 : CALL para_env%bcast(mos(ispin)%nelectron)
839 9663 : CALL para_env%bcast(mos(ispin)%eigenvalues)
840 9663 : CALL para_env%bcast(mos(ispin)%occupation_numbers)
841 :
842 1017 : IF (PRESENT(rt_mos)) THEN
843 24 : DO imat = 2*ispin - 1, 2*ispin
844 40 : DO i = 1, nmo
845 16 : IF (para_env%is_source()) THEN
846 168 : READ (rst_unit) vecbuffer
847 : ELSE
848 88 : vecbuffer(1, :) = 0.0_dp
849 : END IF
850 656 : CALL para_env%bcast(vecbuffer)
851 : CALL cp_fm_set_submatrix(rt_mos(imat), &
852 32 : vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
853 : END DO
854 : END DO
855 : ELSE
856 5288 : DO i = 1, nmo
857 4279 : IF (para_env%is_source()) THEN
858 146193 : READ (rst_unit) vecbuffer_read
859 : ! now, try to assign the read to the real vector
860 : ! in case the basis set changed this involves some guessing
861 2185 : irow = 1
862 9774 : DO iatom = 1, natom
863 7589 : NULLIFY (orb_basis_set, dftb_parameter, l, nshell)
864 7589 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
865 : CALL get_qs_kind(qs_kind_set(ikind), &
866 7589 : basis_set=orb_basis_set, dftb_parameter=dftb_parameter)
867 7589 : IF (ASSOCIATED(orb_basis_set)) THEN
868 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
869 : nset=nset, &
870 : nshell=nshell, &
871 7541 : l=l)
872 7541 : minbas = .FALSE.
873 48 : ELSEIF (ASSOCIATED(dftb_parameter)) THEN
874 48 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
875 48 : nset = 1
876 48 : minbas = .TRUE.
877 : ELSE
878 : ! assume an atom without basis set
879 : ! CPABORT("Unknown basis set type. ")
880 0 : nset = 0
881 : END IF
882 :
883 7589 : use_this = .TRUE.
884 7589 : iset_read = 1
885 35528 : DO iset = 1, nset
886 18165 : ishell_read = 1
887 18165 : IF (minbas) THEN
888 48 : nnshell = lmax + 1
889 : ELSE
890 18117 : nnshell = nshell(iset)
891 : END IF
892 57458 : DO ishell = 1, nnshell
893 31704 : IF (minbas) THEN
894 72 : lshell = ishell - 1
895 : ELSE
896 31632 : lshell = l(ishell, iset)
897 : END IF
898 31704 : IF (iset_read > nset_info(iatom)) use_this = .FALSE.
899 : IF (use_this) THEN ! avoids out of bound access of the lower line if false
900 31680 : IF (nso(lshell) == nso_info(ishell_read, iset_read, iatom)) THEN
901 31680 : offset_read = offset_info(ishell_read, iset_read, iatom)
902 31680 : ishell_read = ishell_read + 1
903 31680 : IF (ishell_read > nshell_info(iset, iatom)) THEN
904 18161 : ishell_read = 1
905 18161 : iset_read = iset_read + 1
906 : END IF
907 : ELSE
908 : use_this = .FALSE.
909 : END IF
910 : END IF
911 103852 : DO iso = 1, nso(lshell)
912 72148 : IF (use_this) THEN
913 72004 : IF (offset_read - 1 + iso < 1 .OR. offset_read - 1 + iso > nao_read) THEN
914 0 : vecbuffer(1, irow) = 0.0_dp
915 : ELSE
916 72004 : vecbuffer(1, irow) = vecbuffer_read(1, offset_read - 1 + iso)
917 : END IF
918 : ELSE
919 144 : vecbuffer(1, irow) = 0.0_dp
920 : END IF
921 103852 : irow = irow + 1
922 : END DO
923 49869 : use_this = .TRUE.
924 : END DO
925 : END DO
926 : END DO
927 :
928 : ELSE
929 :
930 72764 : vecbuffer(1, :) = 0.0_dp
931 :
932 : END IF
933 :
934 575551 : CALL para_env%bcast(vecbuffer)
935 : CALL cp_fm_set_submatrix(mos(ispin)%mo_coeff, &
936 5288 : vecbuffer, 1, i, nao, 1, transpose=.TRUE.)
937 : END DO
938 : END IF
939 : ! Skip extra MOs if there any
940 1017 : IF (para_env%is_source()) THEN
941 : !ignore nmo = 0
942 539 : IF (nmo > 0) THEN
943 560 : DO i = nmo + 1, nmo_read
944 1504 : READ (rst_unit) vecbuffer_read
945 : END DO
946 : END IF
947 : END IF
948 :
949 1602 : IF (.NOT. PRESENT(rt_mos)) THEN
950 1009 : IF (ispin == 1 .AND. nspin_read < nspin) THEN
951 :
952 0 : mos(ispin + 1)%homo = mos(ispin)%homo
953 0 : mos(ispin + 1)%lfomo = mos(ispin)%lfomo
954 0 : nelectron = mos(ispin)%nelectron
955 0 : IF (my_mult /= 1) THEN
956 : CALL cp_abort(__LOCATION__, &
957 0 : "Restarting an LSD calculation from an LDA wfn only works for multiplicity=1 (singlets).")
958 : END IF
959 0 : IF (mos(ispin + 1)%nelectron < 0) THEN
960 0 : CPABORT("LSD: too few electrons for this multiplisity. ")
961 : END IF
962 0 : mos(ispin + 1)%eigenvalues = mos(ispin)%eigenvalues
963 0 : mos(ispin)%occupation_numbers = mos(ispin)%occupation_numbers/2.0_dp
964 0 : mos(ispin + 1)%occupation_numbers = mos(ispin)%occupation_numbers
965 0 : CALL cp_fm_to_fm(mos(ispin)%mo_coeff, mos(ispin + 1)%mo_coeff)
966 0 : EXIT
967 : END IF
968 : END IF
969 : END DO ! ispin
970 :
971 585 : DEALLOCATE (vecbuffer)
972 :
973 585 : IF (para_env%is_source()) THEN
974 311 : DEALLOCATE (vecbuffer_read)
975 311 : DEALLOCATE (offset_info)
976 311 : DEALLOCATE (nso_info)
977 311 : DEALLOCATE (nshell_info)
978 311 : DEALLOCATE (nset_info)
979 : END IF
980 :
981 1170 : END SUBROUTINE read_mos_restart_low
982 :
983 : ! **************************************************************************************************
984 : !> \brief Write MO information to output file (eigenvalues, occupation numbers, coefficients)
985 : !> \param mo_set ...
986 : !> \param qs_kind_set ...
987 : !> \param particle_set ...
988 : !> \param dft_section ...
989 : !> \param before Digits before the dot
990 : !> \param kpoint An integer that labels the current k point, e.g. its index
991 : !> \param final_mos ...
992 : !> \param spin ...
993 : !> \param solver_method ...
994 : !> \param rtp ...
995 : !> \param cpart ...
996 : !> \param sim_step ...
997 : !> \param umo_set ...
998 : !> \date 15.05.2001
999 : !> \par History:
1000 : !> - Optionally print Cartesian MOs (20.04.2005, MK)
1001 : !> - Revise printout of MO information (05.05.2021, MK)
1002 : !> \par Variables
1003 : !> - after : Number of digits after point.
1004 : !> - before: Number of digits before point.
1005 : !> \author Matthias Krack (MK)
1006 : !> \version 1.1
1007 : ! **************************************************************************************************
1008 6001 : SUBROUTINE write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, &
1009 : dft_section, before, kpoint, final_mos, spin, &
1010 : solver_method, rtp, cpart, sim_step, umo_set)
1011 :
1012 : TYPE(mo_set_type), INTENT(IN) :: mo_set
1013 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1014 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1015 : TYPE(section_vals_type), POINTER :: dft_section
1016 : INTEGER, INTENT(IN) :: before, kpoint
1017 : LOGICAL, INTENT(IN), OPTIONAL :: final_mos
1018 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: spin
1019 : CHARACTER(LEN=2), INTENT(IN), OPTIONAL :: solver_method
1020 : LOGICAL, INTENT(IN), OPTIONAL :: rtp
1021 : INTEGER, INTENT(IN), OPTIONAL :: cpart, sim_step
1022 : TYPE(mo_set_type), INTENT(IN), OPTIONAL :: umo_set
1023 :
1024 : CHARACTER(LEN=12) :: symbol
1025 6001 : CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
1026 : CHARACTER(LEN=14) :: fmtstr5
1027 : CHARACTER(LEN=15) :: energy_str, orbital_str, step_string
1028 : CHARACTER(LEN=2) :: element_symbol, my_solver_method
1029 : CHARACTER(LEN=2*default_string_length) :: name
1030 : CHARACTER(LEN=21) :: vector_str
1031 : CHARACTER(LEN=22) :: fmtstr4
1032 : CHARACTER(LEN=24) :: fmtstr2
1033 : CHARACTER(LEN=25) :: fmtstr1
1034 : CHARACTER(LEN=29) :: fmtstr6
1035 : CHARACTER(LEN=4) :: reim
1036 : CHARACTER(LEN=40) :: fmtstr3
1037 6001 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
1038 : INTEGER :: after, first_mo, from, homo, iatom, icgf, ico, icol, ikind, imo, irow, iset, &
1039 : isgf, ishell, iso, iw, jcol, last_mo, left, lmax, lshell, nao, natom, ncgf, ncol, nmo, &
1040 : nset, nsgf, numo, right, scf_step, to, width
1041 6001 : INTEGER, DIMENSION(:), POINTER :: mo_index_range, nshell
1042 6001 : INTEGER, DIMENSION(:, :), POINTER :: l
1043 : LOGICAL :: ionode, my_final, my_rtp, &
1044 : print_cartesian, print_eigvals, &
1045 : print_eigvecs, print_occup, &
1046 : should_output
1047 : REAL(KIND=dp) :: gap, maxocc
1048 6001 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mo_eigenvalues, mo_occupation_numbers
1049 6001 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
1050 6001 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
1051 : TYPE(cp_fm_type), POINTER :: mo_coeff, umo_coeff
1052 : TYPE(cp_logger_type), POINTER :: logger
1053 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1054 : TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter
1055 :
1056 6001 : NULLIFY (bcgf_symbol)
1057 6001 : NULLIFY (bsgf_symbol)
1058 6001 : NULLIFY (logger)
1059 6001 : NULLIFY (mo_index_range)
1060 6001 : NULLIFY (nshell)
1061 6001 : NULLIFY (mo_coeff)
1062 :
1063 12002 : logger => cp_get_default_logger()
1064 6001 : ionode = logger%para_env%is_source()
1065 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVALUES", l_val=print_eigvals)
1066 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
1067 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occup)
1068 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%CARTESIAN", l_val=print_cartesian)
1069 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%MO_INDEX_RANGE", i_vals=mo_index_range)
1070 6001 : CALL section_vals_val_get(dft_section, "PRINT%MO%NDIGITS", i_val=after)
1071 6001 : after = MIN(MAX(after, 1), 16)
1072 :
1073 : ! Do we print the final MO information after SCF convergence is reached (default: no)
1074 6001 : IF (PRESENT(final_mos)) THEN
1075 5993 : my_final = final_mos
1076 : ELSE
1077 : my_final = .FALSE.
1078 : END IF
1079 :
1080 : ! complex MOS for RTP, no eigenvalues
1081 6001 : my_rtp = .FALSE.
1082 6001 : IF (PRESENT(rtp)) THEN
1083 8 : my_rtp = rtp
1084 : ! print the first time step if MO print required
1085 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
1086 : "PRINT%MO"), cp_p_file) &
1087 8 : .OR. (sim_step == 1)
1088 : ELSE
1089 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
1090 6914 : "PRINT%MO"), cp_p_file) .OR. my_final
1091 : END IF
1092 :
1093 6001 : IF ((.NOT. should_output) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup))) RETURN
1094 :
1095 5078 : IF (my_rtp) THEN
1096 8 : CPASSERT(PRESENT(sim_step))
1097 8 : CPASSERT(PRESENT(cpart))
1098 8 : scf_step = sim_step
1099 8 : IF (cpart == 0) THEN
1100 4 : reim = "IMAG"
1101 : ELSE
1102 4 : reim = "REAL"
1103 : END IF
1104 8 : print_eigvals = .FALSE.
1105 : ELSE
1106 5070 : scf_step = MAX(0, logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1)
1107 : END IF
1108 :
1109 5078 : IF (.NOT. my_final) THEN
1110 3944 : IF (.NOT. my_rtp) THEN
1111 3936 : step_string = " AFTER SCF STEP"
1112 : ELSE
1113 8 : step_string = " AFTER RTP STEP"
1114 : END IF
1115 : END IF
1116 :
1117 5078 : IF (PRESENT(solver_method)) THEN
1118 4924 : my_solver_method = solver_method
1119 : ELSE
1120 : ! Traditional diagonalization is assumed as default solver method
1121 154 : my_solver_method = "TD"
1122 : END IF
1123 :
1124 : ! Retrieve MO information
1125 : CALL get_mo_set(mo_set=mo_set, &
1126 : mo_coeff=mo_coeff, &
1127 : eigenvalues=eigenvalues, &
1128 : occupation_numbers=occupation_numbers, &
1129 : homo=homo, &
1130 : maxocc=maxocc, &
1131 : nao=nao, &
1132 5078 : nmo=nmo)
1133 5078 : IF (PRESENT(umo_set)) THEN
1134 : CALL get_mo_set(mo_set=umo_set, &
1135 : mo_coeff=umo_coeff, &
1136 20 : nmo=numo)
1137 20 : nmo = nmo + numo
1138 : ELSE
1139 5058 : numo = 0
1140 : END IF
1141 15184 : ALLOCATE (mo_eigenvalues(nmo))
1142 43484 : mo_eigenvalues(:) = 0.0_dp
1143 43374 : mo_eigenvalues(1:nmo - numo) = eigenvalues(1:nmo - numo)
1144 10106 : ALLOCATE (mo_occupation_numbers(nmo))
1145 43484 : mo_occupation_numbers(:) = 0.0_dp
1146 43374 : mo_occupation_numbers(1:nmo - numo) = occupation_numbers(1:nmo - numo)
1147 5078 : IF (numo > 0) THEN
1148 : CALL get_mo_set(mo_set=umo_set, &
1149 20 : eigenvalues=eigenvalues)
1150 130 : mo_eigenvalues(nmo - numo + 1:nmo) = eigenvalues(1:numo)
1151 : END IF
1152 :
1153 5078 : IF (print_eigvecs) THEN
1154 13878 : ALLOCATE (smatrix(nao, nmo))
1155 3482 : CALL cp_fm_get_submatrix(mo_coeff, smatrix(1:nao, 1:nmo - numo))
1156 3482 : IF (numo > 0) THEN
1157 14 : CALL cp_fm_get_submatrix(umo_coeff, smatrix(1:nao, nmo - numo + 1:nmo))
1158 : END IF
1159 3482 : IF (.NOT. ionode) THEN
1160 1741 : DEALLOCATE (smatrix)
1161 : END IF
1162 : END IF
1163 :
1164 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO", &
1165 : ignore_should_output=should_output, &
1166 5078 : extension=".MOLog")
1167 :
1168 5078 : IF (iw > 0) THEN
1169 :
1170 2539 : natom = SIZE(particle_set)
1171 2539 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1172 :
1173 : ! Definition of the variable formats
1174 :
1175 2539 : fmtstr1 = "(T2,A,21X, ( X,I5, X))"
1176 2539 : fmtstr2 = "(T2,A,21X, (1X,F . ))"
1177 2539 : fmtstr3 = "(T2,A,I5,1X,I5,1X,A,1X,A6, (1X,F . ))"
1178 :
1179 2539 : width = before + after + 3
1180 2539 : ncol = INT(56/width)
1181 :
1182 2539 : right = MAX((after - 2), 1)
1183 2539 : left = width - right - 5
1184 :
1185 2539 : WRITE (UNIT=fmtstr1(11:12), FMT="(I2)") ncol
1186 2539 : WRITE (UNIT=fmtstr1(14:15), FMT="(I2)") left
1187 2539 : WRITE (UNIT=fmtstr1(21:22), FMT="(I2)") right
1188 :
1189 2539 : WRITE (UNIT=fmtstr2(11:12), FMT="(I2)") ncol
1190 2539 : WRITE (UNIT=fmtstr2(18:19), FMT="(I2)") width - 1
1191 2539 : WRITE (UNIT=fmtstr2(21:22), FMT="(I2)") after
1192 :
1193 2539 : WRITE (UNIT=fmtstr3(27:28), FMT="(I2)") ncol
1194 2539 : WRITE (UNIT=fmtstr3(34:35), FMT="(I2)") width - 1
1195 2539 : WRITE (UNIT=fmtstr3(37:38), FMT="(I2)") after
1196 :
1197 2539 : IF (my_final .OR. (my_solver_method == "TD")) THEN
1198 2539 : energy_str = "EIGENVALUES"
1199 2539 : vector_str = "EIGENVECTORS"
1200 : ELSE
1201 0 : energy_str = "ENERGIES"
1202 0 : vector_str = "COEFFICIENTS"
1203 : END IF
1204 :
1205 2539 : IF (my_rtp) THEN
1206 4 : energy_str = "ZEROS"
1207 4 : vector_str = TRIM(reim)//" RTP COEFFICIENTS"
1208 : END IF
1209 :
1210 2539 : IF (print_eigvecs) THEN
1211 :
1212 1741 : IF (print_cartesian) THEN
1213 :
1214 85 : orbital_str = "CARTESIAN"
1215 :
1216 340 : ALLOCATE (cmatrix(ncgf, ncgf))
1217 28495 : cmatrix = 0.0_dp
1218 :
1219 : ! Transform spherical MOs to Cartesian MOs
1220 85 : icgf = 1
1221 85 : isgf = 1
1222 271 : DO iatom = 1, natom
1223 186 : NULLIFY (orb_basis_set, dftb_parameter)
1224 186 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1225 : CALL get_qs_kind(qs_kind_set(ikind), &
1226 : basis_set=orb_basis_set, &
1227 186 : dftb_parameter=dftb_parameter)
1228 457 : IF (ASSOCIATED(orb_basis_set)) THEN
1229 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1230 : nset=nset, &
1231 : nshell=nshell, &
1232 150 : l=l)
1233 466 : DO iset = 1, nset
1234 944 : DO ishell = 1, nshell(iset)
1235 478 : lshell = l(ishell, iset)
1236 : CALL dgemm("T", "N", nco(lshell), nmo, nso(lshell), 1.0_dp, &
1237 : orbtramat(lshell)%c2s, nso(lshell), &
1238 : smatrix(isgf, 1), nsgf, 0.0_dp, &
1239 478 : cmatrix(icgf, 1), ncgf)
1240 478 : icgf = icgf + nco(lshell)
1241 794 : isgf = isgf + nso(lshell)
1242 : END DO
1243 : END DO
1244 36 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
1245 36 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
1246 90 : DO ishell = 1, lmax + 1
1247 54 : lshell = ishell - 1
1248 : CALL dgemm("T", "N", nco(lshell), nsgf, nso(lshell), 1.0_dp, &
1249 : orbtramat(lshell)%c2s, nso(lshell), &
1250 : smatrix(isgf, 1), nsgf, 0.0_dp, &
1251 54 : cmatrix(icgf, 1), ncgf)
1252 54 : icgf = icgf + nco(lshell)
1253 90 : isgf = isgf + nso(lshell)
1254 : END DO
1255 : ELSE
1256 : ! assume atom without basis set
1257 : ! CPABORT("Unknown basis set type")
1258 : END IF
1259 : END DO ! iatom
1260 :
1261 : ELSE
1262 :
1263 1656 : orbital_str = "SPHERICAL"
1264 :
1265 : END IF ! print_cartesian
1266 :
1267 : name = TRIM(energy_str)//", OCCUPATION NUMBERS, AND "// &
1268 1741 : TRIM(orbital_str)//" "//TRIM(vector_str)
1269 :
1270 1741 : IF (.NOT. my_final) &
1271 1397 : WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
1272 :
1273 798 : ELSE IF (print_occup .OR. print_eigvals) THEN
1274 798 : name = TRIM(energy_str)//" AND OCCUPATION NUMBERS"
1275 :
1276 798 : IF (.NOT. my_final) &
1277 575 : WRITE (UNIT=name, FMT="(A,1X,I0)") TRIM(name)//step_string, scf_step
1278 : END IF ! print_eigvecs
1279 :
1280 : ! Print headline
1281 2539 : IF (PRESENT(spin) .AND. (kpoint > 0)) THEN
1282 : WRITE (UNIT=iw, FMT="(/,T2,A,I0)") &
1283 0 : "MO| "//TRIM(spin)//" "//TRIM(name)//" FOR K POINT ", kpoint
1284 : ELSE IF (PRESENT(spin)) THEN
1285 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1286 312 : "MO| "//TRIM(spin)//" "//TRIM(name)
1287 2227 : ELSE IF (kpoint > 0) THEN
1288 : WRITE (UNIT=iw, FMT="(/,T2,A,I0)") &
1289 101 : "MO| "//TRIM(name)//" FOR K POINT ", kpoint
1290 : ELSE
1291 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1292 2126 : "MO| "//TRIM(name)
1293 : END IF
1294 :
1295 : ! Check if only a subset of the MOs has to be printed
1296 : IF ((mo_index_range(1) > 0) .AND. &
1297 2539 : (mo_index_range(2) > 0) .AND. &
1298 : (mo_index_range(2) >= mo_index_range(1))) THEN
1299 142 : first_mo = MAX(1, mo_index_range(1))
1300 142 : last_mo = MIN(nmo, mo_index_range(2))
1301 142 : IF (first_mo > nmo) THEN
1302 0 : CPWARN("The first MO index is larger than the number of MOs.")
1303 : END IF
1304 : ELSE
1305 2397 : first_mo = 1
1306 2397 : last_mo = nmo
1307 : END IF
1308 :
1309 2539 : IF (print_eigvecs) THEN
1310 :
1311 : ! Print full MO information
1312 :
1313 3928 : DO icol = first_mo, last_mo, ncol
1314 :
1315 2187 : from = icol
1316 2187 : to = MIN((from + ncol - 1), last_mo)
1317 :
1318 2187 : WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1319 : WRITE (UNIT=iw, FMT=fmtstr1) &
1320 10146 : "MO|", (jcol, jcol=from, to)
1321 : WRITE (UNIT=iw, FMT=fmtstr2) &
1322 2187 : "MO|", (mo_eigenvalues(jcol), jcol=from, to)
1323 2187 : WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1324 : WRITE (UNIT=iw, FMT=fmtstr2) &
1325 2187 : "MO|", (mo_occupation_numbers(jcol), jcol=from, to)
1326 2187 : WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1327 :
1328 2187 : irow = 1
1329 :
1330 8639 : DO iatom = 1, natom
1331 :
1332 4711 : IF (iatom /= 1) WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1333 :
1334 4711 : NULLIFY (orb_basis_set, dftb_parameter)
1335 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
1336 4711 : element_symbol=element_symbol, kind_number=ikind)
1337 : CALL get_qs_kind(qs_kind_set(ikind), &
1338 : basis_set=orb_basis_set, &
1339 4711 : dftb_parameter=dftb_parameter)
1340 :
1341 11609 : IF (print_cartesian) THEN
1342 :
1343 852 : IF (ASSOCIATED(orb_basis_set)) THEN
1344 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1345 : nset=nset, &
1346 : nshell=nshell, &
1347 : l=l, &
1348 744 : cgf_symbol=bcgf_symbol)
1349 :
1350 744 : icgf = 1
1351 2520 : DO iset = 1, nset
1352 4944 : DO ishell = 1, nshell(iset)
1353 2424 : lshell = l(ishell, iset)
1354 10824 : DO ico = 1, nco(lshell)
1355 : WRITE (UNIT=iw, FMT=fmtstr3) &
1356 6624 : "MO|", irow, iatom, ADJUSTR(element_symbol), bcgf_symbol(icgf), &
1357 13248 : (cmatrix(irow, jcol), jcol=from, to)
1358 6624 : icgf = icgf + 1
1359 9048 : irow = irow + 1
1360 : END DO
1361 : END DO
1362 : END DO
1363 108 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
1364 108 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
1365 108 : icgf = 1
1366 270 : DO ishell = 1, lmax + 1
1367 162 : lshell = ishell - 1
1368 540 : DO ico = 1, nco(lshell)
1369 270 : symbol = cgf_symbol(1, indco(1:3, icgf))
1370 270 : symbol(1:2) = " "
1371 : WRITE (UNIT=iw, FMT=fmtstr3) &
1372 270 : "MO|", irow, iatom, ADJUSTR(element_symbol), symbol, &
1373 540 : (cmatrix(irow, jcol), jcol=from, to)
1374 270 : icgf = icgf + 1
1375 432 : irow = irow + 1
1376 : END DO
1377 : END DO
1378 : ELSE
1379 : ! assume atom without basis set
1380 : ! CPABORT("Unknown basis set type")
1381 : END IF
1382 :
1383 : ELSE
1384 :
1385 3859 : IF (ASSOCIATED(orb_basis_set)) THEN
1386 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1387 : nset=nset, &
1388 : nshell=nshell, &
1389 : l=l, &
1390 3859 : sgf_symbol=bsgf_symbol)
1391 3859 : isgf = 1
1392 11027 : DO iset = 1, nset
1393 20499 : DO ishell = 1, nshell(iset)
1394 9472 : lshell = l(ishell, iset)
1395 39040 : DO iso = 1, nso(lshell)
1396 : WRITE (UNIT=iw, FMT=fmtstr3) &
1397 22400 : "MO|", irow, iatom, ADJUSTR(element_symbol), bsgf_symbol(isgf), &
1398 44800 : (smatrix(irow, jcol), jcol=from, to)
1399 22400 : isgf = isgf + 1
1400 31872 : irow = irow + 1
1401 : END DO
1402 : END DO
1403 : END DO
1404 0 : ELSE IF (ASSOCIATED(dftb_parameter)) THEN
1405 0 : CALL get_dftb_atom_param(dftb_parameter, lmax=lmax)
1406 0 : isgf = 1
1407 0 : DO ishell = 1, lmax + 1
1408 0 : lshell = ishell - 1
1409 0 : DO iso = 1, nso(lshell)
1410 0 : symbol = sgf_symbol(1, lshell, -lshell + iso - 1)
1411 0 : symbol(1:2) = " "
1412 : WRITE (UNIT=iw, FMT=fmtstr3) &
1413 0 : "MO|", irow, iatom, ADJUSTR(element_symbol), symbol, &
1414 0 : (smatrix(irow, jcol), jcol=from, to)
1415 0 : isgf = isgf + 1
1416 0 : irow = irow + 1
1417 : END DO
1418 : END DO
1419 : ELSE
1420 : ! assume atom without basis set
1421 : ! CPABORT("Unknown basis set type")
1422 : END IF
1423 :
1424 : END IF ! print_cartesian
1425 :
1426 : END DO ! iatom
1427 :
1428 : END DO ! icol
1429 :
1430 1741 : WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1431 :
1432 : ! Release work storage
1433 :
1434 1741 : DEALLOCATE (smatrix)
1435 1741 : IF (print_cartesian) THEN
1436 85 : DEALLOCATE (cmatrix)
1437 : END IF
1438 :
1439 798 : ELSE IF (print_occup .OR. print_eigvals) THEN
1440 :
1441 798 : WRITE (UNIT=iw, FMT="(T2,A)") "MO|"
1442 798 : fmtstr4 = "(T2,A,I7,3(1X,F22. ))"
1443 798 : WRITE (UNIT=fmtstr4(19:20), FMT="(I2)") after
1444 798 : IF (my_final .OR. (my_solver_method == "TD")) THEN
1445 : WRITE (UNIT=iw, FMT="(A)") &
1446 798 : " MO| Index Eigenvalue [a.u.] Eigenvalue [eV] Occupation"
1447 : ELSE
1448 : WRITE (UNIT=iw, FMT="(A)") &
1449 0 : " MO| Index Energy [a.u.] Energy [eV] Occupation"
1450 : END IF
1451 11998 : DO imo = first_mo, last_mo
1452 : WRITE (UNIT=iw, FMT=fmtstr4) &
1453 11200 : "MO|", imo, mo_eigenvalues(imo), &
1454 11200 : mo_eigenvalues(imo)*evolt, &
1455 23198 : mo_occupation_numbers(imo)
1456 : END DO
1457 798 : fmtstr5 = "(A,T59,F22. )"
1458 798 : WRITE (UNIT=fmtstr5(12:13), FMT="(I2)") after
1459 : WRITE (UNIT=iw, FMT=fmtstr5) &
1460 798 : " MO| Sum:", accurate_sum(mo_occupation_numbers(:))
1461 :
1462 : END IF ! print_eigvecs
1463 :
1464 2539 : IF (.NOT. my_rtp) THEN
1465 2535 : fmtstr6 = "(A,T18,F17. ,A,T41,F17. ,A)"
1466 2535 : WRITE (UNIT=fmtstr6(12:13), FMT="(I2)") after
1467 2535 : WRITE (UNIT=fmtstr6(25:26), FMT="(I2)") after
1468 : WRITE (UNIT=iw, FMT=fmtstr6) &
1469 2535 : " MO| E(Fermi):", mo_set%mu, " a.u.", mo_set%mu*evolt, " eV"
1470 : END IF
1471 2539 : IF ((homo > 0) .AND. .NOT. my_rtp) THEN
1472 2510 : IF ((mo_occupation_numbers(homo) == maxocc) .AND. (last_mo > homo)) THEN
1473 : gap = mo_eigenvalues(homo + 1) - &
1474 351 : mo_eigenvalues(homo)
1475 : WRITE (UNIT=iw, FMT=fmtstr6) &
1476 351 : " MO| Band gap:", gap, " a.u.", gap*evolt, " eV"
1477 : END IF
1478 : END IF
1479 2539 : WRITE (UNIT=iw, FMT="(A)") ""
1480 :
1481 : END IF ! iw
1482 :
1483 5078 : IF (ALLOCATED(mo_eigenvalues)) DEALLOCATE (mo_eigenvalues)
1484 5078 : IF (ALLOCATED(mo_occupation_numbers)) DEALLOCATE (mo_occupation_numbers)
1485 :
1486 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%MO", &
1487 5078 : ignore_should_output=should_output)
1488 :
1489 23081 : END SUBROUTINE write_mo_set_to_output_unit
1490 :
1491 : END MODULE qs_mo_io
|
