(git:7fc1fe1)
Loading...
Searching...
No Matches
trexio_utils.F
Go to the documentation of this file.
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 The module to read/write TREX IO files for interfacing CP2K with other programs
10!> \par History
11!> 05.2024 created [SB]
12!> 05.2026 improved [KN]
13!> \author Stefano Battaglia
14!> \author Kosuke Nakano
15! **************************************************************************************************
16MODULE trexio_utils
17
18 USE ai_onecenter, ONLY: sg_overlap
19 USE atomic_kind_types, ONLY: get_atomic_kind
20 USE basis_set_types, ONLY: gto_basis_set_type, get_gto_basis_set
21 USE cell_types, ONLY: cell_type, pbc, pbc_stable, real_to_scaled
22 USE cp2k_info, ONLY: cp2k_version
23 USE cp_blacs_env, ONLY: cp_blacs_env_type
24 USE cp_control_types, ONLY: dft_control_type
25 USE cp_cfm_types, ONLY: cp_cfm_to_fm
26 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
27 USE cp_files, ONLY: close_file, file_exists, open_file
28 USE cp_fm_types, ONLY: cp_fm_get_info, cp_fm_type, cp_fm_create, cp_fm_set_all, &
29 cp_fm_get_submatrix, cp_fm_to_fm_submat_general, cp_fm_release, &
30 cp_fm_set_element
31 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
32 cp_fm_struct_release, &
33 cp_fm_struct_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_dbcsr_api, ONLY: dbcsr_p_type, dbcsr_iterator_type, dbcsr_iterator_start, &
38 dbcsr_iterator_stop, dbcsr_iterator_blocks_left, &
39 dbcsr_iterator_next_block, dbcsr_copy, dbcsr_set, &
40 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
41 dbcsr_type_symmetric, dbcsr_get_matrix_type
42 USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
43 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
44 USE cp_output_handling, ONLY: medium_print_level
45 USE external_potential_types, ONLY: sgp_potential_type, get_potential
46 USE input_section_types, ONLY: section_vals_type, section_vals_get, &
47 section_vals_val_get
48 USE kinds, ONLY: default_path_length, dp
49 USE kpoint_types, ONLY: kpoint_type
50 USE mathconstants, ONLY: fourpi, pi, fac
51 USE mathlib, ONLY: symmetrize_matrix
52 USE message_passing, ONLY: mp_para_env_type
53 USE orbital_pointers, ONLY: nco, nso
54 USE orbital_transformation_matrices, ONLY: orbtramat
55 USE particle_types, ONLY: particle_type
56 USE qs_energy_types, ONLY: qs_energy_type
57 USE qs_environment_types, ONLY: get_qs_env, &
58 qs_environment_type
59 USE qs_kind_types, ONLY: get_qs_kind, get_qs_kind_set, &
60 qs_kind_type
61 USE qs_mo_types, ONLY: mo_set_type, get_mo_set, init_mo_set, allocate_mo_set
62#ifdef __TREXIO
63 USE kinds, ONLY: default_string_length
64 USE kpoint_methods, ONLY: kpoint_env_initialize, kpoint_init_cell_index, &
65 kpoint_initialize, kpoint_initialize_mo_set, &
66 kpoint_initialize_mos
67 USE kpoint_types, ONLY: kpoint_env_p_type, &
68 get_kpoint_info, get_kpoint_env, kpoint_create, kpoint_release
69 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
70 USE qs_scf_diagonalization, ONLY: do_general_diag_kp
71 USE qs_scf_types, ONLY: qs_scf_env_type
72 USE qs_wannier90, ONLY: prepare_wannier90_scf_mos
73 USE scf_control_types, ONLY: scf_control_type
74 USE trexio, ONLY: trexio_open, trexio_close, &
75 trexio_hdf5, trexio_success, &
76 trexio_string_of_error, trexio_t, trexio_exit_code, &
77 trexio_write_metadata_code, trexio_write_metadata_code_num, &
78 trexio_write_nucleus_coord, trexio_read_nucleus_coord, &
79 trexio_write_nucleus_num, trexio_read_nucleus_num, &
80 trexio_write_nucleus_charge, trexio_read_nucleus_charge, &
81 trexio_write_nucleus_label, trexio_read_nucleus_label, &
82 trexio_write_nucleus_repulsion, &
83 trexio_write_cell_a, trexio_write_cell_b, trexio_write_cell_c, &
84 trexio_write_cell_g_a, trexio_write_cell_g_b, &
85 trexio_write_cell_g_c, trexio_write_cell_two_pi, &
86 trexio_write_pbc_periodic, trexio_write_pbc_k_point_num, &
87 trexio_write_pbc_k_point, trexio_write_pbc_k_point_weight, &
88 trexio_write_electron_num, trexio_read_electron_num, &
89 trexio_write_electron_up_num, trexio_read_electron_up_num, &
90 trexio_write_electron_dn_num, trexio_read_electron_dn_num, &
91 trexio_write_state_num, trexio_write_state_id, &
92 trexio_write_state_energy, &
93 trexio_write_basis_type, trexio_write_basis_prim_num, &
94 trexio_write_basis_shell_num, trexio_read_basis_shell_num, &
95 trexio_write_basis_nucleus_index, &
96 trexio_write_basis_shell_ang_mom, trexio_read_basis_shell_ang_mom, &
97 trexio_write_basis_shell_factor, &
98 trexio_write_basis_r_power, trexio_write_basis_shell_index, &
99 trexio_write_basis_exponent, trexio_write_basis_coefficient, &
100 trexio_write_basis_prim_factor, &
101 trexio_write_ecp_z_core, trexio_write_ecp_max_ang_mom_plus_1, &
102 trexio_write_ecp_num, trexio_write_ecp_ang_mom, &
103 trexio_write_ecp_nucleus_index, trexio_write_ecp_exponent, &
104 trexio_write_ecp_coefficient, trexio_write_ecp_power, &
105 trexio_write_ao_cartesian, trexio_write_ao_num, &
106 trexio_read_ao_cartesian, trexio_read_ao_num, &
107 trexio_write_ao_shell, trexio_write_ao_normalization, &
108 trexio_read_ao_shell, trexio_read_ao_normalization, &
109 trexio_write_mo_num, trexio_write_mo_energy, &
110 trexio_read_mo_num, trexio_read_mo_energy, &
111 trexio_write_mo_occupation, trexio_write_mo_spin, &
112 trexio_read_mo_occupation, trexio_read_mo_spin, &
113 trexio_write_mo_class, trexio_write_mo_coefficient, &
114 trexio_read_mo_class, trexio_read_mo_coefficient, &
115 trexio_write_mo_coefficient_im, trexio_write_mo_k_point, &
116 trexio_write_mo_type
117#endif
118#include "./base/base_uses.f90"
119
120 IMPLICIT NONE
121
122 PRIVATE
123
124 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'trexio_utils'
125
126 PUBLIC :: write_trexio, read_trexio
127
128CONTAINS
129
130! **************************************************************************************************
131!> \brief Write a trexio file
132!> \param qs_env the qs environment with all the info of the computation
133!> \param trexio_section the section with the trexio info
134!> \param energy_derivative ...
135! **************************************************************************************************
136 SUBROUTINE write_trexio(qs_env, trexio_section, energy_derivative)
137 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
138 TYPE(section_vals_type), INTENT(IN), POINTER :: trexio_section
139 TYPE(dbcsr_p_type), INTENT(IN), DIMENSION(:), POINTER, OPTIONAL :: energy_derivative
140
141#ifdef __TREXIO
142 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_trexio'
143
144 INTEGER :: handle, output_unit, unit_trexio
145 CHARACTER(len=default_path_length) :: filename, filename_de
146 INTEGER(trexio_t) :: f ! The TREXIO file handle
147 INTEGER(trexio_exit_code) :: rc ! TREXIO return code
148 LOGICAL :: explicit, do_kpoints, ecp_semi_local, &
149 ecp_local, sgp_potential_present, ionode, &
150 use_real_wfn, save_cartesian, &
151 trexio_kpoints_created
152 REAL(kind=dp) :: e_nn, zeff, expzet, prefac, zeta, gcca, &
153 prim_cart_fac, nsgto
154 TYPE(cell_type), POINTER :: cell
155 TYPE(cp_logger_type), POINTER :: logger
156 TYPE(dft_control_type), POINTER :: dft_control
157 TYPE(gto_basis_set_type), POINTER :: basis_set
158 TYPE(kpoint_type), POINTER :: kpoints, trexio_kpoints
159 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
160 TYPE(qs_energy_type), POINTER :: energy
161 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
162 TYPE(sgp_potential_type), POINTER :: sgp_potential
163 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
164 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_kp
165 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
166 TYPE(mp_para_env_type), POINTER :: para_env, para_env_inter_kp
167 TYPE(cp_blacs_env_type), POINTER :: blacs_env
168 TYPE(cp_fm_struct_type), POINTER :: fm_struct
169 TYPE(cp_fm_type) :: fm_mo_coeff, fm_dummy, fm_mo_coeff_im, &
170 kp_coeff_re, kp_coeff_im
171 TYPE(dbcsr_iterator_type) :: iter
172
173 CHARACTER(LEN=2) :: element_symbol
174 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: label
175 INTEGER :: iatom, natoms, periodic, nkp, nel_tot, &
176 nspins, ikind, ishell_loc, ishell, &
177 shell_num, prim_num, nset, iset, ipgf, z, &
178 sl_lmax, ecp_num, nloc, nsemiloc, sl_l, iecp, &
179 iao, icgf_atom, ncgf, nao_shell, ao_num, nmo, &
180 mo_num, ispin, ikp, imo, ikp_loc, nsgf, ncgf_atom, &
181 i, j, k, l, m, unit_de, &
182 row, col, row_size, col_size, &
183 row_offset, col_offset
184 INTEGER, DIMENSION(2) :: nel_spin, kp_range, nmo_spin
185 INTEGER, DIMENSION(0:10) :: npot
186 INTEGER, DIMENSION(:), ALLOCATABLE :: nucleus_index, shell_ang_mom, r_power, &
187 shell_index, z_core, max_ang_mom_plus_1, &
188 ang_mom, powers, ao_shell, mo_spin, mo_kpoint, &
189 cp2k_to_trexio_ang_mom, ao_to_atom
190 ! Per-atom Bloch-gauge correction:
191 ! CP2K's k-space matrix builder (rskp_transform) Bloch-sums real-space blocks with
192 ! lattice vectors R supplied by the neighbour list. The neighbour list, in turn,
193 ! wraps interatomic vectors through subsys/cell_types.F :: pbc or, for symmetry reduction,
194 ! pbc_stable. Hence the effective
195 ! per-atom gauge is the image shift between the raw and wrapped coordinates, as used by
196 ! kpoint_methods.F :: kpoint_initialize. Derive it through the matching wrapper instead of
197 ! duplicating its boundary convention. Since nucleus_coord is written as the raw particle_set(i)%r, we
198 ! rephase each atom's MO block by
199 ! exp(-i 2*pi * k * agauge_i) so the (coord, MO) pair is self-consistent in the
200 ! standard Bloch convention used by TREXIO consumers.
201 INTEGER, DIMENSION(:, :), ALLOCATABLE :: agauge
202 REAL(kind=dp) :: scoord(3), scoord_pbc(3), r_pbc(3), &
203 kdotg, cval, sval, re_old, im_old
204 INTEGER, DIMENSION(:), POINTER :: nshell, npgf
205 INTEGER, DIMENSION(:, :), POINTER :: l_shell_set
206 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: charge, shell_factor, exponents, coefficients, &
207 prim_factor, ao_normalization, mo_energy, &
208 mo_occupation, sgcc, ecp_coefficients, &
209 sgf_coefficients
210 REAL(kind=dp), DIMENSION(:), POINTER :: wkp, norm_cgf
211 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: coord, mo_coefficient, mo_coefficient_im, &
212 mos_sgf, dedp, sloc
213 REAL(kind=dp), DIMENSION(:, :), POINTER :: zetas, data_block, xkp
214 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
215
216 CALL timeset(routinen, handle)
217
218 NULLIFY (cell, logger, dft_control, basis_set, kpoints, trexio_kpoints, particle_set, &
219 energy, kind_set)
220 NULLIFY (sgp_potential, mos, mos_kp, kp_env, para_env, para_env_inter_kp, blacs_env)
221 NULLIFY (fm_struct, nshell, npgf, l_shell_set, wkp, norm_cgf, zetas, data_block, gcc)
222
223 logger => cp_get_default_logger()
224 output_unit = cp_logger_get_default_io_unit(logger)
225
226 cpassert(ASSOCIATED(qs_env))
227
228 ! get filename
229 CALL section_vals_val_get(trexio_section, "FILENAME", c_val=filename, explicit=explicit)
230 IF (.NOT. explicit) THEN
231 filename = trim(logger%iter_info%project_name)//'-TREXIO.h5'
232 ELSE
233 filename = trim(filename)//'.h5'
234 END IF
235
236 CALL get_qs_env(qs_env, para_env=para_env)
237 ionode = para_env%is_source()
238 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints)
239 trexio_kpoints => kpoints
240 trexio_kpoints_created = .false.
241 CALL prepare_trexio_kpoint_grid(qs_env, trexio_section, do_kpoints, kpoints, &
242 trexio_kpoints, trexio_kpoints_created)
243
244 ! inquire whether a file with the same name already exists, if yes, delete it
245 IF (ionode) THEN
246 IF (file_exists(filename)) THEN
247 CALL open_file(filename, unit_number=unit_trexio)
248 CALL close_file(unit_number=unit_trexio, file_status="DELETE")
249 END IF
250
251 !========================================================================================!
252 ! Open the TREXIO file
253 !========================================================================================!
254 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing trexio file ', trim(filename)
255 f = trexio_open(filename, 'w', trexio_hdf5, rc)
256 CALL trexio_error(rc)
257
258 !========================================================================================!
259 ! Metadata group
260 !========================================================================================!
261 rc = trexio_write_metadata_code_num(f, 1)
262 CALL trexio_error(rc)
263
264 rc = trexio_write_metadata_code(f, cp2k_version, len_trim(cp2k_version) + 1)
265 CALL trexio_error(rc)
266
267 !========================================================================================!
268 ! Nucleus group
269 !========================================================================================!
270 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
271
272 rc = trexio_write_nucleus_num(f, natoms)
273 CALL trexio_error(rc)
274
275 ALLOCATE (coord(3, natoms))
276 ALLOCATE (label(natoms))
277 ALLOCATE (charge(natoms))
278 DO iatom = 1, natoms
279 ! store the coordinates
280 coord(:, iatom) = particle_set(iatom)%r(1:3)
281 ! figure out the element symbol and to which kind_set entry this atomic_kind corresponds to
282 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
283 ! store the element symbol
284 label(iatom) = element_symbol
285 ! get and store the effective nuclear charge of this kind_type (ikind)
286 CALL get_qs_kind(kind_set(ikind), zeff=zeff)
287 charge(iatom) = zeff
288 END DO
289
290 rc = trexio_write_nucleus_coord(f, coord)
291 CALL trexio_error(rc)
292 DEALLOCATE (coord)
293
294 rc = trexio_write_nucleus_charge(f, charge)
295 CALL trexio_error(rc)
296 DEALLOCATE (charge)
297
298 rc = trexio_write_nucleus_label(f, label, 3)
299 CALL trexio_error(rc)
300 DEALLOCATE (label)
301
302 ! nuclear repulsion energy well-defined for molecules only
303 IF (sum(cell%perd) == 0) THEN
304 CALL nuclear_repulsion_energy(particle_set, kind_set, e_nn)
305 rc = trexio_write_nucleus_repulsion(f, e_nn)
306 CALL trexio_error(rc)
307 END IF
308
309 !========================================================================================!
310 ! Cell group
311 !========================================================================================!
312 rc = trexio_write_cell_a(f, cell%hmat(:, 1))
313 CALL trexio_error(rc)
314
315 rc = trexio_write_cell_b(f, cell%hmat(:, 2))
316 CALL trexio_error(rc)
317
318 rc = trexio_write_cell_c(f, cell%hmat(:, 3))
319 CALL trexio_error(rc)
320
321 rc = trexio_write_cell_g_a(f, cell%h_inv(:, 1))
322 CALL trexio_error(rc)
323
324 rc = trexio_write_cell_g_b(f, cell%h_inv(:, 2))
325 CALL trexio_error(rc)
326
327 rc = trexio_write_cell_g_c(f, cell%h_inv(:, 3))
328 CALL trexio_error(rc)
329
330 rc = trexio_write_cell_two_pi(f, 0)
331 CALL trexio_error(rc)
332
333 !========================================================================================!
334 ! PBC group
335 !========================================================================================!
336 periodic = 0
337 IF (sum(cell%perd) /= 0) periodic = 1
338 rc = trexio_write_pbc_periodic(f, periodic)
339 CALL trexio_error(rc)
340
341 IF (do_kpoints) THEN
342 CALL get_kpoint_info(trexio_kpoints, nkp=nkp, xkp=xkp, wkp=wkp)
343
344 rc = trexio_write_pbc_k_point_num(f, nkp)
345 CALL trexio_error(rc)
346
347 rc = trexio_write_pbc_k_point(f, xkp)
348 CALL trexio_error(rc)
349
350 rc = trexio_write_pbc_k_point_weight(f, wkp)
351 CALL trexio_error(rc)
352 END IF
353
354 !========================================================================================!
355 ! Electron group
356 !========================================================================================!
357 CALL get_qs_env(qs_env, dft_control=dft_control, nelectron_total=nel_tot)
358
359 rc = trexio_write_electron_num(f, nel_tot)
360 CALL trexio_error(rc)
361
362 nspins = dft_control%nspins
363 IF (nspins == 1) THEN
364 ! it is a spin-restricted calculation and we need to split the electrons manually,
365 ! because in CP2K they are all otherwise weirdly stored in nelectron_spin(1)
366 nel_spin(1) = nel_tot/2
367 nel_spin(2) = nel_tot/2
368 ELSE
369 ! for UKS/ROKS, the two spin channels are populated correctly and according to
370 ! the multiplicity
371 CALL get_qs_env(qs_env, nelectron_spin=nel_spin)
372 END IF
373 rc = trexio_write_electron_up_num(f, nel_spin(1))
374 CALL trexio_error(rc)
375 rc = trexio_write_electron_dn_num(f, nel_spin(2))
376 CALL trexio_error(rc)
377
378 !========================================================================================!
379 ! State group
380 !========================================================================================!
381 CALL get_qs_env(qs_env, energy=energy)
382
383 rc = trexio_write_state_num(f, 1)
384 CALL trexio_error(rc)
385
386 rc = trexio_write_state_id(f, 1)
387 CALL trexio_error(rc)
388
389 ! rc = trexio_write_state_energy(f, energy%total)
390 CALL trexio_error(rc)
391
392 END IF ! ionode
393
394 !========================================================================================!
395 ! Basis group
396 !========================================================================================!
397 CALL get_qs_env(qs_env, qs_kind_set=kind_set, natom=natoms, particle_set=particle_set)
398 CALL get_qs_kind_set(kind_set, nshell=shell_num, npgf_seg=prim_num)
399
400 CALL section_vals_val_get(trexio_section, "CARTESIAN", l_val=save_cartesian)
401
402 IF (ionode) THEN
403 rc = trexio_write_basis_type(f, 'Gaussian', len_trim('Gaussian') + 1)
404 CALL trexio_error(rc)
405
406 rc = trexio_write_basis_shell_num(f, shell_num)
407 CALL trexio_error(rc)
408
409 rc = trexio_write_basis_prim_num(f, prim_num)
410 CALL trexio_error(rc)
411 END IF ! ionode
412
413 ! one-to-one mapping between shells and ...
414 ALLOCATE (nucleus_index(shell_num)) ! ...atomic indices
415 ALLOCATE (shell_ang_mom(shell_num)) ! ...angular momenta
416 ALLOCATE (shell_index(prim_num)) ! ...indices of primitive functions
417 ALLOCATE (exponents(prim_num)) ! ...primitive exponents
418 ALLOCATE (coefficients(prim_num)) ! ...contraction coefficients
419 ALLOCATE (prim_factor(prim_num)) ! ...primitive normalization factors
420
421 ! needed in AO group
422 IF (.NOT. save_cartesian) THEN
423 ALLOCATE (sgf_coefficients(prim_num)) ! ...contraction coefficients
424 END IF
425
426 ishell = 0 ! global shell index
427 ipgf = 0 ! global primitives index
428 DO iatom = 1, natoms
429 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
430 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
431 ! get the primary (orbital) basis set associated to this qs_kind
432 CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
433 ! get the info from the basis set
434 CALL get_gto_basis_set(basis_set, &
435 nset=nset, &
436 nshell=nshell, &
437 npgf=npgf, &
438 zet=zetas, &
439 gcc=gcc, &
440 l=l_shell_set)
441
442 DO iset = 1, nset
443 DO ishell_loc = 1, nshell(iset)
444 ishell = ishell + 1
445
446 ! nucleus_index array
447 nucleus_index(ishell) = iatom
448
449 ! shell_ang_mom array
450 l = l_shell_set(ishell_loc, iset)
451 shell_ang_mom(ishell) = l
452
453 ! shell_index array
454 shell_index(ipgf + 1:ipgf + npgf(iset)) = ishell
455
456 ! exponents array
457 exponents(ipgf + 1:ipgf + npgf(iset)) = zetas(1:npgf(iset), iset)
458
459 ! compute on the fly the normalization factor as in normalise_gcc_orb
460 ! and recover the original contraction coefficients to store them separately
461 expzet = 0.25_dp*real(2*l + 3, dp)
462 prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
463 DO i = 1, npgf(iset)
464 gcca = gcc(i, ishell_loc, iset)
465 zeta = zetas(i, iset)
466 prim_cart_fac = prefac*zeta**expzet
467
468 ! contraction coefficients array
469 coefficients(ipgf + i) = gcca/prim_cart_fac
470
471 IF (save_cartesian) THEN
472 ! primitives normalization factors array
473 prim_factor(ipgf + i) = prim_cart_fac
474 ELSE
475 ! for spherical harmonics we have a different factor
476 prim_factor(ipgf + i) = sgf_norm(l, exponents(ipgf + i))
477 ! we need these later in the AO group
478 sgf_coefficients(ipgf + i) = coefficients(ipgf + i)*prim_factor(ipgf + i)
479 END IF
480
481 END DO
482
483 ipgf = ipgf + npgf(iset)
484 END DO
485 END DO
486 END DO
487 ! just a failsafe check
488 cpassert(ishell == shell_num)
489 cpassert(ipgf == prim_num)
490
491 IF (ionode) THEN
492 rc = trexio_write_basis_nucleus_index(f, nucleus_index)
493 CALL trexio_error(rc)
494
495 rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
496 CALL trexio_error(rc)
497
498 ! Normalization factors are shoved in the AO group
499 ALLOCATE (shell_factor(shell_num)) ! 1-to-1 map bw shells and normalization factors
500 shell_factor(:) = 1.0_dp
501 rc = trexio_write_basis_shell_factor(f, shell_factor)
502 CALL trexio_error(rc)
503 DEALLOCATE (shell_factor)
504
505 ! This is always 0 for Gaussian basis sets
506 ALLOCATE (r_power(shell_num)) ! 1-to-1 map bw shells radial function powers
507 r_power(:) = 0
508 rc = trexio_write_basis_r_power(f, r_power)
509 CALL trexio_error(rc)
510 DEALLOCATE (r_power)
511
512 rc = trexio_write_basis_shell_index(f, shell_index)
513 CALL trexio_error(rc)
514
515 rc = trexio_write_basis_exponent(f, exponents)
516 CALL trexio_error(rc)
517
518 rc = trexio_write_basis_coefficient(f, coefficients)
519 CALL trexio_error(rc)
520
521 ! Normalization factors are shoved in the AO group
522 rc = trexio_write_basis_prim_factor(f, prim_factor)
523 CALL trexio_error(rc)
524 END IF
525
526 DEALLOCATE (nucleus_index)
527 DEALLOCATE (shell_index)
528 DEALLOCATE (exponents)
529 DEALLOCATE (coefficients)
530 DEALLOCATE (prim_factor)
531 ! shell_ang_mom is needed in the MO group, so will be deallocated there
532
533 !========================================================================================!
534 ! ECP group
535 !========================================================================================!
536 IF (ionode) THEN
537 CALL get_qs_kind_set(kind_set, sgp_potential_present=sgp_potential_present)
538
539 ! figure out whether we actually have ECP potentials
540 ecp_num = 0
541 IF (sgp_potential_present) THEN
542 DO iatom = 1, natoms
543 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
544 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
545 ! get the the sgp_potential associated to this qs_kind
546 CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential)
547
548 ! get the info on the potential
549 IF (ASSOCIATED(sgp_potential)) THEN
550 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
551 IF (ecp_local) THEN
552 ! get number of local terms
553 CALL get_potential(potential=sgp_potential, nloc=nloc)
554 ecp_num = ecp_num + nloc
555 END IF
556 IF (ecp_semi_local) THEN
557 ! get number of semilocal terms
558 CALL get_potential(potential=sgp_potential, npot=npot)
559 ecp_num = ecp_num + sum(npot)
560 END IF
561 END IF
562 END DO
563 END IF
564
565 ! if we have ECP potentials, populate the ECP group
566 IF (ecp_num > 0) THEN
567 ALLOCATE (z_core(natoms))
568 ALLOCATE (max_ang_mom_plus_1(natoms))
569 max_ang_mom_plus_1(:) = 0
570
571 ALLOCATE (ang_mom(ecp_num))
572 ALLOCATE (nucleus_index(ecp_num))
573 ALLOCATE (exponents(ecp_num))
574 ALLOCATE (ecp_coefficients(ecp_num))
575 ALLOCATE (powers(ecp_num))
576
577 iecp = 0
578 DO iatom = 1, natoms
579 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
580 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind, z=z)
581 ! get the the sgp_potential associated to this qs_kind
582 CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential, zeff=zeff)
583
584 ! number of core electrons removed by the ECP
585 z_core(iatom) = z - int(zeff)
586
587 ! get the info on the potential
588 IF (ASSOCIATED(sgp_potential)) THEN
589 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
590
591 ! deal with the local part
592 IF (ecp_local) THEN
593 CALL get_potential(potential=sgp_potential, nloc=nloc, sl_lmax=sl_lmax)
594 ang_mom(iecp + 1:iecp + nloc) = sl_lmax + 1
595 nucleus_index(iecp + 1:iecp + nloc) = iatom
596 exponents(iecp + 1:iecp + nloc) = sgp_potential%bloc(1:nloc)
597 ecp_coefficients(iecp + 1:iecp + nloc) = sgp_potential%aloc(1:nloc)
598 powers(iecp + 1:iecp + nloc) = sgp_potential%nrloc(1:nloc) - 2
599 iecp = iecp + nloc
600 END IF
601
602 ! deal with the semilocal part
603 IF (ecp_semi_local) THEN
604 CALL get_potential(potential=sgp_potential, npot=npot, sl_lmax=sl_lmax)
605 max_ang_mom_plus_1(iatom) = sl_lmax + 1
606
607 DO sl_l = 0, sl_lmax
608 nsemiloc = npot(sl_l)
609 ang_mom(iecp + 1:iecp + nsemiloc) = sl_l
610 nucleus_index(iecp + 1:iecp + nsemiloc) = iatom
611 exponents(iecp + 1:iecp + nsemiloc) = sgp_potential%bpot(1:nsemiloc, sl_l)
612 ecp_coefficients(iecp + 1:iecp + nsemiloc) = sgp_potential%apot(1:nsemiloc, sl_l)
613 powers(iecp + 1:iecp + nsemiloc) = sgp_potential%nrpot(1:nsemiloc, sl_l) - 2
614 iecp = iecp + nsemiloc
615 END DO
616 END IF
617 END IF
618 END DO
619
620 ! fail-safe check
621 cpassert(iecp == ecp_num)
622
623 rc = trexio_write_ecp_num(f, ecp_num)
624 CALL trexio_error(rc)
625
626 rc = trexio_write_ecp_z_core(f, z_core)
627 CALL trexio_error(rc)
628 DEALLOCATE (z_core)
629
630 rc = trexio_write_ecp_max_ang_mom_plus_1(f, max_ang_mom_plus_1)
631 CALL trexio_error(rc)
632 DEALLOCATE (max_ang_mom_plus_1)
633
634 rc = trexio_write_ecp_ang_mom(f, ang_mom)
635 CALL trexio_error(rc)
636 DEALLOCATE (ang_mom)
637
638 rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
639 CALL trexio_error(rc)
640 DEALLOCATE (nucleus_index)
641
642 rc = trexio_write_ecp_exponent(f, exponents)
643 CALL trexio_error(rc)
644 DEALLOCATE (exponents)
645
646 rc = trexio_write_ecp_coefficient(f, ecp_coefficients)
647 CALL trexio_error(rc)
648 DEALLOCATE (ecp_coefficients)
649
650 rc = trexio_write_ecp_power(f, powers)
651 CALL trexio_error(rc)
652 DEALLOCATE (powers)
653 END IF
654
655 END IF ! ionode
656
657 !========================================================================================!
658 ! Grid group
659 !========================================================================================!
660 ! TODO
661
662 !========================================================================================!
663 ! AO group
664 !========================================================================================!
665 CALL get_qs_env(qs_env, qs_kind_set=kind_set)
666 CALL get_qs_kind_set(kind_set, ncgf=ncgf, nsgf=nsgf)
667
668 IF (save_cartesian) THEN
669 ao_num = ncgf
670 ELSE
671 ao_num = nsgf
672 END IF
673
674 IF (ionode) THEN
675 IF (save_cartesian) THEN
676 rc = trexio_write_ao_cartesian(f, 1)
677 ELSE
678 rc = trexio_write_ao_cartesian(f, 0)
679 END IF
680 CALL trexio_error(rc)
681
682 rc = trexio_write_ao_num(f, ao_num)
683 CALL trexio_error(rc)
684 END IF
685
686 ! one-to-one mapping between AOs and ...
687 ALLOCATE (ao_shell(ao_num)) ! ..shells
688 ALLOCATE (ao_normalization(ao_num)) ! ..normalization factors
689 ALLOCATE (ao_to_atom(ao_num)) ! ..parent atom (needed for the k-point gauge fix)
690
691 IF (.NOT. save_cartesian) THEN
692 ! AO order map from CP2K to TREXIO convention
693 ! from m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
694 ! to m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
695 ALLOCATE (cp2k_to_trexio_ang_mom(ao_num))
696 i = 0
697 DO ishell = 1, shell_num
698 l = shell_ang_mom(ishell)
699 DO k = 1, 2*l + 1
700 m = (-1)**k*floor(real(k, kind=dp)/2.0_dp)
701 cp2k_to_trexio_ang_mom(i + k) = i + l + 1 + m
702 END DO
703 i = i + 2*l + 1
704 END DO
705 cpassert(i == ao_num)
706 END IF
707
708 ! we need to be consistent with the basis group on the shell indices
709 ishell = 0 ! global shell index
710 iao = 0 ! global AO index
711 ipgf = 0 ! global primitives index
712 DO iatom = 1, natoms
713 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
714 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
715 ! get the primary (orbital) basis set associated to this qs_kind
716 CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
717 ! get the info from the basis set
718 CALL get_gto_basis_set(basis_set, &
719 nset=nset, &
720 nshell=nshell, &
721 norm_cgf=norm_cgf, &
722 ncgf=ncgf_atom, &
723 npgf=npgf, &
724 zet=zetas, &
725 l=l_shell_set)
726
727 icgf_atom = 0
728 DO iset = 1, nset
729 DO ishell_loc = 1, nshell(iset)
730 ! global shell index
731 ishell = ishell + 1
732 ! angular momentum l of this shell
733 l = l_shell_set(ishell_loc, iset)
734
735 ! number of AOs in this shell
736 IF (save_cartesian) THEN
737 nao_shell = nco(l)
738 ELSE
739 nao_shell = nso(l)
740 END IF
741
742 ! one-to-one mapping between AOs and shells
743 ao_shell(iao + 1:iao + nao_shell) = ishell
744
745 ! one-to-one mapping between AOs and parent atoms
746 ao_to_atom(iao + 1:iao + nao_shell) = iatom
747
748 ! one-to-one mapping between AOs and normalization factors
749 IF (save_cartesian) THEN
750 ao_normalization(iao + 1:iao + nao_shell) = norm_cgf(icgf_atom + 1:icgf_atom + nao_shell)
751 ELSE
752 ! for each shell, compute the overlap between spherical primitives
753 ALLOCATE (sloc(npgf(iset), npgf(iset)))
754 ALLOCATE (sgcc(npgf(iset)))
755 CALL sg_overlap(sloc, l, zetas(1:npgf(iset), iset), zetas(1:npgf(iset), iset))
756
757 ! and compute the normalizaztion factor for contracted spherical GTOs
758 sgcc(:) = matmul(sloc, sgf_coefficients(ipgf + 1:ipgf + npgf(iset)))
759 nsgto = 1.0_dp/sqrt(dot_product(sgf_coefficients(ipgf + 1:ipgf + npgf(iset)), sgcc))
760
761 DEALLOCATE (sloc)
762 DEALLOCATE (sgcc)
763
764 ! TREXIO employs solid harmonics and not spherical harmonics like cp2k
765 ! so we need the opposite of Racah normalization, multiplied by the Nsgto
766 ! just computed above
767 ao_normalization(iao + 1:iao + nao_shell) = nsgto*sqrt((2*l + 1)/(4*pi))
768 END IF
769
770 ipgf = ipgf + npgf(iset)
771 iao = iao + nao_shell
772 icgf_atom = icgf_atom + nco(l)
773 END DO
774 END DO
775 ! just a failsafe check
776 cpassert(icgf_atom == ncgf_atom)
777 END DO
778
779 IF (ionode) THEN
780 rc = trexio_write_ao_shell(f, ao_shell)
781 CALL trexio_error(rc)
782
783 rc = trexio_write_ao_normalization(f, ao_normalization)
784 CALL trexio_error(rc)
785 END IF
786
787 DEALLOCATE (ao_shell)
788 DEALLOCATE (ao_normalization)
789 IF (ALLOCATED(sgf_coefficients)) DEALLOCATE (sgf_coefficients)
790
791 !========================================================================================!
792 ! MO group
793 !========================================================================================!
794 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints, dft_control=dft_control, &
795 particle_set=particle_set, qs_kind_set=kind_set, blacs_env=blacs_env, &
796 cell=cell)
797 nspins = dft_control%nspins
798 CALL get_qs_kind_set(kind_set, nsgf=nsgf, ncgf=ncgf)
799 nmo_spin = 0
800
801 ! figure out that total number of MOs
802 mo_num = 0
803 IF (do_kpoints) THEN
804 CALL get_kpoint_info(trexio_kpoints, kp_env=kp_env, nkp=nkp, use_real_wfn=use_real_wfn)
805 CALL get_kpoint_env(kp_env(1)%kpoint_env, mos=mos_kp)
806 DO ispin = 1, nspins
807 CALL get_mo_set(mos_kp(ispin), nmo=nmo)
808 nmo_spin(ispin) = nmo
809 END DO
810 mo_num = nkp*sum(nmo_spin)
811
812 ! we create a distributed fm matrix to gather the MOs from everywhere (in sph basis)
813 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
814 nrow_global=nsgf, ncol_global=mo_num)
815 CALL cp_fm_create(fm_mo_coeff, fm_struct)
816 CALL cp_fm_set_all(fm_mo_coeff, 0.0_dp)
817 IF (.NOT. use_real_wfn) THEN
818 CALL cp_fm_create(fm_mo_coeff_im, fm_struct)
819 CALL cp_fm_set_all(fm_mo_coeff_im, 0.0_dp)
820 CALL cp_fm_create(kp_coeff_re, mos_kp(1)%cmo_coeff%matrix_struct)
821 CALL cp_fm_create(kp_coeff_im, mos_kp(1)%cmo_coeff%matrix_struct)
822 END IF
823 CALL cp_fm_struct_release(fm_struct)
824 ELSE
825 CALL get_qs_env(qs_env, mos=mos)
826 DO ispin = 1, nspins
827 CALL get_mo_set(mos(ispin), nmo=nmo)
828 nmo_spin(ispin) = nmo
829 END DO
830 mo_num = sum(nmo_spin)
831 END IF
832
833 ! allocate all the arrays
834 ALLOCATE (mo_coefficient(ao_num, mo_num))
835 mo_coefficient(:, :) = 0.0_dp
836 ALLOCATE (mo_energy(mo_num))
837 mo_energy(:) = 0.0_dp
838 ALLOCATE (mo_occupation(mo_num))
839 mo_occupation(:) = 0.0_dp
840 ALLOCATE (mo_spin(mo_num))
841 mo_spin(:) = 0
842 ! extra arrays for kpoints
843 IF (do_kpoints) THEN
844 ALLOCATE (mo_coefficient_im(ao_num, mo_num))
845 mo_coefficient_im(:, :) = 0.0_dp
846 ALLOCATE (mo_kpoint(mo_num))
847 mo_kpoint(:) = 0
848 END IF
849
850 ! in case of kpoints, we do this in 2 steps:
851 ! 1. we gather the MOs of each kpt and pipe them into a single large distributed fm matrix;
852 ! 2. we possibly transform the MOs of each kpt to Cartesian AOs and write them in the single large local array;
853 IF (do_kpoints) THEN
854 CALL get_kpoint_info(trexio_kpoints, kp_env=kp_env, nkp=nkp, kp_range=kp_range)
855
856 DO ispin = 1, nspins
857 DO ikp = 1, nkp
858 nmo = nmo_spin(ispin)
859 ! global index to store the MOs
860 imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
861
862 ! do I have this kpoint on this rank?
863 IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
864 ikp_loc = ikp - kp_range(1) + 1
865 ! get the mo set for this kpoint
866 CALL get_kpoint_env(kp_env(ikp_loc)%kpoint_env, mos=mos_kp)
867
868 IF (use_real_wfn) THEN
869 IF (mos_kp(ispin)%use_mo_coeff_b) THEN
870 CALL copy_dbcsr_to_fm(mos_kp(ispin)%mo_coeff_b, mos_kp(ispin)%mo_coeff)
871 END IF
872 ELSE
873 CALL cp_cfm_to_fm(mos_kp(ispin)%cmo_coeff, kp_coeff_re, kp_coeff_im)
874 END IF
875 IF (use_real_wfn) THEN
876 CALL cp_fm_to_fm_submat_general(mos_kp(ispin)%mo_coeff, fm_mo_coeff, &
877 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
878 ELSE
879 CALL cp_fm_to_fm_submat_general(kp_coeff_re, fm_mo_coeff, &
880 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
881 END IF
882
883 ! copy MO energies to local arrays
884 mo_energy(imo + 1:imo + nmo) = mos_kp(ispin)%eigenvalues(1:nmo)
885
886 ! copy MO occupations to local arrays
887 mo_occupation(imo + 1:imo + nmo) = mos_kp(ispin)%occupation_numbers(1:nmo)
888
889 ! same for the imaginary part of MO coefficients
890 IF (.NOT. use_real_wfn) THEN
891 CALL cp_fm_to_fm_submat_general(kp_coeff_im, fm_mo_coeff_im, &
892 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
893 END IF
894 ELSE
895 ! call with a dummy fm for receiving the data
896 CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff, &
897 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
898 IF (.NOT. use_real_wfn) THEN
899 CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff_im, &
900 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
901 END IF
902 END IF
903 END DO
904 END DO
905 END IF
906
907 ! reduce MO energies and occupations to the master node
908 IF (do_kpoints) THEN
909 CALL get_kpoint_info(trexio_kpoints, para_env_inter_kp=para_env_inter_kp)
910 CALL para_env_inter_kp%sum(mo_energy)
911 CALL para_env_inter_kp%sum(mo_occupation)
912 END IF
913
914 ! Bloch-gauge correction (k-points, complex wfn only):
915 ! Build per-atom agauge matching kpoint_methods.F :: kpoint_initialize. The MO
916 ! coefficients gathered above are referenced to atoms wrapped into [-1/2, 1/2),
917 ! while nucleus_coord was written using the raw particle_set(i)%r. We compensate
918 ! by multiplying each AO column block by exp(-i 2*pi * k * agauge_i) per k-point.
919 IF (do_kpoints .AND. .NOT. use_real_wfn) THEN
920 CALL get_kpoint_info(trexio_kpoints, xkp=xkp)
921 ALLOCATE (agauge(3, natoms))
922 DO iatom = 1, natoms
923 CALL real_to_scaled(scoord, particle_set(iatom)%r(1:3), cell)
924 IF (kpoints%symmetry) THEN
925 r_pbc = pbc_stable(particle_set(iatom)%r(1:3), cell)
926 ELSE
927 r_pbc = pbc(particle_set(iatom)%r(1:3), cell)
928 END IF
929 CALL real_to_scaled(scoord_pbc, r_pbc, cell)
930 agauge(:, iatom) = nint(scoord_pbc - scoord)
931 END DO
932 END IF
933
934 ! second step: here we actually put everything in the local arrays for writing to trexio
935 DO ispin = 1, nspins
936 ! get number of MOs for this spin
937 nmo = nmo_spin(ispin)
938 ! allocate local temp array to transform the MOs of each kpoint/spin
939 ALLOCATE (mos_sgf(nsgf, nmo))
940 mos_sgf(:, :) = 0.0_dp
941
942 IF (do_kpoints) THEN
943 DO ikp = 1, nkp
944 ! global index to store the MOs
945 imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
946
947 ! store kpoint index
948 mo_kpoint(imo + 1:imo + nmo) = ikp
949 ! store the MO spins
950 mo_spin(imo + 1:imo + nmo) = ispin - 1
951
952 ! transform and store the MO coefficients
953 CALL cp_fm_get_submatrix(fm_mo_coeff, mos_sgf, 1, imo + 1, nsgf, nmo)
954 IF (save_cartesian) THEN
955 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
956 ELSE
957 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
958 DO i = 1, nsgf
959 mo_coefficient(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
960 END DO
961 END IF
962
963 ! we have to do it for the imaginary part as well
964 IF (.NOT. use_real_wfn) THEN
965 CALL cp_fm_get_submatrix(fm_mo_coeff_im, mos_sgf, 1, imo + 1, nsgf, nmo)
966 IF (save_cartesian) THEN
967 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient_im(:, imo + 1:imo + nmo))
968 ELSE
969 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
970 DO i = 1, nsgf
971 mo_coefficient_im(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
972 END DO
973 END IF
974
975 ! Apply per-atom Bloch-gauge phase factor exp(-i 2*pi * k_ikp * agauge_iatom)
976 ! to remove the spurious phase the consumer would otherwise pick up from
977 ! the (raw nucleus_coord, agauge-gauge MO) mismatch.
978 DO iao = 1, ao_num
979 iatom = ao_to_atom(iao)
980 kdotg = 2.0_dp*pi*dot_product(xkp(:, ikp), real(agauge(:, iatom), kind=dp))
981 cval = cos(kdotg)
982 sval = sin(kdotg)
983 DO j = imo + 1, imo + nmo
984 re_old = mo_coefficient(iao, j)
985 im_old = mo_coefficient_im(iao, j)
986 mo_coefficient(iao, j) = cval*re_old + sval*im_old
987 mo_coefficient_im(iao, j) = -sval*re_old + cval*im_old
988 END DO
989 END DO
990 END IF
991 END DO
992 ELSE ! no k-points
993 ! global index to store the MOs
994 imo = (ispin - 1)*nmo_spin(1)
995 ! store the MO energies
996 mo_energy(imo + 1:imo + nmo) = mos(ispin)%eigenvalues
997 ! store the MO occupations
998 mo_occupation(imo + 1:imo + nmo) = mos(ispin)%occupation_numbers
999 ! store the MO spins
1000 mo_spin(imo + 1:imo + nmo) = ispin - 1
1001
1002 ! check if we are using the dbcsr mo_coeff and copy them to fm if needed
1003 IF (mos(ispin)%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff)
1004
1005 ! allocate a normal fortran array to store the spherical MO coefficients
1006 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, mos_sgf)
1007
1008 IF (save_cartesian) THEN
1009 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
1010 ELSE
1011 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
1012 DO i = 1, nsgf
1013 mo_coefficient(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
1014 END DO
1015 END IF
1016 END IF
1017
1018 DEALLOCATE (mos_sgf)
1019 END DO
1020
1021 IF (ionode) THEN
1022 rc = trexio_write_mo_type(f, 'Canonical', len_trim('Canonical') + 1)
1023 CALL trexio_error(rc)
1024
1025 rc = trexio_write_mo_num(f, mo_num)
1026 CALL trexio_error(rc)
1027
1028 rc = trexio_write_mo_coefficient(f, mo_coefficient)
1029 CALL trexio_error(rc)
1030
1031 rc = trexio_write_mo_energy(f, mo_energy)
1032 CALL trexio_error(rc)
1033
1034 rc = trexio_write_mo_occupation(f, mo_occupation)
1035 CALL trexio_error(rc)
1036
1037 rc = trexio_write_mo_spin(f, mo_spin)
1038 CALL trexio_error(rc)
1039
1040 IF (do_kpoints) THEN
1041 rc = trexio_write_mo_coefficient_im(f, mo_coefficient_im)
1042 CALL trexio_error(rc)
1043
1044 rc = trexio_write_mo_k_point(f, mo_kpoint)
1045 CALL trexio_error(rc)
1046 END IF
1047 END IF
1048
1049 DEALLOCATE (mo_coefficient)
1050 DEALLOCATE (mo_energy)
1051 DEALLOCATE (mo_occupation)
1052 DEALLOCATE (mo_spin)
1053 IF (do_kpoints) THEN
1054 DEALLOCATE (mo_coefficient_im)
1055 DEALLOCATE (mo_kpoint)
1056 CALL cp_fm_release(fm_mo_coeff)
1057 IF (.NOT. use_real_wfn) THEN
1058 CALL cp_fm_release(fm_mo_coeff_im)
1059 CALL cp_fm_release(kp_coeff_re)
1060 CALL cp_fm_release(kp_coeff_im)
1061 END IF
1062 END IF
1063 IF (ALLOCATED(ao_to_atom)) DEALLOCATE (ao_to_atom)
1064 IF (ALLOCATED(agauge)) DEALLOCATE (agauge)
1065 IF (trexio_kpoints_created) CALL kpoint_release(trexio_kpoints)
1066
1067 !========================================================================================!
1068 ! RDM group
1069 !========================================================================================!
1070 !TODO
1071
1072 !========================================================================================!
1073 ! Energy derivative group
1074 !========================================================================================!
1075 IF (PRESENT(energy_derivative)) THEN
1076 filename_de = trim(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
1077
1078 ALLOCATE (dedp(nsgf, nsgf))
1079 dedp(:, :) = 0.0_dp
1080
1081 DO ispin = 1, nspins
1082 CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
1083 DO WHILE (dbcsr_iterator_blocks_left(iter))
1084 ! the offsets tell me the global index of the matrix, not the index of the block
1085 CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
1086 row_size=row_size, col_size=col_size, &
1087 row_offset=row_offset, col_offset=col_offset)
1088
1089 ! Copy data from block to array
1090 DO i = 1, row_size
1091 DO j = 1, col_size
1092 dedp(row_offset + i - 1, col_offset + j - 1) = data_block(i, j)
1093 END DO
1094 END DO
1095 END DO
1096 CALL dbcsr_iterator_stop(iter)
1097
1098 ! symmetrize the matrix if needed
1099 SELECT CASE (dbcsr_get_matrix_type(energy_derivative(ispin)%matrix))
1100 CASE (dbcsr_type_symmetric)
1101 CALL symmetrize_matrix(dedp, "upper_to_lower")
1102 CASE (dbcsr_type_antisymmetric)
1103 CALL symmetrize_matrix(dedp, "anti_upper_to_lower")
1104 CASE (dbcsr_type_no_symmetry)
1105 CASE DEFAULT
1106 cpabort("Unknown matrix type for energy derivative")
1107 END SELECT
1108 END DO
1109
1110 ! reduce the dEdP matrix to the master node
1111 CALL para_env%sum(dedp)
1112
1113 ! print the dEdP matrix to a file
1114 IF (ionode) THEN
1115 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing derivative file ', trim(filename_de)
1116
1117 unit_de = 10
1118 CALL open_file(file_name=filename_de, &
1119 file_action="WRITE", &
1120 file_status="UNKNOWN", &
1121 unit_number=unit_de)
1122 WRITE (unit_de, '(I0, 1X, I0)') nsgf, nsgf
1123 DO i = 1, nsgf
1124 WRITE (unit_de, '(*(1X, F15.8))') (dedp(cp2k_to_trexio_ang_mom(i), &
1125 cp2k_to_trexio_ang_mom(j)), j=1, nsgf)
1126 END DO
1127 CALL close_file(unit_number=unit_de)
1128 END IF
1129
1130 DEALLOCATE (dedp)
1131 END IF
1132
1133 ! Deallocate arrays used throughout the subroutine
1134 IF (ALLOCATED(shell_ang_mom)) DEALLOCATE (shell_ang_mom)
1135 IF (ALLOCATED(cp2k_to_trexio_ang_mom)) DEALLOCATE (cp2k_to_trexio_ang_mom)
1136
1137 !========================================================================================!
1138 ! Close the TREXIO file
1139 !========================================================================================!
1140 IF (ionode) THEN
1141 rc = trexio_close(f)
1142 CALL trexio_error(rc)
1143 END IF
1144
1145 CALL timestop(handle)
1146#else
1147 mark_used(qs_env)
1148 mark_used(trexio_section)
1149 mark_used(energy_derivative)
1150 cpwarn('TREXIO support has not been enabled in this build.')
1151#endif
1152
1153 END SUBROUTINE write_trexio
1154
1155! **************************************************************************************************
1156!> \brief Prepare the k-point object used for TREXIO export.
1157!> \param qs_env the QS environment
1158!> \param trexio_section the TREXIO print section
1159!> \param do_kpoints true when the SCF used k-points
1160!> \param kpoints_scf the converged SCF k-point object
1161!> \param kpoints_out the k-point object to write
1162!> \param created true if kpoints_out must be released by the caller
1163! **************************************************************************************************
1164 SUBROUTINE prepare_trexio_kpoint_grid(qs_env, trexio_section, do_kpoints, kpoints_scf, &
1165 kpoints_out, created)
1166 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1167 TYPE(section_vals_type), INTENT(IN), POINTER :: trexio_section
1168 LOGICAL, INTENT(IN) :: do_kpoints
1169 TYPE(kpoint_type), POINTER :: kpoints_scf, kpoints_out
1170 LOGICAL, INTENT(OUT) :: created
1171
1172#ifdef __TREXIO
1173 CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_trexio_kpoint_grid'
1174
1175 CHARACTER(LEN=default_string_length) :: kp_scheme, reuse_reason
1176 INTEGER :: aligned_blocks, aligned_max_size, handle, &
1177 nfull, output_unit
1178 INTEGER, DIMENSION(3) :: nkp_grid
1179 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1180 LOGICAL :: diis_step, full_grid, full_kpoint_grid, &
1181 gamma_centered, reuse_scf_mos, &
1182 reused_scf_mos, symmetry
1183 REAL(kind=dp) :: aligned_min_svalue, eps_geo, wsum
1184 REAL(kind=dp), DIMENSION(3) :: kp_shift
1185 REAL(kind=dp), DIMENSION(:), POINTER :: wkp_source
1186 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp_source
1187 TYPE(cell_type), POINTER :: cell
1188 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1189 TYPE(cp_logger_type), POINTER :: logger
1190 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
1191 TYPE(dft_control_type), POINTER :: dft_control
1192 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1193 TYPE(mp_para_env_type), POINTER :: para_env
1194 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1195 POINTER :: sab_nl
1196 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1197 TYPE(qs_scf_env_type), POINTER :: scf_env
1198 TYPE(scf_control_type), POINTER :: scf_control
1199
1200 CALL timeset(routinen, handle)
1201
1202 created = .false.
1203 kpoints_out => kpoints_scf
1204 NULLIFY (blacs_env, cell, cell_to_index, dft_control, logger, matrix_ks, matrix_s, mos, &
1205 para_env, particle_set, sab_nl, scf_control, scf_env, wkp_source, xkp_source)
1206
1207 CALL section_vals_val_get(trexio_section, "FULL_KPOINT_GRID", l_val=full_kpoint_grid)
1208 IF (.NOT. do_kpoints .OR. .NOT. full_kpoint_grid) THEN
1209 CALL timestop(handle)
1210 RETURN
1211 END IF
1212 cpassert(ASSOCIATED(kpoints_scf))
1213
1214 CALL get_kpoint_info(kpoints_scf, kp_scheme=kp_scheme, symmetry=symmetry, &
1215 full_grid=full_grid, nkp_grid=nkp_grid, kp_shift=kp_shift, &
1216 gamma_centered=gamma_centered, eps_geo=eps_geo)
1217 IF (.NOT. symmetry .OR. full_grid) THEN
1218 CALL timestop(handle)
1219 RETURN
1220 END IF
1221
1222 SELECT CASE (trim(kp_scheme))
1223 CASE ("MONKHORST-PACK", "MACDONALD", "GENERAL")
1224 ! supported below
1225 CASE DEFAULT
1226 cpabort("TREXIO%FULL_KPOINT_GRID supports only MONKHORST-PACK, MACDONALD, and GENERAL k-points.")
1227 END SELECT
1228
1229 logger => cp_get_default_logger()
1230 output_unit = cp_logger_get_default_io_unit(logger)
1231 CALL section_vals_val_get(trexio_section, "REUSE_SCF_MOS", l_val=reuse_scf_mos)
1232 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, cell=cell, &
1233 particle_set=particle_set, mos=mos, dft_control=dft_control, &
1234 sab_orb=sab_nl, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
1235 scf_env=scf_env, scf_control=scf_control)
1236 cpassert(ASSOCIATED(para_env))
1237 cpassert(ASSOCIATED(blacs_env))
1238 cpassert(ASSOCIATED(cell))
1239 cpassert(ASSOCIATED(particle_set))
1240 cpassert(ASSOCIATED(mos))
1241 cpassert(ASSOCIATED(dft_control))
1242 cpassert(ASSOCIATED(sab_nl))
1243 cpassert(ASSOCIATED(matrix_ks))
1244 cpassert(ASSOCIATED(matrix_s))
1245 cpassert(ASSOCIATED(scf_env))
1246 cpassert(ASSOCIATED(scf_control))
1247
1248 NULLIFY (kpoints_out)
1249 CALL kpoint_create(kpoints_out)
1250 kpoints_out%kp_scheme = kp_scheme
1251 kpoints_out%symmetry = .false.
1252 kpoints_out%full_grid = .true.
1253 kpoints_out%verbose = .false.
1254 kpoints_out%use_real_wfn = .false.
1255 kpoints_out%eps_geo = eps_geo
1256 kpoints_out%parallel_group_size = para_env%num_pe
1257
1258 SELECT CASE (trim(kp_scheme))
1259 CASE ("MONKHORST-PACK", "MACDONALD")
1260 kpoints_out%nkp_grid(1:3) = nkp_grid(1:3)
1261 kpoints_out%kp_shift(1:3) = kp_shift(1:3)
1262 kpoints_out%gamma_centered = gamma_centered
1263 CALL kpoint_initialize(kpoints_out, particle_set, cell)
1264 CASE ("GENERAL")
1265 IF (.NOT. ASSOCIATED(kpoints_scf%xkp_input) .OR. &
1266 .NOT. ASSOCIATED(kpoints_scf%wkp_input)) THEN
1267 cpabort("TREXIO%FULL_KPOINT_GRID cannot recover the unreduced GENERAL k-point set.")
1268 END IF
1269 xkp_source => kpoints_scf%xkp_input
1270 wkp_source => kpoints_scf%wkp_input
1271 nfull = SIZE(wkp_source)
1272 wsum = sum(wkp_source)
1273 IF (wsum <= 0.0_dp) cpabort("TREXIO%FULL_KPOINT_GRID found invalid GENERAL k-point weights.")
1274 kpoints_out%nkp = nfull
1275 ALLOCATE (kpoints_out%xkp(3, nfull), kpoints_out%wkp(nfull))
1276 kpoints_out%xkp(1:3, 1:nfull) = xkp_source(1:3, 1:nfull)
1277 kpoints_out%wkp(1:nfull) = wkp_source(1:nfull)/wsum
1278 END SELECT
1279
1280 CALL kpoint_env_initialize(kpoints_out, para_env, blacs_env)
1281 CALL kpoint_initialize_mos(kpoints_out, mos)
1282 CALL kpoint_initialize_mo_set(kpoints_out)
1283 CALL kpoint_init_cell_index(kpoints_out, sab_nl, para_env, dft_control%nimages)
1284
1285 reused_scf_mos = .false.
1286 reuse_reason = ""
1287 aligned_blocks = 0
1288 aligned_max_size = 0
1289 aligned_min_svalue = 0.0_dp
1290 diis_step = .false.
1291 IF (reuse_scf_mos) THEN
1292 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints_scf, scf_env, scf_control, .false., &
1293 diis_step)
1294 CALL get_kpoint_info(kpoints_out, cell_to_index=cell_to_index)
1295 CALL prepare_wannier90_scf_mos(kpoints_out, kpoints_scf, matrix_s, matrix_ks, &
1296 cell_to_index, sab_nl, para_env, reused_scf_mos, &
1297 reuse_reason, aligned_blocks, aligned_max_size, &
1298 aligned_min_svalue)
1299 END IF
1300 IF (reused_scf_mos) THEN
1301 IF (output_unit > 0) THEN
1302 WRITE (output_unit, '(T2,A)') &
1303 "TREXIO| Reused SCF MO coefficients for the full k-point grid."
1304 IF (aligned_blocks > 0) THEN
1305 WRITE (output_unit, '(T2,A,I0,A,I0,A,ES10.3)') &
1306 "TREXIO| Ritz-stabilized ", aligned_blocks, &
1307 " degenerate SCF MO subspace(s); largest block has ", aligned_max_size, &
1308 " band(s), min metric eigenvalue ", aligned_min_svalue
1309 END IF
1310 END IF
1311 ELSE
1312 IF (output_unit > 0) THEN
1313 IF (reuse_scf_mos) THEN
1314 WRITE (output_unit, '(T2,A,A)') &
1315 "TREXIO| Could not reuse SCF MOs: ", trim(reuse_reason)
1316 END IF
1317 WRITE (output_unit, '(T2,A)') &
1318 "TREXIO| Diagonalizing the full k-point grid for export."
1319 END IF
1320 diis_step = .false.
1321 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoints_out, scf_env, scf_control, .false., &
1322 diis_step)
1323 END IF
1324 created = .true.
1325
1326 CALL timestop(handle)
1327#else
1328 mark_used(qs_env)
1329 mark_used(trexio_section)
1330 mark_used(do_kpoints)
1331 mark_used(kpoints_scf)
1332 NULLIFY (kpoints_out)
1333 created = .false.
1334#endif
1335
1336 END SUBROUTINE prepare_trexio_kpoint_grid
1337
1338! **************************************************************************************************
1339!> \brief Read a trexio file
1340!> \param qs_env the qs environment with all the info of the computation
1341!> \param trexio_filename the trexio filename without the extension
1342!> \param mo_set_trexio the MO set to read from the trexio file
1343!> \param energy_derivative the energy derivative to read from the trexio file
1344! **************************************************************************************************
1345 SUBROUTINE read_trexio(qs_env, trexio_filename, mo_set_trexio, energy_derivative)
1346 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1347 CHARACTER(len=*), INTENT(IN), OPTIONAL :: trexio_filename
1348 TYPE(mo_set_type), INTENT(OUT), DIMENSION(:), POINTER, OPTIONAL :: mo_set_trexio
1349 TYPE(dbcsr_p_type), INTENT(OUT), DIMENSION(:), POINTER, OPTIONAL :: energy_derivative
1350
1351#ifdef __TREXIO
1352
1353 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_trexio'
1354
1355 INTEGER :: handle, output_unit, unit_de
1356 CHARACTER(len=default_path_length) :: filename, filename_de
1357 INTEGER(trexio_t) :: f ! The TREXIO file handle
1358 INTEGER(trexio_exit_code) :: rc ! TREXIO return code
1359
1360 LOGICAL :: ionode
1361
1362 CHARACTER(LEN=2) :: element_symbol
1363 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: label
1364
1365 INTEGER :: ao_num, mo_num, nmo, nspins, ispin, nsgf, &
1366 save_cartesian, i, j, k, l, m, imo, ishell, &
1367 nshell, shell_num, nucleus_num, natoms, ikind, &
1368 iatom, nelectron, nrows, ncols, &
1369 row, col, row_size, col_size, &
1370 row_offset, col_offset, myprint
1371 INTEGER, DIMENSION(2) :: nmo_spin, electron_num
1372 INTEGER, DIMENSION(:), ALLOCATABLE :: mo_spin, shell_ang_mom, trexio_to_cp2k_ang_mom
1373
1374 REAL(kind=dp) :: zeff, maxocc
1375 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: mo_energy, mo_occupation, charge
1376 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: mo_coefficient, mos_sgf, coord, dedp, temp
1377 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
1378
1379 TYPE(cp_logger_type), POINTER :: logger
1380 TYPE(cp_fm_type), POINTER :: mo_coeff_ref, mo_coeff_target
1381 TYPE(mp_para_env_type), POINTER :: para_env
1382 TYPE(dft_control_type), POINTER :: dft_control
1383 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1384 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1385 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1386 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1387 TYPE(dbcsr_iterator_type) :: iter
1388
1389 CALL timeset(routinen, handle)
1390
1391 NULLIFY (logger, mo_coeff_ref, mo_coeff_target, para_env, dft_control, matrix_s, kind_set, mos, particle_set)
1392
1393 logger => cp_get_default_logger()
1394 output_unit = cp_logger_get_default_io_unit(logger)
1395 myprint = logger%iter_info%print_level
1396
1397 cpassert(ASSOCIATED(qs_env))
1398
1399 ! get filename
1400 IF (.NOT. PRESENT(trexio_filename)) THEN
1401 filename = trim(logger%iter_info%project_name)//'-TREXIO.h5'
1402 filename_de = trim(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
1403 ELSE
1404 filename = trim(trexio_filename)//'.h5'
1405 filename_de = trim(trexio_filename)//'.dEdP.dat'
1406 END IF
1407
1408 CALL get_qs_env(qs_env, para_env=para_env)
1409 ionode = para_env%is_source()
1410
1411 ! Open the TREXIO file and check that we have the same molecule as in qs_env
1412 IF (ionode) THEN
1413 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Opening file named ', trim(filename)
1414 f = trexio_open(filename, 'r', trexio_hdf5, rc)
1415 CALL trexio_error(rc)
1416
1417 IF (myprint > medium_print_level) THEN
1418 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading molecule information...'
1419 END IF
1420 rc = trexio_read_nucleus_num(f, nucleus_num)
1421 CALL trexio_error(rc)
1422
1423 IF (myprint > medium_print_level) THEN
1424 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear coordinates...'
1425 END IF
1426 ALLOCATE (coord(3, nucleus_num))
1427 rc = trexio_read_nucleus_coord(f, coord)
1428 CALL trexio_error(rc)
1429
1430 IF (myprint > medium_print_level) THEN
1431 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear labels...'
1432 END IF
1433 ALLOCATE (label(nucleus_num))
1434 rc = trexio_read_nucleus_label(f, label, 3)
1435 CALL trexio_error(rc)
1436
1437 IF (myprint > medium_print_level) THEN
1438 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear charges...'
1439 END IF
1440 ALLOCATE (charge(nucleus_num))
1441 rc = trexio_read_nucleus_charge(f, charge)
1442 CALL trexio_error(rc)
1443
1444 ! get the same info from qs_env
1445 CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
1446
1447 ! check that we have the same number of atoms
1448 cpassert(nucleus_num == natoms)
1449
1450 DO iatom = 1, natoms
1451 ! compare the coordinates within a certain tolerance
1452 DO i = 1, 3
1453 cpassert(abs(coord(i, iatom) - particle_set(iatom)%r(i)) < 1.0e-6_dp)
1454 END DO
1455
1456 ! figure out the element symbol and to which kind_set entry this atomic_kind corresponds to
1457 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
1458 ! check that the element symbol is the same
1459 cpassert(trim(element_symbol) == trim(label(iatom)))
1460
1461 ! get the effective nuclear charge for this kind
1462 CALL get_qs_kind(kind_set(ikind), zeff=zeff)
1463 ! check that the nuclear charge is also the same
1464 cpassert(charge(iatom) == zeff)
1465 END DO
1466
1467 WRITE (output_unit, "((T2,A))") 'TREXIO| Molecule is the same as in qs_env'
1468 ! if we get here, we have the same molecule
1469 DEALLOCATE (coord)
1470 DEALLOCATE (label)
1471 DEALLOCATE (charge)
1472
1473 ! get info from trexio to map cp2k and trexio AOs
1474 rc = trexio_read_ao_cartesian(f, save_cartesian)
1475 CALL trexio_error(rc)
1476
1477 rc = trexio_read_ao_num(f, ao_num)
1478 CALL trexio_error(rc)
1479
1480 rc = trexio_read_basis_shell_num(f, shell_num)
1481 CALL trexio_error(rc)
1482 END IF
1483
1484 CALL para_env%bcast(save_cartesian, para_env%source)
1485 CALL para_env%bcast(ao_num, para_env%source)
1486 CALL para_env%bcast(shell_num, para_env%source)
1487
1488 IF (save_cartesian == 1) THEN
1489 cpabort('Reading Cartesian AOs is not yet supported.')
1490 END IF
1491
1492 ! check that the number of AOs and shells is the same
1493 CALL get_qs_env(qs_env, qs_kind_set=kind_set)
1494 CALL get_qs_kind_set(kind_set, nsgf=nsgf, nshell=nshell)
1495 cpassert(ao_num == nsgf)
1496 cpassert(shell_num == nshell)
1497
1498 ALLOCATE (shell_ang_mom(shell_num))
1499 shell_ang_mom(:) = 0
1500
1501 IF (ionode) THEN
1502 IF (myprint > medium_print_level) THEN
1503 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading shell angular momenta...'
1504 END IF
1505 rc = trexio_read_basis_shell_ang_mom(f, shell_ang_mom)
1506 CALL trexio_error(rc)
1507 END IF
1508
1509 CALL para_env%bcast(shell_ang_mom, para_env%source)
1510
1511 ! AO order map from TREXIO to CP2K convention
1512 ! from m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
1513 ! to m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
1514 ALLOCATE (trexio_to_cp2k_ang_mom(nsgf))
1515 i = 0
1516 DO ishell = 1, shell_num
1517 l = shell_ang_mom(ishell)
1518 DO k = 1, 2*l + 1
1519 m = (-1)**k*floor(real(k, kind=dp)/2.0_dp)
1520 trexio_to_cp2k_ang_mom(i + l + 1 + m) = i + k
1521 END DO
1522 i = i + 2*l + 1
1523 END DO
1524 cpassert(i == nsgf)
1525
1526 ! check whether we want to read MOs
1527 IF (PRESENT(mo_set_trexio)) THEN
1528 IF (output_unit > 1) THEN
1529 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading molecular orbitals...'
1530 END IF
1531
1532 ! at the moment, we assume that the basis set is the same
1533 ! first we read all arrays lengths we need from the trexio file
1534 IF (ionode) THEN
1535 rc = trexio_read_mo_num(f, mo_num)
1536 CALL trexio_error(rc)
1537
1538 rc = trexio_read_electron_up_num(f, electron_num(1))
1539 CALL trexio_error(rc)
1540
1541 rc = trexio_read_electron_dn_num(f, electron_num(2))
1542 CALL trexio_error(rc)
1543 END IF
1544
1545 ! broadcast information to all processors and allocate arrays
1546 CALL para_env%bcast(mo_num, para_env%source)
1547 CALL para_env%bcast(electron_num, para_env%source)
1548
1549 ! check that the number of MOs is the same
1550 CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1551 nspins = dft_control%nspins
1552 nmo_spin(:) = 0
1553 DO ispin = 1, nspins
1554 CALL get_mo_set(mos(ispin), nmo=nmo)
1555 nmo_spin(ispin) = nmo
1556 END DO
1557 cpassert(mo_num == sum(nmo_spin))
1558
1559 ALLOCATE (mo_coefficient(ao_num, mo_num))
1560 ALLOCATE (mo_energy(mo_num))
1561 ALLOCATE (mo_occupation(mo_num))
1562 ALLOCATE (mo_spin(mo_num))
1563
1564 mo_coefficient(:, :) = 0.0_dp
1565 mo_energy(:) = 0.0_dp
1566 mo_occupation(:) = 0.0_dp
1567 mo_spin(:) = 0
1568
1569 ! read the MOs info
1570 IF (ionode) THEN
1571 IF (myprint > medium_print_level) THEN
1572 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO coefficients...'
1573 END IF
1574 rc = trexio_read_mo_coefficient(f, mo_coefficient)
1575 CALL trexio_error(rc)
1576
1577 IF (myprint > medium_print_level) THEN
1578 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO energies...'
1579 END IF
1580 rc = trexio_read_mo_energy(f, mo_energy)
1581 CALL trexio_error(rc)
1582
1583 IF (myprint > medium_print_level) THEN
1584 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO occupations...'
1585 END IF
1586 rc = trexio_read_mo_occupation(f, mo_occupation)
1587 CALL trexio_error(rc)
1588
1589 IF (myprint > medium_print_level) THEN
1590 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO spins...'
1591 END IF
1592 rc = trexio_read_mo_spin(f, mo_spin)
1593 CALL trexio_error(rc)
1594 END IF
1595
1596 ! broadcast the data to all processors
1597 CALL para_env%bcast(mo_coefficient, para_env%source)
1598 CALL para_env%bcast(mo_energy, para_env%source)
1599 CALL para_env%bcast(mo_occupation, para_env%source)
1600 CALL para_env%bcast(mo_spin, para_env%source)
1601
1602 ! assume nspins and nmo_spin match the ones in the trexio file
1603 ! reorder magnetic quantum number
1604 DO ispin = 1, nspins
1605 ! global MOs index
1606 imo = (ispin - 1)*nmo_spin(1)
1607 ! get number of MOs for this spin
1608 nmo = nmo_spin(ispin)
1609 ! allocate local temp array to read MOs
1610 ALLOCATE (mos_sgf(nsgf, nmo))
1611 mos_sgf(:, :) = 0.0_dp
1612
1613 ! we need to reorder the MOs according to CP2K convention
1614 DO i = 1, nsgf
1615 mos_sgf(i, :) = mo_coefficient(trexio_to_cp2k_ang_mom(i), imo + 1:imo + nmo)
1616 END DO
1617
1618 IF (nspins == 1) THEN
1619 maxocc = 2.0_dp
1620 nelectron = electron_num(1) + electron_num(2)
1621 ELSE
1622 maxocc = 1.0_dp
1623 nelectron = electron_num(ispin)
1624 END IF
1625 ! the right number of active electrons per spin channel is initialized further down
1626 CALL allocate_mo_set(mo_set_trexio(ispin), nsgf, nmo, nelectron, 0.0_dp, maxocc, 0.0_dp)
1627
1628 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff_ref)
1629 CALL init_mo_set(mo_set_trexio(ispin), fm_ref=mo_coeff_ref, name="TREXIO MOs")
1630
1631 CALL get_mo_set(mo_set_trexio(ispin), mo_coeff=mo_coeff_target)
1632 DO j = 1, nmo
1633 ! make sure I copy the right spin channel
1634 cpassert(mo_spin(j) == ispin - 1)
1635 mo_set_trexio(ispin)%eigenvalues(j) = mo_energy(imo + j)
1636 mo_set_trexio(ispin)%occupation_numbers(j) = mo_occupation(imo + j)
1637 DO i = 1, nsgf
1638 CALL cp_fm_set_element(mo_coeff_target, i, j, mos_sgf(i, j))
1639 END DO
1640 END DO
1641
1642 DEALLOCATE (mos_sgf)
1643 END DO
1644
1645 DEALLOCATE (mo_coefficient)
1646 DEALLOCATE (mo_energy)
1647 DEALLOCATE (mo_occupation)
1648 DEALLOCATE (mo_spin)
1649
1650 END IF ! if MOs should be read
1651
1652 ! check whether we want to read derivatives
1653 IF (PRESENT(energy_derivative)) THEN
1654 IF (output_unit > 1) THEN
1655 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading energy derivatives...'
1656 END IF
1657
1658 ! Temporary solution: allocate here the energy derivatives matrix here
1659 ! assuming that nsgf is the same as the number read from the dEdP file
1660 ! TODO: once available in TREXIO, first read size and then allocate
1661 ! in the same way done for the MOs
1662 ALLOCATE (temp(nsgf, nsgf))
1663 temp(:, :) = 0.0_dp
1664
1665 ! check if file exists and open it
1666 IF (ionode) THEN
1667 IF (file_exists(filename_de)) THEN
1668 CALL open_file(file_name=filename_de, file_status="OLD", unit_number=unit_de)
1669 ELSE
1670 cpabort("Energy derivatives file "//trim(filename_de)//" not found")
1671 END IF
1672
1673 ! read the header and check everything is fine
1674 IF (myprint > medium_print_level) THEN
1675 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading header information...'
1676 END IF
1677 READ (unit_de, *) nrows, ncols
1678 IF (myprint > medium_print_level) THEN
1679 WRITE (output_unit, "((T2,A))") 'TREXIO| Check size of dEdP matrix...'
1680 END IF
1681 cpassert(nrows == nsgf)
1682 cpassert(ncols == nsgf)
1683
1684 ! read the data
1685 IF (myprint > medium_print_level) THEN
1686 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading dEdP matrix...'
1687 END IF
1688 ! Read the data matrix
1689 DO i = 1, nrows
1690 READ (unit_de, *) (temp(i, j), j=1, ncols)
1691 END DO
1692
1693 CALL close_file(unit_number=unit_de)
1694 END IF
1695
1696 ! send data to all processes
1697 CALL para_env%bcast(temp, para_env%source)
1698
1699 ! Reshuffle
1700 ALLOCATE (dedp(nsgf, nsgf))
1701 dedp(:, :) = 0.0_dp
1702
1703 ! Reorder rows and columns according to trexio_to_cp2k_ang_mom mapping
1704 DO j = 1, nsgf
1705 DO i = 1, nsgf
1706 ! either this
1707 dedp(i, j) = temp(trexio_to_cp2k_ang_mom(i), trexio_to_cp2k_ang_mom(j))
1708 ! or this
1709 ! dEdP(cp2k_to_trexio_ang_mom(i), cp2k_to_trexio_ang_mom(j)) = temp(i, j)
1710 END DO
1711 END DO
1712
1713 DEALLOCATE (temp)
1714
1715 CALL get_qs_env(qs_env, matrix_s=matrix_s)
1716 DO ispin = 1, nspins
1717 ALLOCATE (energy_derivative(ispin)%matrix)
1718
1719 ! we use the overlap matrix as a template, copying it but removing the sparsity
1720 CALL dbcsr_copy(energy_derivative(ispin)%matrix, matrix_s(1)%matrix, &
1721 name='Energy Derivative', keep_sparsity=.false.)
1722 CALL dbcsr_set(energy_derivative(ispin)%matrix, 0.0_dp)
1723
1724 CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
1725 DO WHILE (dbcsr_iterator_blocks_left(iter))
1726 CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
1727 row_size=row_size, col_size=col_size, &
1728 row_offset=row_offset, col_offset=col_offset)
1729
1730 ! Copy data from array to block
1731 DO i = 1, row_size
1732 DO j = 1, col_size
1733 data_block(i, j) = dedp(row_offset + i - 1, col_offset + j - 1)
1734 END DO
1735 END DO
1736 END DO
1737 CALL dbcsr_iterator_stop(iter)
1738 END DO
1739
1740 DEALLOCATE (dedp)
1741 END IF ! finished reading energy derivatives
1742
1743 ! Clean up
1744 IF (ALLOCATED(shell_ang_mom)) DEALLOCATE (shell_ang_mom)
1745 IF (ALLOCATED(trexio_to_cp2k_ang_mom)) DEALLOCATE (trexio_to_cp2k_ang_mom)
1746
1747 ! Close the TREXIO file
1748 IF (ionode) THEN
1749 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Closing file named ', trim(filename)
1750 rc = trexio_close(f)
1751 CALL trexio_error(rc)
1752 END IF
1753
1754 CALL timestop(handle)
1755
1756#else
1757 mark_used(qs_env)
1758 mark_used(trexio_filename)
1759 mark_used(mo_set_trexio)
1760 mark_used(energy_derivative)
1761 cpwarn('TREXIO support has not been enabled in this build.')
1762 cpabort('TREXIO Not Available')
1763#endif
1764
1765 END SUBROUTINE read_trexio
1766
1767#ifdef __TREXIO
1768! **************************************************************************************************
1769!> \brief Handles TREXIO errors
1770!> \param rc the TREXIO return code
1771! **************************************************************************************************
1772 SUBROUTINE trexio_error(rc)
1773 INTEGER(trexio_exit_code), INTENT(IN) :: rc
1774
1775 CHARACTER(LEN=128) :: err_msg
1776
1777 IF (rc /= trexio_success) THEN
1778 CALL trexio_string_of_error(rc, err_msg)
1779 cpabort('TREXIO Error: '//trim(err_msg))
1780 END IF
1781
1782 END SUBROUTINE trexio_error
1783
1784! **************************************************************************************************
1785!> \brief Computes the nuclear repulsion energy of a molecular system
1786!> \param particle_set the set of particles in the system
1787!> \param kind_set the set of qs_kinds in the system
1788!> \param e_nn the nuclear repulsion energy
1789! **************************************************************************************************
1790 SUBROUTINE nuclear_repulsion_energy(particle_set, kind_set, e_nn)
1791 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1792 POINTER :: particle_set
1793 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1794 POINTER :: kind_set
1795 REAL(kind=dp), INTENT(OUT) :: e_nn
1796
1797 INTEGER :: i, ikind, j, jkind, natoms
1798 REAL(kind=dp) :: r_ij, zeff_i, zeff_j
1799
1800 natoms = SIZE(particle_set)
1801 e_nn = 0.0_dp
1802 DO i = 1, natoms
1803 CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
1804 CALL get_qs_kind(kind_set(ikind), zeff=zeff_i)
1805 DO j = i + 1, natoms
1806 r_ij = norm2(particle_set(i)%r - particle_set(j)%r)
1807
1808 CALL get_atomic_kind(particle_set(j)%atomic_kind, kind_number=jkind)
1809 CALL get_qs_kind(kind_set(jkind), zeff=zeff_j)
1810
1811 e_nn = e_nn + zeff_i*zeff_j/r_ij
1812 END DO
1813 END DO
1814
1815 END SUBROUTINE nuclear_repulsion_energy
1816
1817! **************************************************************************************************
1818!> \brief Returns the normalization coefficient for a spherical GTO
1819!> \param l the angular momentum quantum number
1820!> \param expnt the exponent of the Gaussian function
1821!> \return ...
1822! **************************************************************************************************
1823 FUNCTION sgf_norm(l, expnt) RESULT(norm)
1824 INTEGER, INTENT(IN) :: l
1825 REAL(kind=dp), INTENT(IN) :: expnt
1826 REAL(kind=dp) :: norm
1827
1828 IF (l >= 0) THEN
1829 norm = sqrt(2**(2*l + 3)*fac(l + 1)*(2*expnt)**(l + 1.5)/(fac(2*l + 2)*sqrt(pi)))
1830 ELSE
1831 cpabort("The angular momentum should be >= 0!")
1832 END IF
1833
1834 END FUNCTION sgf_norm
1835
1836! **************************************************************************************************
1837!> \brief Computes a spherical to cartesian MO transformation (solid harmonics in reality)
1838!> \param mos_sgf the MO coefficients in spherical AO basis
1839!> \param particle_set the set of particles in the system
1840!> \param qs_kind_set the set of qs_kinds in the system
1841!> \param mos_cgf the transformed MO coefficients in Cartesian AO basis
1842! **************************************************************************************************
1843 SUBROUTINE spherical_to_cartesian_mo(mos_sgf, particle_set, qs_kind_set, mos_cgf)
1844 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: mos_sgf
1845 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1846 POINTER :: particle_set
1847 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1848 POINTER :: qs_kind_set
1849 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: mos_cgf
1850
1851 INTEGER :: iatom, icgf, ikind, iset, isgf, ishell, &
1852 lshell, ncgf, nmo, nset, nsgf
1853 INTEGER, DIMENSION(:), POINTER :: nshell
1854 INTEGER, DIMENSION(:, :), POINTER :: l
1855 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1856
1857 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1858
1859 mos_cgf(:, :) = 0.0_dp
1860 nmo = SIZE(mos_sgf, 2)
1861
1862 ! Transform spherical MOs to Cartesian MOs
1863 icgf = 1
1864 isgf = 1
1865 DO iatom = 1, SIZE(particle_set)
1866 NULLIFY (orb_basis_set)
1867 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1868 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1869
1870 IF (ASSOCIATED(orb_basis_set)) THEN
1871 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1872 nset=nset, &
1873 nshell=nshell, &
1874 l=l)
1875 DO iset = 1, nset
1876 DO ishell = 1, nshell(iset)
1877 lshell = l(ishell, iset)
1878 CALL dgemm("T", "N", nco(lshell), nmo, nso(lshell), 1.0_dp, &
1879 orbtramat(lshell)%c2s, nso(lshell), &
1880 mos_sgf(isgf, 1), nsgf, 0.0_dp, &
1881 mos_cgf(icgf, 1), ncgf)
1882 icgf = icgf + nco(lshell)
1883 isgf = isgf + nso(lshell)
1884 END DO
1885 END DO
1886 ELSE
1887 ! assume atom without basis set
1888 cpabort("Unknown basis set type")
1889 END IF
1890 END DO ! iatom
1891
1892 END SUBROUTINE spherical_to_cartesian_mo
1893
1894! **************************************************************************************************
1895!> \brief Computes a cartesian to spherical MO transformation
1896!> \param mos_cgf the transformed MO coefficients in Cartesian AO basis
1897!> \param particle_set the set of particles in the system
1898!> \param qs_kind_set the set of qs_kinds in the system
1899!> \param mos_sgf the MO coefficients in spherical AO basis
1900! **************************************************************************************************
1901 SUBROUTINE cartesian_to_spherical_mo(mos_cgf, particle_set, qs_kind_set, mos_sgf)
1902 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: mos_cgf
1903 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1904 POINTER :: particle_set
1905 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1906 POINTER :: qs_kind_set
1907 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: mos_sgf
1908
1909 INTEGER :: iatom, icgf, ikind, iset, isgf, ishell, &
1910 lshell, ncgf, nmo, nset, nsgf
1911 INTEGER, DIMENSION(:), POINTER :: nshell
1912 INTEGER, DIMENSION(:, :), POINTER :: l
1913 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1914
1915 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1916
1917 mos_sgf(:, :) = 0.0_dp
1918 nmo = SIZE(mos_cgf, 2)
1919
1920 ! Transform Cartesian MOs to spherical MOs
1921 icgf = 1
1922 isgf = 1
1923 DO iatom = 1, SIZE(particle_set)
1924 NULLIFY (orb_basis_set)
1925 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1926 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1927
1928 IF (ASSOCIATED(orb_basis_set)) THEN
1929 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1930 nset=nset, &
1931 nshell=nshell, &
1932 l=l)
1933 DO iset = 1, nset
1934 DO ishell = 1, nshell(iset)
1935 lshell = l(ishell, iset)
1936 CALL dgemm("N", "N", nso(lshell), nmo, nco(lshell), 1.0_dp, &
1937 orbtramat(lshell)%s2c, nso(lshell), &
1938 mos_cgf(icgf, 1), ncgf, 0.0_dp, &
1939 mos_sgf(isgf, 1), nsgf)
1940 icgf = icgf + nco(lshell)
1941 isgf = isgf + nso(lshell)
1942 END DO
1943 END DO
1944 ELSE
1945 ! assume atom without basis set
1946 cpabort("Unknown basis set type")
1947 END IF
1948 END DO ! iatom
1949
1950 END SUBROUTINE cartesian_to_spherical_mo
1951#endif
1952
1953END MODULE trexio_utils
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cell_types::pbc
Definition cell_types.F:105
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
external_potential_types::get_potential
Definition external_potential_types.F:279
ai_onecenter
Definition ai_onecenter.F:28
ai_onecenter::sg_overlap
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition ai_onecenter.F:65
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:131
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:631
cell_types::pbc_stable
real(kind=dp) function, dimension(3), public pbc_stable(r, cell)
Apply a stable periodic-image convention for k-point Bloch gauges.
Definition cell_types.F:458
cp2k_info
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:22
cp2k_info::cp2k_version
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:50
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:793
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:300
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if a file exists, considering also the file discovery mechanism.
Definition cp_files.F:518
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1046
cp_fm_types::cp_fm_to_fm_submat_general
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
Definition cp_fm_types.F:2049
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:536
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:933
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:692
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:77
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_initialize_mo_set
subroutine, public kpoint_initialize_mo_set(kpoint)
...
Definition kpoint_methods.F:988
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:1023
kpoint_methods::kpoint_initialize_mos
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group).
Definition kpoint_methods.F:832
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:377
kpoint_methods::kpoint_env_initialize
subroutine, public kpoint_env_initialize(kpoint, para_env, blacs_env, with_aux_fit)
Initialize the kpoint environment.
Definition kpoint_methods.F:668
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_env
subroutine, public get_kpoint_env(kpoint_env, nkpoint, wkp, xkp, is_local, mos)
Get information from a single kpoint environment.
Definition kpoint_types.F:1202
kpoint_types::kpoint_release
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Definition kpoint_types.F:321
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:613
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:268
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::symmetrize_matrix
subroutine, public symmetrize_matrix(a, option)
Symmetrize the matrix a.
Definition mathlib.F:1206
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:483
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:976
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter, complex_coeff)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:261
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:219
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:440
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_scf_diagonalization
Different diagonalization schemes that can be used for the iterative solution of the eigenvalue probl...
Definition qs_scf_diagonalization.F:15
qs_scf_diagonalization::do_general_diag_kp
subroutine, public do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, diis_step, diis_error, qs_env, probe, added_mos_auto_grow, potential)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
Definition qs_scf_diagonalization.F:707
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_wannier90
Interface to Wannier90 code.
Definition qs_wannier90.F:14
qs_wannier90::prepare_wannier90_scf_mos
subroutine, public prepare_wannier90_scf_mos(kpoint, qs_kpoint, matrix_s, matrix_ks, cell_to_index, sab_nl, para_env, success, reason, aligned_degenerate_blocks, aligned_degenerate_max_size, aligned_degenerate_min_svalue, spin_channel)
Reconstruct a full Wannier90 k-point MO set from the SCF k-point MOs.
Definition qs_wannier90.F:1650
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
trexio_utils
The module to read/write TREX IO files for interfacing CP2K with other programs.
Definition trexio_utils.F:16
trexio_utils::write_trexio
subroutine, public write_trexio(qs_env, trexio_section, energy_derivative)
Write a trexio file.
Definition trexio_utils.F:137
trexio_utils::read_trexio
subroutine, public read_trexio(qs_env, trexio_filename, mo_set_trexio, energy_derivative)
Read a trexio file.
Definition trexio_utils.F:1346
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_p_type
Definition kpoint_types.F:119
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:120
scf_control_types::scf_control_type
Definition scf_control_types.F:135