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
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, thresh, zeff
113 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval_buf, occ_buf
114 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Cim_buf, Cre_buf, Sim_buf, Sre_buf
115 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers, wkp
116 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp, zet
117 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
118 : TYPE(cell_type), POINTER :: cell
119 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
120 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_S
121 : TYPE(cp_fm_type) :: fm_S_global
122 : TYPE(cp_fm_type), POINTER :: mo_coeff_im, mo_coeff_re
123 : TYPE(cp_logger_type), POINTER :: logger
124 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
125 : TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
126 : TYPE(dft_control_type), POINTER :: dft_control
127 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
128 : TYPE(kpoint_env_type), POINTER :: kp
129 : TYPE(kpoint_type), POINTER :: kpoints
130 0 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
131 : TYPE(mp_para_env_type), POINTER :: para_env_global, para_env_inter_kp
132 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
133 0 : POINTER :: sab_nl
134 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
135 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
136 :
137 0 : CALL timeset(routineN, handle)
138 :
139 0 : logger => cp_get_default_logger()
140 0 : output_unit = cp_logger_get_default_io_unit(logger)
141 :
142 : ! ================================================================
143 : ! Gather all required data from qs_env
144 : ! ================================================================
145 0 : NULLIFY (cell, dft_control, kpoints, particle_set, qs_kind_set, matrix_s, &
146 0 : sab_nl, para_env_global, blacs_env_global)
147 :
148 : CALL get_qs_env(qs_env, &
149 : cell=cell, &
150 : dft_control=dft_control, &
151 : kpoints=kpoints, &
152 : particle_set=particle_set, &
153 : qs_kind_set=qs_kind_set, &
154 : matrix_s_kp=matrix_s, &
155 : para_env=para_env_global, &
156 0 : blacs_env=blacs_env_global)
157 :
158 0 : nspins = dft_control%nspins
159 0 : natom = SIZE(particle_set)
160 :
161 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp, &
162 : use_real_wfn=use_real_wfn, kp_range=kp_range, &
163 : sab_nl=sab_nl, cell_to_index=cell_to_index, &
164 0 : para_env_inter_kp=para_env_inter_kp)
165 0 : CPASSERT(ASSOCIATED(sab_nl))
166 :
167 : ! Total number of AOs and MOs
168 0 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao)
169 :
170 0 : nmo = 0
171 0 : IF (SIZE(kpoints%kp_env) > 0) THEN
172 0 : mos_kp => kpoints%kp_env(1)%kpoint_env%mos
173 0 : CALL get_mo_set(mo_set=mos_kp(1, 1), nmo=nmo)
174 : END IF
175 0 : CALL para_env_inter_kp%max(nmo)
176 :
177 0 : IF (output_unit > 0) THEN
178 0 : WRITE (output_unit, '(/,T3,A)') "KPOINT_MO_DUMP| Writing k-point wavefunction data"
179 : WRITE (output_unit, '(T3,A,I6,A,I6,A,I6,A,I4)') &
180 0 : "KPOINT_MO_DUMP| nao=", nao, " nmo=", nmo, " nkp=", nkp, " nspins=", nspins
181 : END IF
182 :
183 : ! ================================================================
184 : ! Read keywords
185 : ! ================================================================
186 0 : CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
187 0 : ndigits = MIN(MAX(3, ndigits), 30)
188 0 : CALL section_vals_val_get(print_section, "OVERLAP_EXPORT_TYPE", i_val=overlap_data)
189 0 : write_overlap = (overlap_data == mokp_overlap_matrix)
190 :
191 : ! Build format strings controlled by NDIGITS
192 0 : WRITE (UNIT=fmtstr_sparse, FMT='("(2I6,1X,ES",I0,".",I0,")")') ndigits + 10, ndigits
193 0 : WRITE (UNIT=fmtstr_vec, FMT='("(5ES",I0,".",I0,")")') ndigits + 10, ndigits
194 0 : WRITE (UNIT=fmtstr_gto, FMT='("(2ES",I0,".",I0,")")') ndigits + 10, ndigits
195 0 : thresh = 10.0_dp**(-ndigits)
196 :
197 : ! ================================================================
198 : ! Build atom-to-AO mapping
199 : ! ================================================================
200 0 : ALLOCATE (first_sgf(natom), last_sgf(natom))
201 : CALL get_particle_set(particle_set, qs_kind_set, &
202 0 : first_sgf=first_sgf, last_sgf=last_sgf)
203 :
204 : ! AO angular momentum map: only needed when overlap is written explicitly,
205 : ! since in GTO mode l info is derivable from the basis set definition.
206 0 : IF (write_overlap) THEN
207 0 : ALLOCATE (ao_l(nao))
208 0 : ao_l(:) = -1
209 0 : DO iatom = 1, natom
210 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
211 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
212 : CALL get_gto_basis_set(orb_basis_set, nset=nset_atom, &
213 : nshell=nshell_set, l=l_set, &
214 0 : first_sgf=first_sgf_set, last_sgf=last_sgf_set)
215 0 : DO iset = 1, nset_atom
216 0 : DO ishell = 1, nshell_set(iset)
217 0 : lshell = l_set(ishell, iset)
218 0 : DO i = first_sgf_set(ishell, iset), last_sgf_set(ishell, iset)
219 0 : isgf_global = first_sgf(iatom) - 1 + i
220 0 : ao_l(isgf_global) = lshell
221 : END DO
222 : END DO
223 : END DO
224 : END DO
225 : END IF
226 :
227 : ! ================================================================
228 : ! Allocate work buffers
229 : ! ================================================================
230 0 : ALLOCATE (Cre_buf(nao, nmo), Cim_buf(nao, nmo))
231 0 : ALLOCATE (eval_buf(nmo), occ_buf(nmo))
232 :
233 : ! S(k) infrastructure: only needed in MATRIX mode
234 0 : IF (write_overlap) THEN
235 0 : ALLOCATE (Sre_buf(nao, nao), Sim_buf(nao, nao))
236 :
237 0 : ALLOCATE (rmatrix, cmatrix, tmpmat)
238 : CALL dbcsr_create(rmatrix, template=matrix_s(1, 1)%matrix, &
239 0 : matrix_type=dbcsr_type_symmetric)
240 : CALL dbcsr_create(cmatrix, template=matrix_s(1, 1)%matrix, &
241 0 : matrix_type=dbcsr_type_antisymmetric)
242 : CALL dbcsr_create(tmpmat, template=matrix_s(1, 1)%matrix, &
243 0 : matrix_type=dbcsr_type_no_symmetry)
244 0 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
245 0 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
246 :
247 0 : NULLIFY (fm_struct_S)
248 : CALL cp_fm_struct_create(fm_struct_S, context=blacs_env_global, &
249 : nrow_global=nao, ncol_global=nao, &
250 0 : para_env=para_env_global)
251 0 : CALL cp_fm_create(fm_S_global, fm_struct_S, name="S(k) work")
252 0 : CALL cp_fm_struct_release(fm_struct_S)
253 : END IF
254 :
255 : ! ================================================================
256 : ! Open output file via print key
257 : ! ================================================================
258 : iw = cp_print_key_unit_nr(logger, print_section, "", &
259 0 : extension=".mokp", file_status="REPLACE")
260 :
261 : ! ================================================================
262 : ! WRITE HEADER
263 : ! ================================================================
264 0 : IF (iw > 0) THEN
265 0 : WRITE (iw, '(A)') "# CP2K_KPOINT_MO_DUMP, Version 1.0"
266 0 : IF (write_overlap) THEN
267 0 : WRITE (iw, '(A)') "# EXPORT_MODE: Coefficients + Overlap matrices"
268 : ELSE
269 0 : WRITE (iw, '(A)') "# EXPORT_MODE: GTO + Coefficients"
270 : END IF
271 0 : WRITE (iw, '(A)') "# All energies in Hartree, lengths in Angstrom, k-points in fractional coords"
272 :
273 0 : WRITE (iw, '(A)') "# CELL_VECTORS [Angstrom]"
274 : WRITE (iw, '(3ES16.6)') &
275 0 : cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom
276 : WRITE (iw, '(3ES16.6)') &
277 0 : cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom
278 : WRITE (iw, '(3ES16.6)') &
279 0 : cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom
280 :
281 0 : WRITE (iw, '(A,4I8)') "# DIMENSIONS: natom nspins nao nkp =", natom, nspins, nao, nkp
282 0 : WRITE (iw, '(A,I8)') "# NMO =", nmo
283 0 : WRITE (iw, '(A,L2)') "# USE_REAL_WFN =", use_real_wfn
284 0 : IF (write_overlap) THEN
285 0 : WRITE (iw, '(A)') "# OVERLAP_EXPORT_TYPE = MATRIX"
286 : ELSE
287 0 : WRITE (iw, '(A)') "# OVERLAP_EXPORT_TYPE = GTO"
288 : END IF
289 :
290 : ! Atom table
291 : WRITE (iw, '(A)') &
292 0 : "# ATOM_LIST: Atom_ID Element Z_eff X [Ang] Y [Ang] Z [Ang] First_AO Last_AO"
293 0 : DO iatom = 1, natom
294 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
295 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
296 0 : element_symbol=element_symbol)
297 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
298 : WRITE (iw, '(I6,2X,A2,I6,3F15.6,2I10)') &
299 0 : iatom, element_symbol, NINT(zeff), particle_set(iatom)%r(:)*angstrom, &
300 0 : first_sgf(iatom), last_sgf(iatom)
301 : END DO
302 :
303 : ! K-point list
304 0 : WRITE (iw, '(A)') "# KPOINT_LIST: ikp kx ky kz weight"
305 0 : DO ikp = 1, nkp
306 0 : WRITE (iw, '(I6,4ES18.8)') ikp, xkp(1:3, ikp), wkp(ikp)
307 : END DO
308 :
309 0 : IF (write_overlap) THEN
310 : ! MATRIX mode: write AO angular momentum map (needed since no GTO info)
311 0 : WRITE (iw, '(A)') "# AO_L_MAP: iao l (1-based)"
312 0 : DO i = 1, nao
313 0 : WRITE (iw, '(2I6)') i, ao_l(i)
314 : END DO
315 : ELSE
316 : ! GTO mode: write basis set definition
317 : ! Contraction coefficients denormalized to MOLDEN convention
318 : ! (i.e. for raw unnormalized primitive Gaussians).
319 : ! Shell angular momentum ordering: spherical harmonics,
320 : ! CP2K convention m = -l, -l+1, ..., 0, ..., +l
321 0 : WRITE (iw, '(A)') "# GTO_BASIS (spherical, MOLDEN convention for contraction coefficients)"
322 0 : DO iatom = 1, natom
323 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
324 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
325 0 : IF (ASSOCIATED(orb_basis_set)) THEN
326 0 : WRITE (iw, '(I6,I4)') iatom, 0
327 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
328 : nset=nset_atom, npgf=npgf_set, &
329 : nshell=nshell_set, l=l_set, &
330 0 : zet=zet, gcc=gcc)
331 0 : DO iset = 1, nset_atom
332 0 : DO ishell = 1, nshell_set(iset)
333 0 : lshell = l_set(ishell, iset)
334 0 : IF (lshell + 1 > LEN(angmom)) THEN
335 0 : CPWARN("MOKP: Angular momentum l > 5 not supported in GTO output, skipping")
336 0 : CYCLE
337 : END IF
338 : WRITE (iw, '(A2,I6,F8.2)') &
339 0 : angmom(lshell + 1:lshell + 1)//" ", npgf_set(iset), 1.0_dp
340 : ! Denormalize gcc: undo CP2K's internal normalization
341 : ! (same formula as molden_utils.F, reverse of normalise_gcc_orb)
342 0 : prefac = 2.0_dp**lshell*(2.0_dp/pi)**0.75_dp
343 0 : expzet = 0.25_dp*(2*lshell + 3.0_dp)
344 0 : DO ipgf = 1, npgf_set(iset)
345 : WRITE (iw, fmtstr_gto) &
346 0 : zet(ipgf, iset), &
347 0 : gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet)
348 : END DO
349 : END DO
350 : END DO
351 0 : WRITE (iw, '(A)') ""
352 : END IF
353 : END DO
354 : END IF
355 : END IF
356 :
357 : ! ================================================================
358 : ! WRITE PER-KPOINT DATA
359 : ! ================================================================
360 0 : DO ikp = 1, nkp
361 :
362 : ! ==============================================================
363 : ! A) MO coefficients, eigenvalues, occupations (per spin)
364 : ! ==============================================================
365 0 : DO ispin = 1, nspins
366 :
367 0 : Cre_buf(:, :) = 0.0_dp
368 0 : Cim_buf(:, :) = 0.0_dp
369 0 : eval_buf(:) = 0.0_dp
370 0 : occ_buf(:) = 0.0_dp
371 :
372 0 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
373 0 : ikp_loc = ikp - kp_range(1) + 1
374 0 : kp => kpoints%kp_env(ikp_loc)%kpoint_env
375 0 : mos_kp => kp%mos
376 :
377 : CALL get_mo_set(mos_kp(1, ispin), &
378 : eigenvalues=eigenvalues, &
379 0 : occupation_numbers=occupation_numbers)
380 0 : eval_buf(1:nmo) = eigenvalues(1:nmo)
381 0 : occ_buf(1:nmo) = occupation_numbers(1:nmo)
382 :
383 0 : CALL get_mo_set(mos_kp(1, ispin), mo_coeff=mo_coeff_re)
384 0 : CALL cp_fm_get_submatrix(mo_coeff_re, Cre_buf)
385 :
386 0 : IF (.NOT. use_real_wfn) THEN
387 0 : CALL get_mo_set(mos_kp(2, ispin), mo_coeff=mo_coeff_im)
388 0 : CALL cp_fm_get_submatrix(mo_coeff_im, Cim_buf)
389 : END IF
390 : END IF
391 :
392 0 : CALL para_env_inter_kp%sum(eval_buf)
393 0 : CALL para_env_inter_kp%sum(occ_buf)
394 0 : CALL para_env_inter_kp%sum(Cre_buf)
395 0 : IF (.NOT. use_real_wfn) THEN
396 0 : CALL para_env_inter_kp%sum(Cim_buf)
397 : END IF
398 :
399 0 : IF (iw > 0) THEN
400 0 : WRITE (iw, '(A,2I6)') "# BEGIN_KPOINT_SPIN ikp ispin =", ikp, ispin
401 :
402 0 : WRITE (iw, '(A)') "# EIGENVALUES"
403 0 : WRITE (iw, fmtstr_vec) (eval_buf(n), n=1, nmo)
404 :
405 0 : WRITE (iw, '(A)') "# OCCUPATIONS"
406 0 : WRITE (iw, fmtstr_vec) (occ_buf(n), n=1, nmo)
407 :
408 0 : WRITE (iw, '(A)') "# MO_COEFF_RE (Sparse: i_mo j_ao value)"
409 0 : DO i = 1, nmo
410 0 : DO j = 1, nao
411 0 : IF (ABS(Cre_buf(j, i)) >= thresh) THEN
412 0 : WRITE (iw, fmtstr_sparse) i, j, Cre_buf(j, i)
413 : END IF
414 : END DO
415 : END DO
416 :
417 0 : IF (.NOT. use_real_wfn) THEN
418 0 : WRITE (iw, '(A)') "# MO_COEFF_IM (Sparse: i_mo j_ao value)"
419 0 : DO i = 1, nmo
420 0 : DO j = 1, nao
421 0 : IF (ABS(Cim_buf(j, i)) >= thresh) THEN
422 0 : WRITE (iw, fmtstr_sparse) i, j, Cim_buf(j, i)
423 : END IF
424 : END DO
425 : END DO
426 : END IF
427 0 : WRITE (iw, '(A,2I6)') "# END_KPOINT_SPIN ikp ispin =", ikp, ispin
428 : END IF
429 :
430 : END DO ! ispin
431 :
432 : ! ==============================================================
433 : ! B) Overlap matrix S(k) — only in MATRIX mode
434 : ! ==============================================================
435 0 : IF (write_overlap) THEN
436 :
437 0 : CALL dbcsr_set(rmatrix, 0.0_dp)
438 0 : IF (use_real_wfn) THEN
439 : CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
440 : xkp=xkp(1:3, ikp), cell_to_index=cell_to_index, &
441 0 : sab_nl=sab_nl)
442 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
443 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
444 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
445 0 : Sim_buf(:, :) = 0.0_dp
446 : ELSE
447 0 : CALL dbcsr_set(cmatrix, 0.0_dp)
448 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, &
449 : ispin=1, xkp=xkp(1:3, ikp), &
450 0 : cell_to_index=cell_to_index, sab_nl=sab_nl)
451 0 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
452 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
453 0 : CALL cp_fm_get_submatrix(fm_S_global, Sre_buf)
454 0 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
455 0 : CALL copy_dbcsr_to_fm(tmpmat, fm_S_global)
456 0 : CALL cp_fm_get_submatrix(fm_S_global, Sim_buf)
457 : END IF
458 :
459 0 : IF (iw > 0) THEN
460 0 : WRITE (iw, '(A,I6)') "# BEGIN_OVERLAP ikp =", ikp
461 :
462 0 : WRITE (iw, '(A)') "# OVERLAP_RE (Sparse: i_ao j_ao value)"
463 0 : DO i = 1, nao
464 0 : DO j = 1, nao
465 0 : IF (ABS(Sre_buf(i, j)) >= thresh) THEN
466 0 : WRITE (iw, fmtstr_sparse) i, j, Sre_buf(i, j)
467 : END IF
468 : END DO
469 : END DO
470 0 : IF (.NOT. use_real_wfn) THEN
471 0 : WRITE (iw, '(A)') "# OVERLAP_IM (Sparse: i_ao j_ao value)"
472 0 : DO i = 1, nao
473 0 : DO j = 1, nao
474 0 : IF (ABS(Sim_buf(i, j)) >= thresh) THEN
475 0 : WRITE (iw, fmtstr_sparse) i, j, Sim_buf(i, j)
476 : END IF
477 : END DO
478 : END DO
479 : END IF
480 0 : WRITE (iw, '(A,I6)') "# END_OVERLAP ikp =", ikp
481 : END IF
482 :
483 : END IF ! write_overlap
484 :
485 : END DO ! ikp
486 :
487 : ! ================================================================
488 : ! Finalize
489 : ! ================================================================
490 0 : IF (iw > 0) THEN
491 0 : WRITE (iw, '(A)') "# END_OF_FILE"
492 : END IF
493 :
494 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "")
495 :
496 0 : IF (output_unit > 0) THEN
497 : WRITE (output_unit, '(T3,A)') &
498 0 : "KPOINT_MO_DUMP| Data written to .mokp file"
499 : END IF
500 :
501 : ! Release work arrays
502 0 : IF (write_overlap) THEN
503 0 : CALL cp_fm_release(fm_S_global)
504 0 : CALL dbcsr_deallocate_matrix(rmatrix)
505 0 : CALL dbcsr_deallocate_matrix(cmatrix)
506 0 : CALL dbcsr_deallocate_matrix(tmpmat)
507 0 : DEALLOCATE (Sre_buf, Sim_buf)
508 0 : DEALLOCATE (ao_l)
509 : END IF
510 0 : DEALLOCATE (Cre_buf, Cim_buf)
511 0 : DEALLOCATE (eval_buf, occ_buf)
512 0 : DEALLOCATE (first_sgf, last_sgf)
513 :
514 0 : CALL timestop(handle)
515 :
516 0 : END SUBROUTINE write_kpoint_mo_data
517 :
518 : END MODULE kpoint_mo_dump
|
