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