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