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 AO_EXPORT_TYPE keyword
67 : INTEGER, PARAMETER, PUBLIC :: mokp_ao_gto_basis = 1, &
68 : mokp_ao_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 AO_EXPORT_TYPE = OVERLAP_MATRIX)
81 : !> Plus a header with cell, atom info, and either GTO basis set definition
82 : !> (AO_EXPORT_TYPE = GTO_BASIS) or AO metadata list (AO_EXPORT_TYPE = OVERLAP_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, AO_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 :: ao_export_type, atom_shell_index, handle, i, iao, iatom, ikind, ikp, ikp_loc, &
104 : imo, ipgf, iset, isgf_global, ishell, ispin, iw, j, lshell, m_ao, n, nao, natom, ndigits, &
105 : nkp, nmo, nset_atom, nspins, output_unit, unit_choice, z_nuc
106 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: 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, "AO_EXPORT_TYPE", i_val=ao_export_type)
196 0 : write_overlap = (ao_export_type == mokp_ao_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 : ! ================================================================
212 : ! Allocate work buffers
213 : ! ================================================================
214 0 : ALLOCATE (Cre_buf(nao, nmo), Cim_buf(nao, nmo))
215 0 : ALLOCATE (eval_buf(nmo), occ_buf(nmo))
216 :
217 : ! S(k) infrastructure: only needed in MATRIX mode
218 0 : IF (write_overlap) THEN
219 0 : ALLOCATE (Sre_buf(nao, nao), Sim_buf(nao, nao))
220 :
221 0 : ALLOCATE (rmatrix, cmatrix, tmpmat)
222 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
223 0 : matrix_type=dbcsr_type_symmetric)
224 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
225 0 : matrix_type=dbcsr_type_antisymmetric)
226 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
227 0 : matrix_type=dbcsr_type_no_symmetry)
228 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
229 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
230 :
231 0 : NULLIFY (fm_struct_S)
232 : CALL cp_fm_struct_create(fm_struct_S, context=blacs_env_global, &
233 : nrow_global=nao, ncol_global=nao, &
234 0 : para_env=para_env_global)
235 0 : CALL cp_fm_create(fm_S_global, fm_struct_S, name="S(k) work")
236 0 : CALL cp_fm_struct_release(fm_struct_S)
237 : END IF
238 :
239 : ! ================================================================
240 : ! Open output file via print key
241 : ! ================================================================
242 : iw = cp_print_key_unit_nr(logger, print_section, "", &
243 0 : extension=".mokp", file_status="REPLACE")
244 :
245 : ! ================================================================
246 : ! WRITE HEADER
247 : ! ================================================================
248 0 : IF (iw > 0) THEN
249 0 : WRITE (iw, '(A)') "# CP2K_KPOINT_MO_DUMP, Version 2.0"
250 0 : IF (unit_choice == 2) THEN
251 0 : WRITE (iw, '(A)') "# All energies in Hartree, lengths in Angstrom, k-points in fractional reciprocal coords"
252 : ELSE
253 0 : WRITE (iw, '(A)') "# All energies in Hartree, lengths in Bohr, k-points in fractional reciprocal coords"
254 : END IF
255 0 : WRITE (iw, '(A,4I8)') "# DIMENSIONS: natom nspins nao nkp =", natom, nspins, nao, nkp
256 0 : WRITE (iw, '(A,I8)') "# NMO =", nmo
257 0 : WRITE (iw, '(A,L2)') "# USE_REAL_WFN =", use_real_wfn
258 0 : IF (write_overlap) THEN
259 0 : WRITE (iw, '(A)') "# AO_EXPORT_TYPE = OVERLAP_MATRIX"
260 : ELSE
261 0 : WRITE (iw, '(A)') "# AO_EXPORT_TYPE = GTO_BASIS"
262 : END IF
263 0 : WRITE (iw, '(A,ES12.5)') "# SPARSE_THRESHOLD =", thresh
264 0 : IF (unit_choice == 2) THEN
265 0 : WRITE (iw, '(A)') "# CELL_VECTORS [Angstrom]"
266 : ELSE
267 0 : WRITE (iw, '(A)') "# CELL_VECTORS [Bohr]"
268 : END IF
269 : WRITE (iw, '(3(F12.6,3X))') &
270 0 : cell%hmat(1, 1)*scale_factor, cell%hmat(2, 1)*scale_factor, cell%hmat(3, 1)*scale_factor
271 : WRITE (iw, '(3(F12.6,3X))') &
272 0 : cell%hmat(1, 2)*scale_factor, cell%hmat(2, 2)*scale_factor, cell%hmat(3, 2)*scale_factor
273 : WRITE (iw, '(3(F12.6,3X))') &
274 0 : cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
275 : ! Atom table
276 0 : IF (unit_choice == 2) THEN
277 : WRITE (iw, '(A)') &
278 0 : "# ATOM_LIST: Atom_ID Element Z_nuc Z_eff X [Ang] Y [Ang] Z [Ang] First_AO Last_AO"
279 : ELSE
280 : WRITE (iw, '(A)') &
281 0 : "# ATOM_LIST: Atom_ID Element Z_nuc Z_eff X [Bohr] Y [Bohr] Z [Bohr] First_AO Last_AO"
282 : END IF
283 0 : DO iatom = 1, natom
284 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
285 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
286 0 : element_symbol=element_symbol, z=z_nuc)
287 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
288 : WRITE (iw, '(I6,2X,A2,2I6,3F15.6,2I10)') &
289 0 : iatom, element_symbol, z_nuc, NINT(zeff), particle_set(iatom)%r(:)*scale_factor, &
290 0 : first_sgf(iatom), last_sgf(iatom)
291 : END DO
292 :
293 : ! K-point list
294 0 : WRITE (iw, '(A)') "# KPOINT_LIST: ikp kx ky kz weight"
295 0 : DO ikp = 1, nkp
296 0 : WRITE (iw, '(I6,4ES18.8)') ikp, xkp(1:3, ikp), wkp(ikp)
297 : END DO
298 :
299 0 : IF (write_overlap) THEN
300 : ! MATRIX mode: write AO metadata since no GTO basis definition is written.
301 : ! Shell angular momentum ordering: spherical harmonics,
302 : ! CP2K convention m = -l, -l+1, ..., 0, ..., +l
303 0 : WRITE (iw, '(A)') "# AO_LIST: ao_index atom_id atom_shell_index l m"
304 0 : DO iatom = 1, natom
305 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
306 0 : atom_shell_index = 0
307 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
308 0 : IF (ASSOCIATED(orb_basis_set)) THEN
309 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
310 : nset=nset_atom, nshell=nshell_set, l=l_set, &
311 0 : first_sgf=first_sgf_set, last_sgf=last_sgf_set)
312 0 : DO iset = 1, nset_atom
313 0 : DO ishell = 1, nshell_set(iset)
314 0 : atom_shell_index = atom_shell_index + 1
315 0 : lshell = l_set(ishell, iset)
316 0 : DO i = first_sgf_set(ishell, iset), last_sgf_set(ishell, iset)
317 0 : isgf_global = first_sgf(iatom) - 1 + i
318 0 : m_ao = i - first_sgf_set(ishell, iset) - lshell
319 0 : WRITE (iw, '(5I8)') isgf_global, iatom, atom_shell_index, lshell, m_ao
320 : END DO
321 : END DO
322 : END DO
323 : END IF
324 : END DO
325 : ELSE
326 : ! GTO mode: write basis set definition
327 : ! Contraction coefficients denormalized to MOLDEN convention
328 : ! (i.e. for raw unnormalized primitive Gaussians).
329 : ! Shell angular momentum ordering: spherical harmonics,
330 : ! CP2K convention m = -l, -l+1, ..., 0, ..., +l
331 0 : WRITE (iw, '(A)') "# GTO_BASIS (spherical, MOLDEN convention for contraction coefficients)"
332 0 : DO iatom = 1, natom
333 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
334 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
335 0 : IF (ASSOCIATED(orb_basis_set)) THEN
336 0 : WRITE (iw, '(I6,I4)') iatom, 0
337 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
338 : nset=nset_atom, npgf=npgf_set, &
339 : nshell=nshell_set, l=l_set, &
340 0 : zet=zet, gcc=gcc)
341 0 : DO iset = 1, nset_atom
342 0 : DO ishell = 1, nshell_set(iset)
343 0 : lshell = l_set(ishell, iset)
344 0 : IF (lshell + 1 > LEN(angmom)) THEN
345 0 : CPWARN("MOKP: Angular momentum l > 5 not supported in GTO output, skipping")
346 0 : CYCLE
347 : END IF
348 : WRITE (iw, '(A2,I6,F8.2)') &
349 0 : angmom(lshell + 1:lshell + 1)//" ", npgf_set(iset), 1.0_dp
350 : ! Denormalize gcc: undo CP2K's internal normalization
351 : ! (same formula as molden_utils.F, reverse of normalise_gcc_orb)
352 0 : prefac = 2.0_dp**lshell*(2.0_dp/pi)**0.75_dp
353 0 : expzet = 0.25_dp*(2*lshell + 3.0_dp)
354 0 : DO ipgf = 1, npgf_set(iset)
355 : WRITE (iw, fmtstr_gto) &
356 0 : zet(ipgf, iset), &
357 0 : gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet)
358 : END DO
359 : END DO
360 : END DO
361 0 : WRITE (iw, '(A)') ""
362 : END IF
363 : END DO
364 : END IF
365 : END IF
366 :
367 : ! ================================================================
368 : ! WRITE PER-KPOINT DATA
369 : ! ================================================================
370 0 : DO ikp = 1, nkp
371 :
372 : ! ==============================================================
373 : ! A) MO coefficients, eigenvalues, occupations (per spin)
374 : ! ==============================================================
375 0 : DO ispin = 1, nspins
376 :
377 0 : Cre_buf(:, :) = 0.0_dp
378 0 : Cim_buf(:, :) = 0.0_dp
379 0 : eval_buf(:) = 0.0_dp
380 0 : occ_buf(:) = 0.0_dp
381 :
382 0 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
383 0 : ikp_loc = ikp - kp_range(1) + 1
384 0 : kp => kpoints%kp_env(ikp_loc)%kpoint_env
385 0 : mos_kp => kp%mos
386 :
387 : CALL get_mo_set(mos_kp(1, ispin), &
388 : eigenvalues=eigenvalues, &
389 0 : occupation_numbers=occupation_numbers)
390 0 : eval_buf(1:nmo) = eigenvalues(1:nmo)
391 0 : occ_buf(1:nmo) = occupation_numbers(1:nmo)
392 :
393 0 : CALL get_mo_set(mos_kp(1, ispin), mo_coeff=mo_coeff_re)
394 0 : CALL cp_fm_get_submatrix(mo_coeff_re, Cre_buf)
395 :
396 0 : IF (.NOT. use_real_wfn) THEN
397 0 : CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
398 0 : CALL cp_fm_get_submatrix(mo_coeff_im, Cim_buf)
399 : END IF
400 : END IF
401 :
402 0 : CALL para_env_inter_kp%sum(eval_buf)
403 0 : CALL para_env_inter_kp%sum(occ_buf)
404 0 : CALL para_env_inter_kp%sum(Cre_buf)
405 0 : IF (.NOT. use_real_wfn) THEN
406 0 : CALL para_env_inter_kp%sum(Cim_buf)
407 : END IF
408 :
409 0 : IF (iw > 0) THEN
410 0 : WRITE (iw, '(A,2I6)') "# BEGIN_KPOINT_SPIN ikp ispin =", ikp, ispin
411 :
412 0 : WRITE (iw, '(A)') "# EIGENVALUES"
413 0 : WRITE (iw, fmtstr_vec) (eval_buf(n), n=1, nmo)
414 :
415 0 : WRITE (iw, '(A)') "# OCCUPATIONS"
416 0 : WRITE (iw, fmtstr_vec) (occ_buf(n), n=1, nmo)
417 :
418 0 : WRITE (iw, '(A)') "# MO_COEFF_RE (Sparse: mo_index ao_index value)"
419 0 : DO imo = 1, nmo
420 0 : DO iao = 1, nao
421 0 : IF (ABS(Cre_buf(iao, imo)) >= thresh) THEN
422 0 : WRITE (iw, fmtstr_sparse) imo, iao, Cre_buf(iao, imo)
423 : END IF
424 : END DO
425 : END DO
426 :
427 0 : IF (.NOT. use_real_wfn) THEN
428 0 : WRITE (iw, '(A)') "# MO_COEFF_IM (Sparse: mo_index ao_index value)"
429 0 : DO imo = 1, nmo
430 0 : DO iao = 1, nao
431 0 : IF (ABS(Cim_buf(iao, imo)) >= thresh) THEN
432 0 : WRITE (iw, fmtstr_sparse) imo, iao, Cim_buf(iao, imo)
433 : END IF
434 : END DO
435 : END DO
436 : END IF
437 0 : WRITE (iw, '(A,2I6)') "# END_KPOINT_SPIN ikp ispin =", ikp, ispin
438 : END IF
439 :
440 : END DO ! ispin
441 :
442 : ! ==============================================================
443 : ! B) Overlap matrix S(k) — only in MATRIX mode
444 : ! ==============================================================
445 0 : IF (write_overlap) THEN
446 :
447 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
448 0 : IF (use_real_wfn) THEN
449 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
450 : xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, &
451 0 : sab_nl=sab_nl)
452 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
453 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
454 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
455 0 : Sim_buf(:, :) = 0.0_dp
456 : ELSE
457 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
458 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, &
459 : ispin=1, xkp=xkp(1:3, ikp), &
460 0 : cell_to_index=cell_to_index, sab_nl=sab_nl)
461 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
462 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
463 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
464 0 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
465 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
466 0 : CALL cp_fm_get_submatrix(fm_S_global, Sim_buf)
467 : END IF
468 :
469 0 : IF (iw > 0) THEN
470 0 : WRITE (iw, '(A,I6)') "# BEGIN_OVERLAP ikp =", ikp
471 :
472 0 : WRITE (iw, '(A)') "# OVERLAP_RE (Sparse: ao_index_1 ao_index_2 value)"
473 0 : DO i = 1, nao
474 0 : DO j = 1, nao
475 0 : IF (ABS(Sre_buf(i, j)) >= thresh) THEN
476 0 : WRITE (iw, fmtstr_sparse) i, j, Sre_buf(i, j)
477 : END IF
478 : END DO
479 : END DO
480 0 : IF (.NOT. use_real_wfn) THEN
481 0 : WRITE (iw, '(A)') "# OVERLAP_IM (Sparse: ao_index_1 ao_index_2 value)"
482 0 : DO i = 1, nao
483 0 : DO j = 1, nao
484 0 : IF (ABS(Sim_buf(i, j)) >= thresh) THEN
485 0 : WRITE (iw, fmtstr_sparse) i, j, Sim_buf(i, j)
486 : END IF
487 : END DO
488 : END DO
489 : END IF
490 0 : WRITE (iw, '(A,I6)') "# END_OVERLAP ikp =", ikp
491 : END IF
492 :
493 : END IF ! write_overlap
494 :
495 : END DO ! ikp
496 :
497 : ! ================================================================
498 : ! Finalize
499 : ! ================================================================
500 0 : IF (iw > 0) THEN
501 0 : WRITE (iw, '(A)') "# END_OF_FILE"
502 : END IF
503 :
504 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "")
505 :
506 0 : IF (output_unit > 0) THEN
507 : WRITE (output_unit, '(T3,A)') &
508 0 : "KPOINT_MO_DUMP| Data written to .mokp file"
509 : END IF
510 :
511 : ! Release work arrays
512 0 : IF (write_overlap) THEN
513 0 : CALL cp_fm_release(fm_S_global)
514 0 : CALL dbcsr_deallocate_matrix(rmatrix)
515 0 : CALL dbcsr_deallocate_matrix(cmatrix)
516 0 : CALL dbcsr_deallocate_matrix(tmpmat)
517 0 : DEALLOCATE (Sre_buf, Sim_buf)
518 : END IF
519 0 : DEALLOCATE (Cre_buf, Cim_buf)
520 0 : DEALLOCATE (eval_buf, occ_buf)
521 0 : DEALLOCATE (first_sgf, last_sgf)
522 :
523 0 : CALL timestop(handle)
524 :
525 0 : END SUBROUTINE write_kpoint_mo_data
526 :
527 : END MODULE kpoint_mo_dump
|
