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 K-point MO wavefunction dump to TEXT file for post-processing (PDOS, etc.)
10 : !> \par History
11 : !> 2026.02 created
12 : ! **************************************************************************************************
13 :
14 : MODULE kpoint_mo_dump
15 :
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE cell_types, ONLY: cell_type
20 : USE cp_blacs_env, ONLY: cp_blacs_env_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: &
23 : dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_p_type, dbcsr_set, &
24 : dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
25 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
27 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_get_submatrix,&
32 : cp_fm_release,&
33 : cp_fm_type
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_get_default_io_unit,&
36 : cp_logger_type
37 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
38 : cp_print_key_unit_nr
39 : USE input_section_types, ONLY: section_vals_type,&
40 : section_vals_val_get
41 : USE kinds, ONLY: dp
42 : USE kpoint_methods, ONLY: rskp_transform
43 : USE kpoint_types, ONLY: get_kpoint_info,&
44 : kpoint_env_type,&
45 : kpoint_type
46 : USE mathconstants, ONLY: pi
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE particle_methods, ONLY: get_particle_set
49 : USE particle_types, ONLY: particle_type
50 : USE physcon, ONLY: angstrom
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_kind_types, ONLY: get_qs_kind,&
54 : get_qs_kind_set,&
55 : qs_kind_type
56 : USE qs_mo_types, ONLY: get_mo_set,&
57 : mo_set_type
58 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 : PRIVATE
63 :
64 : PUBLIC :: write_kpoint_mo_data
65 :
66 : ! Enum values for OVERLAP_EXPORT_TYPE keyword
67 : INTEGER, PARAMETER, PUBLIC :: mokp_overlap_gto = 1, &
68 : mokp_overlap_matrix = 2
69 :
70 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_mo_dump'
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Write k-point resolved MO data to formatted text file.
76 : !>
77 : !> The output file contains, for each k-point:
78 : !> - Eigenvalues and occupations
79 : !> - MO coefficient matrix C(k) (real and imaginary parts)
80 : !> - Overlap matrix S(k) (if OVERLAP_EXPORT_TYPE = MATRIX)
81 : !> Plus a header with cell, atom info, and either GTO basis set definition
82 : !> (OVERLAP_EXPORT_TYPE = GTO) or AO angular momentum map (OVERLAP_EXPORT_TYPE = MATRIX).
83 : !>
84 : !> Parallel strategy:
85 : !> MO coefficients: cp_fm_get_submatrix (collective on kp-group),
86 : !> then para_env_inter_kp%sum to collect across groups.
87 : !> S(k): Built on global communicator using rskp_transform,
88 : !> then copy_dbcsr_to_fm + cp_fm_get_submatrix.
89 : !> File write: Only on ionode (cp_print_key_unit_nr returns -1 elsewhere).
90 : !>
91 : !> \param qs_env the QS environment (after converged SCF with k-points)
92 : !> \param print_section the &MO_KP print key section (contains NDIGITS, OVERLAP_EXPORT_TYPE)
93 : ! **************************************************************************************************
94 0 : SUBROUTINE write_kpoint_mo_data(qs_env, print_section)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : TYPE(section_vals_type), POINTER :: print_section
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'write_kpoint_mo_data'
99 : CHARACTER(LEN=6), PARAMETER :: angmom = "spdfgh"
100 :
101 : CHARACTER(LEN=2) :: element_symbol
102 : CHARACTER(LEN=20) :: fmtstr_gto, fmtstr_sparse, fmtstr_vec
103 : INTEGER :: handle, i, iatom, ikind, ikp, ikp_loc, ipgf, iset, isgf_global, ishell, ispin, &
104 : iw, j, lshell, n, nao, natom, ndigits, nkp, nmo, nset_atom, nspins, output_unit, &
105 : overlap_data, unit_choice
106 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_l, first_sgf, last_sgf
107 : INTEGER, DIMENSION(2) :: kp_range
108 0 : INTEGER, DIMENSION(:), POINTER :: npgf_set, nshell_set
109 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_set, l_set, last_sgf_set
110 0 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
111 : LOGICAL :: use_real_wfn, write_overlap
112 : REAL(KIND=dp) :: expzet, prefac, scale_factor, thresh, &
113 : zeff
114 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval_buf, occ_buf
115 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Cim_buf, Cre_buf, Sim_buf, Sre_buf
116 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers, wkp
117 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp, zet
118 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
119 : TYPE(cell_type), POINTER :: cell
120 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
121 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_S
122 : TYPE(cp_fm_type) :: fm_S_global
123 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
124 : TYPE(cp_logger_type), POINTER :: logger
125 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
126 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
127 : TYPE(dft_control_type), POINTER :: dft_control
128 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
129 : TYPE(kpoint_env_type), POINTER :: kp
130 : TYPE(kpoint_type), POINTER :: kpoints
131 0 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
132 : TYPE(mp_para_env_type), POINTER :: para_env_global, para_env_inter_kp
133 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
134 0 : POINTER :: sab_nl
135 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
136 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
137 :
138 0 : CALL timeset(routineN, handle)
139 :
140 0 : logger => cp_get_default_logger()
141 0 : output_unit = cp_logger_get_default_io_unit(logger)
142 :
143 : ! ================================================================
144 : ! Gather all required data from qs_env
145 : ! ================================================================
146 0 : NULLIFY (cell, dft_control, kpoints, particle_set, qs_kind_set, matrix_s, &
147 0 : sab_nl, para_env_global, blacs_env_global)
148 :
149 : CALL get_qs_env(qs_env, &
150 : cell=cell, &
151 : dft_control=dft_control, &
152 : kpoints=kpoints, &
153 : particle_set=particle_set, &
154 : qs_kind_set=qs_kind_set, &
155 : matrix_s_kp=matrix_s, &
156 : para_env=para_env_global, &
157 0 : blacs_env=blacs_env_global)
158 :
159 0 : nspins = dft_control%nspins
160 0 : natom = SIZE(particle_set)
161 :
162 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp, &
163 : use_real_wfn=use_real_wfn, kp_range=kp_range, &
164 : sab_nl=sab_nl, cell_to_index=cell_to_index, &
165 0 : para_env_inter_kp=para_env_inter_kp)
166 0 : CPASSERT(ASSOCIATED(sab_nl))
167 :
168 : ! Total number of AOs and MOs
169 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao)
170 :
171 0 : nmo = 0
172 0 : IF (SIZE(kpoints%kp_env) > 0) THEN
173 0 : mos_kp => kpoints%kp_env(1)%kpoint_env%mos
174 0 : CALL get_mo_set(mo_set=mos_kp(1, 1), nmo=nmo)
175 : END IF
176 0 : CALL para_env_inter_kp%max(nmo)
177 :
178 0 : IF (output_unit > 0) THEN
179 0 : WRITE (output_unit, '(/,T3,A)') "KPOINT_MO_DUMP| Writing k-point wavefunction data"
180 : WRITE (output_unit, '(T3,A,I6,A,I6,A,I6,A,I4)') &
181 0 : "KPOINT_MO_DUMP| nao=", nao, " nmo=", nmo, " nkp=", nkp, " nspins=", nspins
182 : END IF
183 :
184 : ! ================================================================
185 : ! Read keywords
186 : ! ================================================================
187 0 : CALL section_vals_val_get(print_section, "UNIT", i_val=unit_choice)
188 0 : IF (unit_choice == 2) THEN
189 : scale_factor = angstrom
190 : ELSE
191 0 : scale_factor = 1.0_dp
192 : END IF
193 0 : CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
194 0 : ndigits = MIN(MAX(3, ndigits), 30)
195 0 : CALL section_vals_val_get(print_section, "OVERLAP_EXPORT_TYPE", i_val=overlap_data)
196 0 : write_overlap = (overlap_data == mokp_overlap_matrix)
197 :
198 : ! Build format strings controlled by NDIGITS
199 0 : WRITE (UNIT=fmtstr_sparse, FMT='("(2I6,1X,ES",I0,".",I0,")")') ndigits + 10, ndigits
200 0 : WRITE (UNIT=fmtstr_vec, FMT='("(5ES",I0,".",I0,")")') ndigits + 10, ndigits
201 0 : WRITE (UNIT=fmtstr_gto, FMT='("(2ES",I0,".",I0,")")') ndigits + 10, ndigits
202 0 : thresh = 10.0_dp**(-ndigits)
203 :
204 : ! ================================================================
205 : ! Build atom-to-AO mapping
206 : ! ================================================================
207 0 : ALLOCATE (first_sgf(natom), last_sgf(natom))
208 : CALL get_particle_set(particle_set, qs_kind_set, &
209 0 : first_sgf=first_sgf, last_sgf=last_sgf)
210 :
211 : ! AO angular momentum map: only needed when overlap is written explicitly,
212 : ! since in GTO mode l info is derivable from the basis set definition.
213 0 : IF (write_overlap) THEN
214 0 : ALLOCATE (ao_l(nao))
215 0 : ao_l(:) = -1
216 0 : DO iatom = 1, natom
217 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
218 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
219 : CALL get_gto_basis_set(orb_basis_set, nset=nset_atom, &
220 : nshell=nshell_set, l=l_set, &
221 0 : first_sgf=first_sgf_set, last_sgf=last_sgf_set)
222 0 : DO iset = 1, nset_atom
223 0 : DO ishell = 1, nshell_set(iset)
224 0 : lshell = l_set(ishell, iset)
225 0 : DO i = first_sgf_set(ishell, iset), last_sgf_set(ishell, iset)
226 0 : isgf_global = first_sgf(iatom) - 1 + i
227 0 : ao_l(isgf_global) = lshell
228 : END DO
229 : END DO
230 : END DO
231 : END DO
232 : END IF
233 :
234 : ! ================================================================
235 : ! Allocate work buffers
236 : ! ================================================================
237 0 : ALLOCATE (Cre_buf(nao, nmo), Cim_buf(nao, nmo))
238 0 : ALLOCATE (eval_buf(nmo), occ_buf(nmo))
239 :
240 : ! S(k) infrastructure: only needed in MATRIX mode
241 0 : IF (write_overlap) THEN
242 0 : ALLOCATE (Sre_buf(nao, nao), Sim_buf(nao, nao))
243 :
244 0 : ALLOCATE (rmatrix, cmatrix, tmpmat)
245 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
246 0 : matrix_type=dbcsr_type_symmetric)
247 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
248 0 : matrix_type=dbcsr_type_antisymmetric)
249 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
250 0 : matrix_type=dbcsr_type_no_symmetry)
251 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
252 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
253 :
254 0 : NULLIFY (fm_struct_S)
255 : CALL cp_fm_struct_create(fm_struct_S, context=blacs_env_global, &
256 : nrow_global=nao, ncol_global=nao, &
257 0 : para_env=para_env_global)
258 0 : CALL cp_fm_create(fm_S_global, fm_struct_S, name="S(k) work")
259 0 : CALL cp_fm_struct_release(fm_struct_S)
260 : END IF
261 :
262 : ! ================================================================
263 : ! Open output file via print key
264 : ! ================================================================
265 : iw = cp_print_key_unit_nr(logger, print_section, "", &
266 0 : extension=".mokp", file_status="REPLACE")
267 :
268 : ! ================================================================
269 : ! WRITE HEADER
270 : ! ================================================================
271 0 : IF (iw > 0) THEN
272 0 : WRITE (iw, '(A)') "# CP2K_KPOINT_MO_DUMP, Version 1.1"
273 0 : IF (write_overlap) THEN
274 0 : WRITE (iw, '(A)') "# EXPORT_MODE: Coefficients + Overlap matrices"
275 : ELSE
276 0 : WRITE (iw, '(A)') "# EXPORT_MODE: GTO + Coefficients"
277 : END IF
278 0 : IF (unit_choice == 2) THEN
279 0 : WRITE (iw, '(A)') "# All energies in Hartree, lengths in Angstrom, k-points in fractional coords"
280 0 : WRITE (iw, '(A)') "# CELL_VECTORS [Angstrom]"
281 : ELSE
282 0 : WRITE (iw, '(A)') "# All energies in Hartree, lengths in Bohr, k-points in fractional coords"
283 0 : WRITE (iw, '(A)') "# CELL_VECTORS [Bohr]"
284 : END IF
285 : WRITE (iw, '(3(F12.6,3X))') &
286 0 : cell%hmat(1, 1)*scale_factor, cell%hmat(2, 1)*scale_factor, cell%hmat(3, 1)*scale_factor
287 : WRITE (iw, '(3(F12.6,3X))') &
288 0 : cell%hmat(1, 2)*scale_factor, cell%hmat(2, 2)*scale_factor, cell%hmat(3, 2)*scale_factor
289 : WRITE (iw, '(3(F12.6,3X))') &
290 0 : cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
291 :
292 0 : WRITE (iw, '(A,4I8)') "# DIMENSIONS: natom nspins nao nkp =", natom, nspins, nao, nkp
293 0 : WRITE (iw, '(A,I8)') "# NMO =", nmo
294 0 : WRITE (iw, '(A,L2)') "# USE_REAL_WFN =", use_real_wfn
295 0 : IF (write_overlap) THEN
296 0 : WRITE (iw, '(A)') "# OVERLAP_EXPORT_TYPE = MATRIX"
297 : ELSE
298 0 : WRITE (iw, '(A)') "# OVERLAP_EXPORT_TYPE = GTO"
299 : END IF
300 :
301 : ! Atom table
302 0 : IF (unit_choice == 2) THEN
303 : WRITE (iw, '(A)') &
304 0 : "# ATOM_LIST: Atom_ID Element Z_eff X [Ang] Y [Ang] Z [Ang] First_AO Last_AO"
305 : ELSE
306 : WRITE (iw, '(A)') &
307 0 : "# ATOM_LIST: Atom_ID Element Z_eff X [Bohr] Y [Bohr] Z [Bohr] First_AO Last_AO"
308 : END IF
309 0 : DO iatom = 1, natom
310 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
311 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
312 0 : element_symbol=element_symbol)
313 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
314 : WRITE (iw, '(I6,2X,A2,I6,3F15.6,2I10)') &
315 0 : iatom, element_symbol, NINT(zeff), particle_set(iatom)%r(:)*scale_factor, &
316 0 : first_sgf(iatom), last_sgf(iatom)
317 : END DO
318 :
319 : ! K-point list
320 0 : WRITE (iw, '(A)') "# KPOINT_LIST: ikp kx ky kz weight"
321 0 : DO ikp = 1, nkp
322 0 : WRITE (iw, '(I6,4ES18.8)') ikp, xkp(1:3, ikp), wkp(ikp)
323 : END DO
324 :
325 0 : IF (write_overlap) THEN
326 : ! MATRIX mode: write AO angular momentum map (needed since no GTO info)
327 0 : WRITE (iw, '(A)') "# AO_L_MAP: iao l (1-based)"
328 0 : DO i = 1, nao
329 0 : WRITE (iw, '(2I6)') i, ao_l(i)
330 : END DO
331 : ELSE
332 : ! GTO mode: write basis set definition
333 : ! Contraction coefficients denormalized to MOLDEN convention
334 : ! (i.e. for raw unnormalized primitive Gaussians).
335 : ! Shell angular momentum ordering: spherical harmonics,
336 : ! CP2K convention m = -l, -l+1, ..., 0, ..., +l
337 0 : WRITE (iw, '(A)') "# GTO_BASIS (spherical, MOLDEN convention for contraction coefficients)"
338 0 : DO iatom = 1, natom
339 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
340 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
341 0 : IF (ASSOCIATED(orb_basis_set)) THEN
342 0 : WRITE (iw, '(I6,I4)') iatom, 0
343 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
344 : nset=nset_atom, npgf=npgf_set, &
345 : nshell=nshell_set, l=l_set, &
346 0 : zet=zet, gcc=gcc)
347 0 : DO iset = 1, nset_atom
348 0 : DO ishell = 1, nshell_set(iset)
349 0 : lshell = l_set(ishell, iset)
350 0 : IF (lshell + 1 > LEN(angmom)) THEN
351 0 : CPWARN("MOKP: Angular momentum l > 5 not supported in GTO output, skipping")
352 0 : CYCLE
353 : END IF
354 : WRITE (iw, '(A2,I6,F8.2)') &
355 0 : angmom(lshell + 1:lshell + 1)//" ", npgf_set(iset), 1.0_dp
356 : ! Denormalize gcc: undo CP2K's internal normalization
357 : ! (same formula as molden_utils.F, reverse of normalise_gcc_orb)
358 0 : prefac = 2.0_dp**lshell*(2.0_dp/pi)**0.75_dp
359 0 : expzet = 0.25_dp*(2*lshell + 3.0_dp)
360 0 : DO ipgf = 1, npgf_set(iset)
361 : WRITE (iw, fmtstr_gto) &
362 0 : zet(ipgf, iset), &
363 0 : gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet)
364 : END DO
365 : END DO
366 : END DO
367 0 : WRITE (iw, '(A)') ""
368 : END IF
369 : END DO
370 : END IF
371 : END IF
372 :
373 : ! ================================================================
374 : ! WRITE PER-KPOINT DATA
375 : ! ================================================================
376 0 : DO ikp = 1, nkp
377 :
378 : ! ==============================================================
379 : ! A) MO coefficients, eigenvalues, occupations (per spin)
380 : ! ==============================================================
381 0 : DO ispin = 1, nspins
382 :
383 0 : Cre_buf(:, :) = 0.0_dp
384 0 : Cim_buf(:, :) = 0.0_dp
385 0 : eval_buf(:) = 0.0_dp
386 0 : occ_buf(:) = 0.0_dp
387 :
388 0 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
389 0 : ikp_loc = ikp - kp_range(1) + 1
390 0 : kp => kpoints%kp_env(ikp_loc)%kpoint_env
391 0 : mos_kp => kp%mos
392 :
393 : CALL get_mo_set(mos_kp(1, ispin), &
394 : eigenvalues=eigenvalues, &
395 0 : occupation_numbers=occupation_numbers)
396 0 : eval_buf(1:nmo) = eigenvalues(1:nmo)
397 0 : occ_buf(1:nmo) = occupation_numbers(1:nmo)
398 :
399 0 : CALL get_mo_set(mos_kp(1, ispin), mo_coeff=mo_coeff_re)
400 0 : CALL cp_fm_get_submatrix(mo_coeff_re, Cre_buf)
401 :
402 0 : IF (.NOT. use_real_wfn) THEN
403 0 : CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
404 0 : CALL cp_fm_get_submatrix(mo_coeff_im, Cim_buf)
405 : END IF
406 : END IF
407 :
408 0 : CALL para_env_inter_kp%sum(eval_buf)
409 0 : CALL para_env_inter_kp%sum(occ_buf)
410 0 : CALL para_env_inter_kp%sum(Cre_buf)
411 0 : IF (.NOT. use_real_wfn) THEN
412 0 : CALL para_env_inter_kp%sum(Cim_buf)
413 : END IF
414 :
415 0 : IF (iw > 0) THEN
416 0 : WRITE (iw, '(A,2I6)') "# BEGIN_KPOINT_SPIN ikp ispin =", ikp, ispin
417 :
418 0 : WRITE (iw, '(A)') "# EIGENVALUES"
419 0 : WRITE (iw, fmtstr_vec) (eval_buf(n), n=1, nmo)
420 :
421 0 : WRITE (iw, '(A)') "# OCCUPATIONS"
422 0 : WRITE (iw, fmtstr_vec) (occ_buf(n), n=1, nmo)
423 :
424 0 : WRITE (iw, '(A)') "# MO_COEFF_RE (Sparse: i_mo j_ao value)"
425 0 : DO i = 1, nmo
426 0 : DO j = 1, nao
427 0 : IF (ABS(Cre_buf(j, i)) >= thresh) THEN
428 0 : WRITE (iw, fmtstr_sparse) i, j, Cre_buf(j, i)
429 : END IF
430 : END DO
431 : END DO
432 :
433 0 : IF (.NOT. use_real_wfn) THEN
434 0 : WRITE (iw, '(A)') "# MO_COEFF_IM (Sparse: i_mo j_ao value)"
435 0 : DO i = 1, nmo
436 0 : DO j = 1, nao
437 0 : IF (ABS(Cim_buf(j, i)) >= thresh) THEN
438 0 : WRITE (iw, fmtstr_sparse) i, j, Cim_buf(j, i)
439 : END IF
440 : END DO
441 : END DO
442 : END IF
443 0 : WRITE (iw, '(A,2I6)') "# END_KPOINT_SPIN ikp ispin =", ikp, ispin
444 : END IF
445 :
446 : END DO ! ispin
447 :
448 : ! ==============================================================
449 : ! B) Overlap matrix S(k) — only in MATRIX mode
450 : ! ==============================================================
451 0 : IF (write_overlap) THEN
452 :
453 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
454 0 : IF (use_real_wfn) THEN
455 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
456 : xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, &
457 0 : sab_nl=sab_nl)
458 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
459 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
460 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
461 0 : Sim_buf(:, :) = 0.0_dp
462 : ELSE
463 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
464 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, &
465 : ispin=1, xkp=xkp(1:3, ikp), &
466 0 : cell_to_index=cell_to_index, sab_nl=sab_nl)
467 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
468 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
469 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
470 0 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
471 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
472 0 : CALL cp_fm_get_submatrix(fm_S_global, Sim_buf)
473 : END IF
474 :
475 0 : IF (iw > 0) THEN
476 0 : WRITE (iw, '(A,I6)') "# BEGIN_OVERLAP ikp =", ikp
477 :
478 0 : WRITE (iw, '(A)') "# OVERLAP_RE (Sparse: i_ao j_ao value)"
479 0 : DO i = 1, nao
480 0 : DO j = 1, nao
481 0 : IF (ABS(Sre_buf(i, j)) >= thresh) THEN
482 0 : WRITE (iw, fmtstr_sparse) i, j, Sre_buf(i, j)
483 : END IF
484 : END DO
485 : END DO
486 0 : IF (.NOT. use_real_wfn) THEN
487 0 : WRITE (iw, '(A)') "# OVERLAP_IM (Sparse: i_ao j_ao value)"
488 0 : DO i = 1, nao
489 0 : DO j = 1, nao
490 0 : IF (ABS(Sim_buf(i, j)) >= thresh) THEN
491 0 : WRITE (iw, fmtstr_sparse) i, j, Sim_buf(i, j)
492 : END IF
493 : END DO
494 : END DO
495 : END IF
496 0 : WRITE (iw, '(A,I6)') "# END_OVERLAP ikp =", ikp
497 : END IF
498 :
499 : END IF ! write_overlap
500 :
501 : END DO ! ikp
502 :
503 : ! ================================================================
504 : ! Finalize
505 : ! ================================================================
506 0 : IF (iw > 0) THEN
507 0 : WRITE (iw, '(A)') "# END_OF_FILE"
508 : END IF
509 :
510 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "")
511 :
512 0 : IF (output_unit > 0) THEN
513 : WRITE (output_unit, '(T3,A)') &
514 0 : "KPOINT_MO_DUMP| Data written to .mokp file"
515 : END IF
516 :
517 : ! Release work arrays
518 0 : IF (write_overlap) THEN
519 0 : CALL cp_fm_release(fm_S_global)
520 0 : CALL dbcsr_deallocate_matrix(rmatrix)
521 0 : CALL dbcsr_deallocate_matrix(cmatrix)
522 0 : CALL dbcsr_deallocate_matrix(tmpmat)
523 0 : DEALLOCATE (Sre_buf, Sim_buf)
524 0 : DEALLOCATE (ao_l)
525 : END IF
526 0 : DEALLOCATE (Cre_buf, Cim_buf)
527 0 : DEALLOCATE (eval_buf, occ_buf)
528 0 : DEALLOCATE (first_sgf, last_sgf)
529 :
530 0 : CALL timestop(handle)
531 :
532 0 : END SUBROUTINE write_kpoint_mo_data
533 :
534 : END MODULE kpoint_mo_dump
|
