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