(git:ed6f26b)
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-2025 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!> \author Stefano Battaglia
13! **************************************************************************************************
14MODULE trexio_utils
15
16 USE atomic_kind_types, ONLY: get_atomic_kind
17 USE basis_set_types, ONLY: gto_basis_set_type, get_gto_basis_set
18 USE cell_types, ONLY: cell_type
19 USE cp2k_info, ONLY: cp2k_version
20 USE cp_blacs_env, ONLY: cp_blacs_env_type
21 USE cp_control_types, ONLY: dft_control_type
22 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
23 USE cp_files, ONLY: close_file, file_exists, open_file
24 USE cp_fm_types, ONLY: cp_fm_get_info, cp_fm_type, cp_fm_create, cp_fm_set_all, &
25 cp_fm_get_submatrix, cp_fm_to_fm_submat_general, cp_fm_release, &
26 cp_fm_set_element
27 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
28 cp_fm_struct_release, &
29 cp_fm_struct_type
30 USE cp_log_handling, ONLY: cp_get_default_logger, &
31 cp_logger_get_default_io_unit, &
32 cp_logger_type
33 USE cp_dbcsr_api, ONLY: dbcsr_p_type, dbcsr_iterator_type, dbcsr_iterator_start, &
34 dbcsr_iterator_stop, dbcsr_iterator_blocks_left, &
35 dbcsr_iterator_next_block, dbcsr_copy, dbcsr_set, &
36 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
37 dbcsr_type_symmetric, dbcsr_get_matrix_type
38 USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
39 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
40 USE cp_output_handling, ONLY: medium_print_level
41 USE external_potential_types, ONLY: sgp_potential_type, get_potential
42 USE input_section_types, ONLY: section_vals_type, section_vals_get, &
43 section_vals_val_get
44 USE kinds, ONLY: default_path_length, dp
45 USE kpoint_types, ONLY: kpoint_type, kpoint_env_p_type, &
46 get_kpoint_info, get_kpoint_env
47 USE mathconstants, ONLY: fourpi, pi
48 USE mathlib, ONLY: symmetrize_matrix
49 USE message_passing, ONLY: mp_para_env_type
50 USE orbital_pointers, ONLY: nco, nso
51 USE orbital_transformation_matrices, ONLY: orbtramat
52 USE particle_types, ONLY: particle_type
53 USE qs_energy_types, ONLY: qs_energy_type
54 USE qs_environment_types, ONLY: get_qs_env, &
55 qs_environment_type
56 USE qs_kind_types, ONLY: get_qs_kind, get_qs_kind_set, &
57 qs_kind_type
58 USE qs_mo_types, ONLY: mo_set_type, get_mo_set, init_mo_set, allocate_mo_set
59#ifdef __TREXIO
60 USE trexio, ONLY: trexio_open, trexio_close, &
61 trexio_hdf5, trexio_success, &
62 trexio_string_of_error, trexio_t, trexio_exit_code, &
63 trexio_write_metadata_code, trexio_write_metadata_code_num, &
64 trexio_write_nucleus_coord, trexio_read_nucleus_coord, &
65 trexio_write_nucleus_num, trexio_read_nucleus_num, &
66 trexio_write_nucleus_charge, trexio_read_nucleus_charge, &
67 trexio_write_nucleus_label, trexio_read_nucleus_label, &
68 trexio_write_nucleus_repulsion, &
69 trexio_write_cell_a, trexio_write_cell_b, trexio_write_cell_c, &
70 trexio_write_cell_g_a, trexio_write_cell_g_b, &
71 trexio_write_cell_g_c, trexio_write_cell_two_pi, &
72 trexio_write_pbc_periodic, trexio_write_pbc_k_point_num, &
73 trexio_write_pbc_k_point, trexio_write_pbc_k_point_weight, &
74 trexio_write_electron_num, trexio_read_electron_num, &
75 trexio_write_electron_up_num, trexio_read_electron_up_num, &
76 trexio_write_electron_dn_num, trexio_read_electron_dn_num, &
77 trexio_write_state_num, trexio_write_state_id, &
78 trexio_write_state_energy, &
79 trexio_write_basis_type, trexio_write_basis_prim_num, &
80 trexio_write_basis_shell_num, trexio_read_basis_shell_num, &
81 trexio_write_basis_nucleus_index, &
82 trexio_write_basis_shell_ang_mom, trexio_read_basis_shell_ang_mom, &
83 trexio_write_basis_shell_factor, &
84 trexio_write_basis_r_power, trexio_write_basis_shell_index, &
85 trexio_write_basis_exponent, trexio_write_basis_coefficient, &
86 trexio_write_basis_prim_factor, &
87 trexio_write_ecp_z_core, trexio_write_ecp_max_ang_mom_plus_1, &
88 trexio_write_ecp_num, trexio_write_ecp_ang_mom, &
89 trexio_write_ecp_nucleus_index, trexio_write_ecp_exponent, &
90 trexio_write_ecp_coefficient, trexio_write_ecp_power, &
91 trexio_write_ao_cartesian, trexio_write_ao_num, &
92 trexio_read_ao_cartesian, trexio_read_ao_num, &
93 trexio_write_ao_shell, trexio_write_ao_normalization, &
94 trexio_read_ao_shell, trexio_read_ao_normalization, &
95 trexio_write_mo_num, trexio_write_mo_energy, &
96 trexio_read_mo_num, trexio_read_mo_energy, &
97 trexio_write_mo_occupation, trexio_write_mo_spin, &
98 trexio_read_mo_occupation, trexio_read_mo_spin, &
99 trexio_write_mo_class, trexio_write_mo_coefficient, &
100 trexio_read_mo_class, trexio_read_mo_coefficient, &
101 trexio_write_mo_coefficient_im, trexio_write_mo_k_point, &
102 trexio_write_mo_type
103#endif
104#include "./base/base_uses.f90"
105
106 IMPLICIT NONE
107
108 PRIVATE
109
110 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'trexio_utils'
111
112 PUBLIC :: write_trexio, read_trexio
113
114CONTAINS
115
116! **************************************************************************************************
117!> \brief Write a trexio file
118!> \param qs_env the qs environment with all the info of the computation
119!> \param trexio_section the section with the trexio info
120!> \param energy_derivative ...
121! **************************************************************************************************
122 SUBROUTINE write_trexio(qs_env, trexio_section, energy_derivative)
123 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
124 TYPE(section_vals_type), INTENT(IN), POINTER :: trexio_section
125 TYPE(dbcsr_p_type), INTENT(IN), DIMENSION(:), POINTER, OPTIONAL :: energy_derivative
126
127 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_trexio'
128
129 INTEGER :: handle
130
131#ifdef __TREXIO
132 INTEGER :: output_unit, unit_trexio
133 CHARACTER(len=default_path_length) :: filename, filename_de
134 INTEGER(trexio_t) :: f ! The TREXIO file handle
135 INTEGER(trexio_exit_code) :: rc ! TREXIO return code
136 LOGICAL :: explicit, do_kpoints, ecp_semi_local, &
137 ecp_local, sgp_potential_present, ionode, &
138 use_real_wfn, save_cartesian
139 REAL(kind=dp) :: e_nn, zeff, expzet, prefac, zeta, gcca
140 TYPE(cell_type), POINTER :: cell => null()
141 TYPE(cp_logger_type), POINTER :: logger => null()
142 TYPE(dft_control_type), POINTER :: dft_control => null()
143 TYPE(gto_basis_set_type), POINTER :: basis_set
144 TYPE(kpoint_type), POINTER :: kpoints => null()
145 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => null()
146 TYPE(qs_energy_type), POINTER :: energy => null()
147 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set => null()
148 TYPE(sgp_potential_type), POINTER :: sgp_potential => null()
149 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos => null()
150 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp => null()
151 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env => null()
152 TYPE(mp_para_env_type), POINTER :: para_env => null(), para_env_inter_kp => null()
153 TYPE(cp_blacs_env_type), POINTER :: blacs_env => null()
154 TYPE(cp_fm_struct_type), POINTER :: fm_struct => null()
155 TYPE(cp_fm_type) :: fm_mo_coeff, fm_dummy, fm_mo_coeff_im
156 TYPE(dbcsr_iterator_type) :: iter
157
158 CHARACTER(LEN=2) :: element_symbol
159 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: label
160 INTEGER :: iatom, natoms, periodic, nkp, nel_tot, &
161 nspins, ikind, ishell_loc, ishell, &
162 shell_num, prim_num, nset, iset, ipgf, z, &
163 sl_lmax, ecp_num, nloc, nsemiloc, sl_l, iecp, &
164 igf, icgf, ncgf, ngf_shell, lshell, ao_num, nmo, &
165 mo_num, ispin, ikp, imo, ikp_loc, nsgf, &
166 i, j, k, l, m, unit_de, &
167 row, col, row_size, col_size, &
168 row_offset, col_offset
169 INTEGER, DIMENSION(2) :: nel_spin, kp_range, nmo_spin
170 INTEGER, DIMENSION(3) :: nkp_grid
171 INTEGER, DIMENSION(0:10) :: npot
172 INTEGER, DIMENSION(:), ALLOCATABLE :: nucleus_index, shell_ang_mom, r_power, &
173 shell_index, z_core, max_ang_mom_plus_1, &
174 ang_mom, powers, ao_shell, mo_spin, mo_kpoint, &
175 cp2k_to_trexio_ang_mom
176 INTEGER, DIMENSION(:), POINTER :: nshell => null(), npgf => null()
177 INTEGER, DIMENSION(:, :), POINTER :: l_shell_set => null()
178 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: charge, shell_factor, exponents, coefficients, &
179 prim_factor, ao_normalization, mo_energy, &
180 mo_occupation
181 REAL(kind=dp), DIMENSION(:), POINTER :: wkp => null(), norm_cgf => null()
182 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: coord, mo_coefficient, mo_coefficient_im, &
183 mos_sgf, diag_nsgf, diag_ncgf, temp, dedp
184 REAL(kind=dp), DIMENSION(:, :), POINTER :: zetas => null(), data_block => null()
185 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc => null()
186#endif
187
188 CALL timeset(routinen, handle)
189
190#ifdef __TREXIO
191 logger => cp_get_default_logger()
192 output_unit = cp_logger_get_default_io_unit(logger)
193
194 cpassert(ASSOCIATED(qs_env))
195
196 ! get filename
197 CALL section_vals_val_get(trexio_section, "FILENAME", c_val=filename, explicit=explicit)
198 IF (.NOT. explicit) THEN
199 filename = trim(logger%iter_info%project_name)//'-TREXIO.h5'
200 ELSE
201 filename = trim(filename)//'.h5'
202 END IF
203
204 CALL get_qs_env(qs_env, para_env=para_env)
205 ionode = para_env%is_source()
206
207 ! inquire whether a file with the same name already exists, if yes, delete it
208 IF (ionode) THEN
209 IF (file_exists(filename)) THEN
210 CALL open_file(filename, unit_number=unit_trexio)
211 CALL close_file(unit_number=unit_trexio, file_status="DELETE")
212 END IF
213
214 !========================================================================================!
215 ! Open the TREXIO file
216 !========================================================================================!
217 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing trexio file ', trim(filename)
218 f = trexio_open(filename, 'w', trexio_hdf5, rc)
219 CALL trexio_error(rc)
220
221 !========================================================================================!
222 ! Metadata group
223 !========================================================================================!
224 rc = trexio_write_metadata_code_num(f, 1)
225 CALL trexio_error(rc)
226
227 rc = trexio_write_metadata_code(f, cp2k_version, len_trim(cp2k_version) + 1)
228 CALL trexio_error(rc)
229
230 !========================================================================================!
231 ! Nucleus group
232 !========================================================================================!
233 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
234
235 rc = trexio_write_nucleus_num(f, natoms)
236 CALL trexio_error(rc)
237
238 ALLOCATE (coord(3, natoms))
239 ALLOCATE (label(natoms))
240 ALLOCATE (charge(natoms))
241 DO iatom = 1, natoms
242 ! store the coordinates
243 coord(:, iatom) = particle_set(iatom)%r(1:3)
244 ! figure out the element symbol and to which kind_set entry this atomic_kind corresponds to
245 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
246 ! store the element symbol
247 label(iatom) = element_symbol
248 ! get and store the effective nuclear charge of this kind_type (ikind)
249 CALL get_qs_kind(kind_set(ikind), zeff=zeff)
250 charge(iatom) = zeff
251 END DO
252
253 rc = trexio_write_nucleus_coord(f, coord)
254 CALL trexio_error(rc)
255 DEALLOCATE (coord)
256
257 rc = trexio_write_nucleus_charge(f, charge)
258 CALL trexio_error(rc)
259 DEALLOCATE (charge)
260
261 rc = trexio_write_nucleus_label(f, label, 3)
262 CALL trexio_error(rc)
263 DEALLOCATE (label)
264
265 ! nuclear repulsion energy well-defined for molecules only
266 IF (sum(cell%perd) == 0) THEN
267 CALL nuclear_repulsion_energy(particle_set, kind_set, e_nn)
268 rc = trexio_write_nucleus_repulsion(f, e_nn)
269 CALL trexio_error(rc)
270 END IF
271
272 !========================================================================================!
273 ! Cell group
274 !========================================================================================!
275 rc = trexio_write_cell_a(f, cell%hmat(:, 1))
276 CALL trexio_error(rc)
277
278 rc = trexio_write_cell_b(f, cell%hmat(:, 2))
279 CALL trexio_error(rc)
280
281 rc = trexio_write_cell_c(f, cell%hmat(:, 3))
282 CALL trexio_error(rc)
283
284 rc = trexio_write_cell_g_a(f, cell%h_inv(:, 1))
285 CALL trexio_error(rc)
286
287 rc = trexio_write_cell_g_b(f, cell%h_inv(:, 2))
288 CALL trexio_error(rc)
289
290 rc = trexio_write_cell_g_c(f, cell%h_inv(:, 3))
291 CALL trexio_error(rc)
292
293 rc = trexio_write_cell_two_pi(f, 0)
294 CALL trexio_error(rc)
295
296 !========================================================================================!
297 ! PBC group
298 !========================================================================================!
299 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints)
300
301 periodic = 0
302 IF (sum(cell%perd) /= 0) periodic = 1
303 rc = trexio_write_pbc_periodic(f, periodic)
304 CALL trexio_error(rc)
305
306 IF (do_kpoints) THEN
307 CALL get_kpoint_info(kpoints, nkp=nkp, nkp_grid=nkp_grid, wkp=wkp)
308
309 rc = trexio_write_pbc_k_point_num(f, nkp)
310 CALL trexio_error(rc)
311
312 rc = trexio_write_pbc_k_point(f, real(nkp_grid, kind=dp))
313 CALL trexio_error(rc)
314
315 rc = trexio_write_pbc_k_point_weight(f, wkp)
316 CALL trexio_error(rc)
317 END IF
318
319 !========================================================================================!
320 ! Electron group
321 !========================================================================================!
322 CALL get_qs_env(qs_env, dft_control=dft_control, nelectron_total=nel_tot)
323
324 rc = trexio_write_electron_num(f, nel_tot)
325 CALL trexio_error(rc)
326
327 nspins = dft_control%nspins
328 IF (nspins == 1) THEN
329 ! it is a spin-restricted calculation and we need to split the electrons manually,
330 ! because in CP2K they are all otherwise weirdly stored in nelectron_spin(1)
331 nel_spin(1) = nel_tot/2
332 nel_spin(2) = nel_tot/2
333 ELSE
334 ! for UKS/ROKS, the two spin channels are populated correctly and according to
335 ! the multiplicity
336 CALL get_qs_env(qs_env, nelectron_spin=nel_spin)
337 END IF
338 rc = trexio_write_electron_up_num(f, nel_spin(1))
339 CALL trexio_error(rc)
340 rc = trexio_write_electron_dn_num(f, nel_spin(2))
341 CALL trexio_error(rc)
342
343 !========================================================================================!
344 ! State group
345 !========================================================================================!
346 CALL get_qs_env(qs_env, energy=energy)
347
348 rc = trexio_write_state_num(f, 1)
349 CALL trexio_error(rc)
350
351 rc = trexio_write_state_id(f, 1)
352 CALL trexio_error(rc)
353
354 rc = trexio_write_state_energy(f, energy%total)
355 CALL trexio_error(rc)
356
357 END IF ! ionode
358
359 !========================================================================================!
360 ! Basis group
361 !========================================================================================!
362 CALL get_qs_env(qs_env, qs_kind_set=kind_set, natom=natoms, particle_set=particle_set)
363 CALL get_qs_kind_set(kind_set, nshell=shell_num, npgf_seg=prim_num)
364
365 IF (ionode) THEN
366 rc = trexio_write_basis_type(f, 'Gaussian', len_trim('Gaussian') + 1)
367 CALL trexio_error(rc)
368
369 rc = trexio_write_basis_shell_num(f, shell_num)
370 CALL trexio_error(rc)
371
372 rc = trexio_write_basis_prim_num(f, prim_num)
373 CALL trexio_error(rc)
374 END IF ! ionode
375
376 ! one-to-one mapping between shells and ...
377 ALLOCATE (nucleus_index(shell_num)) ! ...atomic indices
378 ALLOCATE (shell_ang_mom(shell_num)) ! ...angular momenta
379 ALLOCATE (shell_index(prim_num)) ! ...indices of primitive functions
380 ALLOCATE (exponents(prim_num)) ! ...primitive exponents
381 ALLOCATE (coefficients(prim_num)) ! ...contraction coefficients
382 ALLOCATE (prim_factor(prim_num)) ! ...primitive normalization factors
383
384 ishell = 0
385 ipgf = 0
386 DO iatom = 1, natoms
387 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
388 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
389 ! get the primary (orbital) basis set associated to this qs_kind
390 CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
391 ! get the info from the basis set
392 CALL get_gto_basis_set(basis_set, &
393 nset=nset, &
394 nshell=nshell, &
395 npgf=npgf, &
396 zet=zetas, &
397 gcc=gcc, &
398 l=l_shell_set)
399
400 DO iset = 1, nset
401 DO ishell_loc = 1, nshell(iset)
402 ishell = ishell + 1
403
404 ! nucleus_index array
405 nucleus_index(ishell) = iatom
406
407 ! shell_ang_mom array
408 l = l_shell_set(ishell_loc, iset)
409 shell_ang_mom(ishell) = l
410
411 ! shell_index array
412 shell_index(ipgf + 1:ipgf + npgf(iset)) = ishell
413
414 ! exponents array
415 exponents(ipgf + 1:ipgf + npgf(iset)) = zetas(1:npgf(iset), iset)
416
417 ! compute on the fly the normalization factor as in normalise_gcc_orb
418 ! and recover the original contraction coefficients to store them separately
419 expzet = 0.25_dp*real(2*l + 3, dp)
420 prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
421 DO i = 1, npgf(iset)
422 gcca = gcc(i, ishell_loc, iset)
423 zeta = zetas(i, iset)
424
425 ! primitives normalization factors array
426 prim_factor(i + ipgf) = prefac*zeta**expzet
427
428 ! contractio coefficients array
429 coefficients(i + ipgf) = gcca/prim_factor(i + ipgf)
430 END DO
431
432 ipgf = ipgf + npgf(iset)
433 END DO
434 END DO
435 END DO
436 ! just a failsafe check
437 cpassert(ishell == shell_num)
438 cpassert(ipgf == prim_num)
439
440 IF (ionode) THEN
441 rc = trexio_write_basis_nucleus_index(f, nucleus_index)
442 CALL trexio_error(rc)
443
444 rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
445 CALL trexio_error(rc)
446
447 ! Normalization factors are shoved in the AO group
448 ALLOCATE (shell_factor(shell_num)) ! 1-to-1 map bw shells and normalization factors
449 shell_factor(:) = 1.0_dp
450 rc = trexio_write_basis_shell_factor(f, shell_factor)
451 CALL trexio_error(rc)
452 DEALLOCATE (shell_factor)
453
454 ! This is always 0 for Gaussian basis sets
455 ALLOCATE (r_power(shell_num)) ! 1-to-1 map bw shells radial function powers
456 r_power(:) = 0
457 rc = trexio_write_basis_r_power(f, r_power)
458 CALL trexio_error(rc)
459 DEALLOCATE (r_power)
460
461 rc = trexio_write_basis_shell_index(f, shell_index)
462 CALL trexio_error(rc)
463
464 rc = trexio_write_basis_exponent(f, exponents)
465 CALL trexio_error(rc)
466
467 rc = trexio_write_basis_coefficient(f, coefficients)
468 CALL trexio_error(rc)
469
470 ! Normalization factors are shoved in the AO group
471 rc = trexio_write_basis_prim_factor(f, prim_factor)
472 CALL trexio_error(rc)
473 END IF
474
475 DEALLOCATE (nucleus_index)
476 DEALLOCATE (shell_index)
477 DEALLOCATE (exponents)
478 DEALLOCATE (coefficients)
479 DEALLOCATE (prim_factor)
480 ! shell_ang_mom is needed in the MO group, so will be deallocated there
481
482 !========================================================================================!
483 ! ECP group
484 !========================================================================================!
485 IF (ionode) THEN
486 CALL get_qs_kind_set(kind_set, sgp_potential_present=sgp_potential_present)
487
488 ! figure out whether we actually have ECP potentials
489 ecp_num = 0
490 IF (sgp_potential_present) THEN
491 DO iatom = 1, natoms
492 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
493 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
494 ! get the the sgp_potential associated to this qs_kind
495 CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential)
496
497 ! get the info on the potential
498 IF (ASSOCIATED(sgp_potential)) THEN
499 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
500 IF (ecp_local) THEN
501 ! get number of local terms
502 CALL get_potential(potential=sgp_potential, nloc=nloc)
503 ecp_num = ecp_num + nloc
504 END IF
505 IF (ecp_semi_local) THEN
506 ! get number of semilocal terms
507 CALL get_potential(potential=sgp_potential, npot=npot)
508 ecp_num = ecp_num + sum(npot)
509 END IF
510 END IF
511 END DO
512 END IF
513
514 ! if we have ECP potentials, populate the ECP group
515 IF (ecp_num > 0) THEN
516 ALLOCATE (z_core(natoms))
517 ALLOCATE (max_ang_mom_plus_1(natoms))
518 max_ang_mom_plus_1(:) = 0
519
520 ALLOCATE (ang_mom(ecp_num))
521 ALLOCATE (nucleus_index(ecp_num))
522 ALLOCATE (exponents(ecp_num))
523 ALLOCATE (coefficients(ecp_num))
524 ALLOCATE (powers(ecp_num))
525
526 iecp = 0
527 DO iatom = 1, natoms
528 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
529 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind, z=z)
530 ! get the the sgp_potential associated to this qs_kind
531 CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential, zeff=zeff)
532
533 ! number of core electrons removed by the ECP
534 z_core(iatom) = z - int(zeff)
535
536 ! get the info on the potential
537 IF (ASSOCIATED(sgp_potential)) THEN
538 CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
539
540 ! deal with the local part
541 IF (ecp_local) THEN
542 CALL get_potential(potential=sgp_potential, nloc=nloc, sl_lmax=sl_lmax)
543 ang_mom(iecp + 1:iecp + nloc) = sl_lmax + 1
544 nucleus_index(iecp + 1:iecp + nloc) = iatom
545 exponents(iecp + 1:iecp + nloc) = sgp_potential%bloc(1:nloc)
546 coefficients(iecp + 1:iecp + nloc) = sgp_potential%aloc(1:nloc)
547 powers(iecp + 1:iecp + nloc) = sgp_potential%nrloc(1:nloc) - 2
548 iecp = iecp + nloc
549 END IF
550
551 ! deal with the semilocal part
552 IF (ecp_semi_local) THEN
553 CALL get_potential(potential=sgp_potential, npot=npot, sl_lmax=sl_lmax)
554 max_ang_mom_plus_1(iatom) = sl_lmax + 1
555
556 DO sl_l = 0, sl_lmax
557 nsemiloc = npot(sl_l)
558 ang_mom(iecp + 1:iecp + nsemiloc) = sl_l
559 nucleus_index(iecp + 1:iecp + nsemiloc) = iatom
560 exponents(iecp + 1:iecp + nsemiloc) = sgp_potential%bpot(1:nsemiloc, sl_l)
561 coefficients(iecp + 1:iecp + nsemiloc) = sgp_potential%apot(1:nsemiloc, sl_l)
562 powers(iecp + 1:iecp + nsemiloc) = sgp_potential%nrpot(1:nsemiloc, sl_l) - 2
563 iecp = iecp + nsemiloc
564 END DO
565 END IF
566 END IF
567 END DO
568
569 ! fail-safe check
570 cpassert(iecp == ecp_num)
571
572 rc = trexio_write_ecp_num(f, ecp_num)
573 CALL trexio_error(rc)
574
575 rc = trexio_write_ecp_z_core(f, z_core)
576 CALL trexio_error(rc)
577 DEALLOCATE (z_core)
578
579 rc = trexio_write_ecp_max_ang_mom_plus_1(f, max_ang_mom_plus_1)
580 CALL trexio_error(rc)
581 DEALLOCATE (max_ang_mom_plus_1)
582
583 rc = trexio_write_ecp_ang_mom(f, ang_mom)
584 CALL trexio_error(rc)
585 DEALLOCATE (ang_mom)
586
587 rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
588 CALL trexio_error(rc)
589 DEALLOCATE (nucleus_index)
590
591 rc = trexio_write_ecp_exponent(f, exponents)
592 CALL trexio_error(rc)
593 DEALLOCATE (exponents)
594
595 rc = trexio_write_ecp_coefficient(f, coefficients)
596 CALL trexio_error(rc)
597 DEALLOCATE (coefficients)
598
599 rc = trexio_write_ecp_power(f, powers)
600 CALL trexio_error(rc)
601 DEALLOCATE (powers)
602 END IF
603
604 END IF ! ionode
605
606 !========================================================================================!
607 ! Grid group
608 !========================================================================================!
609 ! TODO
610
611 !========================================================================================!
612 ! AO group
613 !========================================================================================!
614 CALL get_qs_env(qs_env, qs_kind_set=kind_set)
615 CALL get_qs_kind_set(kind_set, ncgf=ncgf, nsgf=nsgf)
616
617 CALL section_vals_val_get(trexio_section, "CARTESIAN", l_val=save_cartesian)
618 IF (save_cartesian) THEN
619 ao_num = ncgf
620 ELSE
621 ao_num = nsgf
622 END IF
623
624 IF (ionode) THEN
625 IF (save_cartesian) THEN
626 rc = trexio_write_ao_cartesian(f, 1)
627 ELSE
628 rc = trexio_write_ao_cartesian(f, 0)
629 END IF
630 CALL trexio_error(rc)
631
632 rc = trexio_write_ao_num(f, ao_num)
633 CALL trexio_error(rc)
634 END IF
635
636 ! one-to-one mapping between AOs and ...
637 ALLOCATE (ao_shell(ao_num)) ! ..shells
638 ALLOCATE (ao_normalization(ao_num)) ! ..normalization factors
639
640 ! we need to be consistent with the basis group on the shell indices
641 ishell = 0 ! global shell index
642 igf = 0 ! global AO index
643 DO iatom = 1, natoms
644 ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
645 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
646 ! get the primary (orbital) basis set associated to this qs_kind
647 CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
648 ! get the info from the basis set
649 CALL get_gto_basis_set(basis_set, &
650 nset=nset, &
651 nshell=nshell, &
652 norm_cgf=norm_cgf, &
653 ncgf=ncgf, &
654 nsgf=nsgf, &
655 l=l_shell_set)
656
657 icgf = 0
658 DO iset = 1, nset
659 DO ishell_loc = 1, nshell(iset)
660 ! global shell index
661 ishell = ishell + 1
662 ! angular momentum l of this shell
663 lshell = l_shell_set(ishell_loc, iset)
664
665 ! number of AOs in this shell
666 IF (save_cartesian) THEN
667 ngf_shell = nco(lshell)
668 ELSE
669 ngf_shell = nso(lshell)
670 END IF
671
672 ! one-to-one mapping between AOs and shells
673 ao_shell(igf + 1:igf + ngf_shell) = ishell
674
675 ! one-to-one mapping between AOs and normalization factors
676 IF (save_cartesian) THEN
677 ao_normalization(igf + 1:igf + ngf_shell) = norm_cgf(icgf + 1:icgf + ngf_shell)
678 ELSE
679 ! allocate some temporary arrays
680 ALLOCATE (diag_ncgf(nco(lshell), nco(lshell)))
681 ALLOCATE (diag_nsgf(nso(lshell), nso(lshell)))
682 ALLOCATE (temp(nso(lshell), nco(lshell)))
683 diag_ncgf = 0.0_dp
684 diag_nsgf = 0.0_dp
685 temp = 0.0_dp
686
687 DO i = 1, nco(lshell)
688 diag_ncgf(i, i) = norm_cgf(icgf + i)
689 END DO
690
691 ! transform the normalization factors from Cartesian to solid harmonics
692 temp(:, :) = matmul(orbtramat(lshell)%c2s, diag_ncgf)
693 diag_nsgf(:, :) = matmul(temp, transpose(orbtramat(lshell)%s2c))
694 DO i = 1, nso(lshell)
695 ao_normalization(igf + i) = diag_nsgf(i, i)
696 END DO
697
698 DEALLOCATE (diag_ncgf)
699 DEALLOCATE (diag_nsgf)
700 DEALLOCATE (temp)
701 END IF
702
703 igf = igf + ngf_shell
704 icgf = icgf + nco(lshell)
705 END DO
706 END DO
707 ! just a failsafe check
708 cpassert(icgf == ncgf)
709 END DO
710
711 IF (ionode) THEN
712 rc = trexio_write_ao_shell(f, ao_shell)
713 CALL trexio_error(rc)
714
715 rc = trexio_write_ao_normalization(f, ao_normalization)
716 CALL trexio_error(rc)
717 END IF
718
719 DEALLOCATE (ao_shell)
720 DEALLOCATE (ao_normalization)
721
722 !========================================================================================!
723 ! MO group
724 !========================================================================================!
725 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints, dft_control=dft_control, &
726 particle_set=particle_set, qs_kind_set=kind_set, blacs_env=blacs_env)
727 nspins = dft_control%nspins
728 CALL get_qs_kind_set(kind_set, nsgf=nsgf)
729 nmo_spin = 0
730
731 ! figure out that total number of MOs
732 mo_num = 0
733 IF (do_kpoints) THEN
734 CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, use_real_wfn=use_real_wfn)
735 CALL get_kpoint_env(kp_env(1)%kpoint_env, mos=mos_kp)
736 DO ispin = 1, nspins
737 CALL get_mo_set(mos_kp(1, ispin), nmo=nmo)
738 nmo_spin(ispin) = nmo
739 END DO
740 mo_num = nkp*sum(nmo_spin)
741
742 ! we create a distributed fm matrix to gather the MOs from everywhere (in sph basis)
743 CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
744 nrow_global=nsgf, ncol_global=mo_num)
745 CALL cp_fm_create(fm_mo_coeff, fm_struct)
746 CALL cp_fm_set_all(fm_mo_coeff, 0.0_dp)
747 IF (.NOT. use_real_wfn) THEN
748 CALL cp_fm_create(fm_mo_coeff_im, fm_struct)
749 CALL cp_fm_set_all(fm_mo_coeff_im, 0.0_dp)
750 END IF
751 CALL cp_fm_struct_release(fm_struct)
752 ELSE
753 CALL get_qs_env(qs_env, mos=mos)
754 DO ispin = 1, nspins
755 CALL get_mo_set(mos(ispin), nmo=nmo)
756 nmo_spin(ispin) = nmo
757 END DO
758 mo_num = sum(nmo_spin)
759 END IF
760
761 ! allocate all the arrays
762 ALLOCATE (mo_coefficient(ao_num, mo_num))
763 mo_coefficient(:, :) = 0.0_dp
764 ALLOCATE (mo_energy(mo_num))
765 mo_energy(:) = 0.0_dp
766 ALLOCATE (mo_occupation(mo_num))
767 mo_occupation(:) = 0.0_dp
768 ALLOCATE (mo_spin(mo_num))
769 mo_spin(:) = 0
770 ! extra arrays for kpoints
771 IF (do_kpoints) THEN
772 ALLOCATE (mo_coefficient_im(ao_num, mo_num))
773 mo_coefficient_im(:, :) = 0.0_dp
774 ALLOCATE (mo_kpoint(mo_num))
775 mo_kpoint(:) = 0
776 END IF
777
778 ! in case of kpoints, we do this in 2 steps:
779 ! 1. we gather the MOs of each kpt and pipe them into a single large distributed fm matrix;
780 ! 2. we possibly transform the MOs of each kpt to Cartesian AOs and write them in the single large local array;
781 IF (do_kpoints) THEN
782 CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, kp_range=kp_range)
783
784 DO ispin = 1, nspins
785 DO ikp = 1, nkp
786 nmo = nmo_spin(ispin)
787 ! global index to store the MOs
788 imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
789
790 ! do I have this kpoint on this rank?
791 IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
792 ikp_loc = ikp - kp_range(1) + 1
793 ! get the mo set for this kpoint
794 CALL get_kpoint_env(kp_env(ikp_loc)%kpoint_env, mos=mos_kp)
795
796 ! if MOs are stored with dbcsr, copy them to fm
797 IF (mos_kp(1, ispin)%use_mo_coeff_b) THEN
798 CALL copy_dbcsr_to_fm(mos_kp(1, ispin)%mo_coeff_b, mos_kp(1, ispin)%mo_coeff)
799 END IF
800 ! copy real part of MO coefficients to large distributed fm matrix
801 CALL cp_fm_to_fm_submat_general(mos_kp(1, ispin)%mo_coeff, fm_mo_coeff, &
802 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
803
804 ! copy MO energies to local arrays
805 mo_energy(imo + 1:imo + nmo) = mos_kp(1, ispin)%eigenvalues(1:nmo)
806
807 ! copy MO occupations to local arrays
808 mo_occupation(imo + 1:imo + nmo) = mos_kp(1, ispin)%occupation_numbers(1:nmo)
809
810 ! same for the imaginary part of MO coefficients
811 IF (.NOT. use_real_wfn) THEN
812 IF (mos_kp(2, ispin)%use_mo_coeff_b) THEN
813 CALL copy_dbcsr_to_fm(mos_kp(2, ispin)%mo_coeff_b, mos_kp(2, ispin)%mo_coeff)
814 END IF
815 CALL cp_fm_to_fm_submat_general(mos_kp(2, ispin)%mo_coeff, fm_mo_coeff_im, &
816 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
817 END IF
818 ELSE
819 ! call with a dummy fm for receiving the data
820 CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff, &
821 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
822 IF (.NOT. use_real_wfn) THEN
823 CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff_im, &
824 nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
825 END IF
826 END IF
827 END DO
828 END DO
829 END IF
830
831 ! reduce MO energies and occupations to the master node
832 IF (do_kpoints) THEN
833 CALL get_kpoint_info(kpoints, para_env_inter_kp=para_env_inter_kp)
834 CALL para_env_inter_kp%sum(mo_energy)
835 CALL para_env_inter_kp%sum(mo_occupation)
836 END IF
837
838 ! AO order map from CP2K to TREXIO convention
839 ! from m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
840 ! to m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
841 ALLOCATE (cp2k_to_trexio_ang_mom(nsgf))
842 i = 0
843 DO ishell = 1, shell_num
844 l = shell_ang_mom(ishell)
845 DO k = 1, 2*l + 1
846 m = (-1)**k*floor(real(k, kind=dp)/2.0_dp)
847 cp2k_to_trexio_ang_mom(i + k) = i + l + 1 + m
848 END DO
849 i = i + 2*l + 1
850 END DO
851 cpassert(i == nsgf)
852
853 ! second step: here we actually put everything in the local arrays for writing to trexio
854 DO ispin = 1, nspins
855 ! get number of MOs for this spin
856 nmo = nmo_spin(ispin)
857 ! allocate local temp array to transform the MOs of each kpoint/spin
858 ALLOCATE (mos_sgf(nsgf, nmo))
859
860 IF (do_kpoints) THEN
861 DO ikp = 1, nkp
862 ! global index to store the MOs
863 imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
864
865 ! store kpoint index
866 mo_kpoint(imo + 1:imo + nmo) = ikp
867 ! store the MO spins
868 mo_spin(imo + 1:imo + nmo) = ispin - 1
869
870 ! transform and store the MO coefficients
871 CALL cp_fm_get_submatrix(fm_mo_coeff, mos_sgf, 1, imo + 1, nsgf, nmo)
872 IF (save_cartesian) THEN
873 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
874 ELSE
875 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
876 DO i = 1, nsgf
877 mo_coefficient(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
878 END DO
879 END IF
880
881 ! we have to do it for the imaginary part as well
882 IF (.NOT. use_real_wfn) THEN
883 CALL cp_fm_get_submatrix(fm_mo_coeff_im, mos_sgf, 1, imo + 1, nsgf, nmo)
884 IF (save_cartesian) THEN
885 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient_im(:, imo + 1:imo + nmo))
886 ELSE
887 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
888 DO i = 1, nsgf
889 mo_coefficient_im(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
890 END DO
891 END IF
892 END IF
893 END DO
894 ELSE ! no k-points
895 ! global index to store the MOs
896 imo = (ispin - 1)*nmo_spin(1)
897 ! store the MO energies
898 mo_energy(imo + 1:imo + nmo) = mos(ispin)%eigenvalues
899 ! store the MO occupations
900 mo_occupation(imo + 1:imo + nmo) = mos(ispin)%occupation_numbers
901 ! store the MO spins
902 mo_spin(imo + 1:imo + nmo) = ispin - 1
903
904 ! check if we are using the dbcsr mo_coeff and copy them to fm if needed
905 IF (mos(ispin)%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff)
906
907 ! allocate a normal fortran array to store the spherical MO coefficients
908 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, mos_sgf)
909
910 IF (save_cartesian) THEN
911 CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
912 ELSE
913 ! we have to reorder the MOs since CP2K and TREXIO have different conventions
914 DO i = 1, nsgf
915 mo_coefficient(i, imo + 1:imo + nmo) = mos_sgf(cp2k_to_trexio_ang_mom(i), :)
916 END DO
917 END IF
918 END IF
919
920 DEALLOCATE (mos_sgf)
921 END DO
922
923 IF (ionode) THEN
924 rc = trexio_write_mo_type(f, 'Canonical', len_trim('Canonical') + 1)
925 CALL trexio_error(rc)
926
927 rc = trexio_write_mo_num(f, mo_num)
928 CALL trexio_error(rc)
929
930 rc = trexio_write_mo_coefficient(f, mo_coefficient)
931 CALL trexio_error(rc)
932
933 rc = trexio_write_mo_energy(f, mo_energy)
934 CALL trexio_error(rc)
935
936 rc = trexio_write_mo_occupation(f, mo_occupation)
937 CALL trexio_error(rc)
938
939 rc = trexio_write_mo_spin(f, mo_spin)
940 CALL trexio_error(rc)
941
942 IF (do_kpoints) THEN
943 rc = trexio_write_mo_coefficient_im(f, mo_coefficient_im)
944 CALL trexio_error(rc)
945
946 rc = trexio_write_mo_k_point(f, mo_kpoint)
947 CALL trexio_error(rc)
948 END IF
949 END IF
950
951 DEALLOCATE (mo_coefficient)
952 DEALLOCATE (mo_energy)
953 DEALLOCATE (mo_occupation)
954 DEALLOCATE (mo_spin)
955 IF (do_kpoints) THEN
956 DEALLOCATE (mo_coefficient_im)
957 DEALLOCATE (mo_kpoint)
958 CALL cp_fm_release(fm_mo_coeff)
959 CALL cp_fm_release(fm_mo_coeff_im)
960 END IF
961
962 !========================================================================================!
963 ! RDM group
964 !========================================================================================!
965 !TODO
966
967 !========================================================================================!
968 ! Energy derivative group
969 !========================================================================================!
970 IF (PRESENT(energy_derivative)) THEN
971 filename_de = trim(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
972
973 ALLOCATE (dedp(nsgf, nsgf))
974 dedp(:, :) = 0.0_dp
975
976 DO ispin = 1, nspins
977 CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
978 DO WHILE (dbcsr_iterator_blocks_left(iter))
979 ! the offsets tell me the global index of the matrix, not the index of the block
980 CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
981 row_size=row_size, col_size=col_size, &
982 row_offset=row_offset, col_offset=col_offset)
983
984 ! Copy data from block to array
985 DO i = 1, row_size
986 DO j = 1, col_size
987 dedp(row_offset + i - 1, col_offset + j - 1) = data_block(i, j)
988 END DO
989 END DO
990 END DO
991 CALL dbcsr_iterator_stop(iter)
992
993 ! symmetrize the matrix if needed
994 SELECT CASE (dbcsr_get_matrix_type(energy_derivative(ispin)%matrix))
995 CASE (dbcsr_type_symmetric)
996 CALL symmetrize_matrix(dedp, "upper_to_lower")
997 CASE (dbcsr_type_antisymmetric)
998 CALL symmetrize_matrix(dedp, "anti_upper_to_lower")
999 CASE (dbcsr_type_no_symmetry)
1000 CASE DEFAULT
1001 cpabort("Unknown matrix type for energy derivative")
1002 END SELECT
1003 END DO
1004
1005 ! reduce the dEdP matrix to the master node
1006 CALL para_env%sum(dedp)
1007
1008 ! print the dEdP matrix to a file
1009 IF (ionode) THEN
1010 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing derivative file ', trim(filename_de)
1011
1012 unit_de = 10
1013 CALL open_file(file_name=filename_de, &
1014 file_action="WRITE", &
1015 file_status="UNKNOWN", &
1016 unit_number=unit_de)
1017 WRITE (unit_de, '(I0, 1X, I0)') nsgf, nsgf
1018 DO i = 1, nsgf
1019 WRITE (unit_de, '(*(1X, F15.8))') (dedp(cp2k_to_trexio_ang_mom(i), &
1020 cp2k_to_trexio_ang_mom(j)), j=1, nsgf)
1021 END DO
1022 CALL close_file(unit_number=unit_de)
1023 END IF
1024
1025 DEALLOCATE (dedp)
1026 END IF
1027
1028 ! Ddeallocate arrays used throughout the subroutine
1029 DEALLOCATE (shell_ang_mom)
1030 DEALLOCATE (cp2k_to_trexio_ang_mom)
1031
1032 !========================================================================================!
1033 ! Close the TREXIO file
1034 !========================================================================================!
1035 IF (ionode) THEN
1036 rc = trexio_close(f)
1037 CALL trexio_error(rc)
1038 END IF
1039#else
1040 mark_used(qs_env)
1041 mark_used(trexio_section)
1042 mark_used(energy_derivative)
1043 cpwarn('TREXIO support has not been enabled in this build.')
1044#endif
1045 CALL timestop(handle)
1046
1047 END SUBROUTINE write_trexio
1048
1049! **************************************************************************************************
1050!> \brief Read a trexio file
1051!> \param qs_env the qs environment with all the info of the computation
1052!> \param trexio_filename the trexio filename without the extension
1053!> \param mo_set_trexio the MO set to read from the trexio file
1054!> \param energy_derivative the energy derivative to read from the trexio file
1055! **************************************************************************************************
1056 SUBROUTINE read_trexio(qs_env, trexio_filename, mo_set_trexio, energy_derivative)
1057 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1058 CHARACTER(len=*), INTENT(IN), OPTIONAL :: trexio_filename
1059 TYPE(mo_set_type), INTENT(OUT), DIMENSION(:), POINTER, OPTIONAL :: mo_set_trexio
1060 TYPE(dbcsr_p_type), INTENT(OUT), DIMENSION(:), POINTER, OPTIONAL :: energy_derivative
1061
1062 CHARACTER(LEN=*), PARAMETER :: routinen = 'read_trexio'
1063
1064 INTEGER :: handle
1065
1066#ifdef __TREXIO
1067 INTEGER :: output_unit, unit_de
1068 CHARACTER(len=default_path_length) :: filename, filename_de
1069 INTEGER(trexio_t) :: f ! The TREXIO file handle
1070 INTEGER(trexio_exit_code) :: rc ! TREXIO return code
1071
1072 LOGICAL :: ionode
1073
1074 CHARACTER(LEN=2) :: element_symbol
1075 CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: label
1076
1077 INTEGER :: ao_num, mo_num, nmo, nspins, ispin, nsgf, &
1078 save_cartesian, i, j, k, l, m, imo, ishell, &
1079 nshell, shell_num, nucleus_num, natoms, ikind, &
1080 iatom, nelectron, nrows, ncols, &
1081 row, col, row_size, col_size, &
1082 row_offset, col_offset, myprint
1083 INTEGER, DIMENSION(2) :: nmo_spin, electron_num
1084 INTEGER, DIMENSION(:), ALLOCATABLE :: mo_spin, shell_ang_mom, trexio_to_cp2k_ang_mom
1085
1086 REAL(kind=dp) :: zeff, maxocc
1087 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: mo_energy, mo_occupation, charge
1088 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: mo_coefficient, mos_sgf, coord, dedp, temp
1089 REAL(kind=dp), DIMENSION(:, :), POINTER :: data_block
1090
1091 TYPE(cp_logger_type), POINTER :: logger => null()
1092 TYPE(cp_fm_type), POINTER :: mo_coeff_ref, mo_coeff_target
1093 TYPE(mp_para_env_type), POINTER :: para_env => null()
1094 TYPE(dft_control_type), POINTER :: dft_control => null()
1095 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s => null()
1096 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set => null()
1097 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos => null()
1098 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => null()
1099 TYPE(dbcsr_iterator_type) :: iter
1100#endif
1101
1102 CALL timeset(routinen, handle)
1103
1104#ifdef __TREXIO
1105 logger => cp_get_default_logger()
1106 output_unit = cp_logger_get_default_io_unit(logger)
1107 myprint = logger%iter_info%print_level
1108
1109 cpassert(ASSOCIATED(qs_env))
1110
1111 ! get filename
1112 IF (.NOT. PRESENT(trexio_filename)) THEN
1113 filename = trim(logger%iter_info%project_name)//'-TREXIO.h5'
1114 filename_de = trim(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
1115 ELSE
1116 filename = trim(trexio_filename)//'.h5'
1117 filename_de = trim(trexio_filename)//'.dEdP.dat'
1118 END IF
1119
1120 CALL get_qs_env(qs_env, para_env=para_env)
1121 ionode = para_env%is_source()
1122
1123 ! Open the TREXIO file and check that we have the same molecule as in qs_env
1124 IF (ionode) THEN
1125 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Opening file named ', trim(filename)
1126 f = trexio_open(filename, 'r', trexio_hdf5, rc)
1127 CALL trexio_error(rc)
1128
1129 IF (myprint > medium_print_level) THEN
1130 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading molecule information...'
1131 END IF
1132 rc = trexio_read_nucleus_num(f, nucleus_num)
1133 CALL trexio_error(rc)
1134
1135 IF (myprint > medium_print_level) THEN
1136 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear coordinates...'
1137 END IF
1138 ALLOCATE (coord(3, nucleus_num))
1139 rc = trexio_read_nucleus_coord(f, coord)
1140 CALL trexio_error(rc)
1141
1142 IF (myprint > medium_print_level) THEN
1143 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear labels...'
1144 END IF
1145 ALLOCATE (label(nucleus_num))
1146 rc = trexio_read_nucleus_label(f, label, 3)
1147 CALL trexio_error(rc)
1148
1149 IF (myprint > medium_print_level) THEN
1150 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear charges...'
1151 END IF
1152 ALLOCATE (charge(nucleus_num))
1153 rc = trexio_read_nucleus_charge(f, charge)
1154 CALL trexio_error(rc)
1155
1156 ! get the same info from qs_env
1157 CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
1158
1159 ! check that we have the same number of atoms
1160 cpassert(nucleus_num == natoms)
1161
1162 DO iatom = 1, natoms
1163 ! compare the coordinates within a certain tolerance
1164 DO i = 1, 3
1165 cpassert(abs(coord(i, iatom) - particle_set(iatom)%r(i)) < 1.0e-6_dp)
1166 END DO
1167
1168 ! figure out the element symbol and to which kind_set entry this atomic_kind corresponds to
1169 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
1170 ! check that the element symbol is the same
1171 cpassert(trim(element_symbol) == trim(label(iatom)))
1172
1173 ! get the effective nuclear charge for this kind
1174 CALL get_qs_kind(kind_set(ikind), zeff=zeff)
1175 ! check that the nuclear charge is also the same
1176 cpassert(charge(iatom) == zeff)
1177 END DO
1178
1179 WRITE (output_unit, "((T2,A))") 'TREXIO| Molecule is the same as in qs_env'
1180 ! if we get here, we have the same molecule
1181 DEALLOCATE (coord)
1182 DEALLOCATE (label)
1183 DEALLOCATE (charge)
1184
1185 ! get info from trexio to map cp2k and trexio AOs
1186 rc = trexio_read_ao_cartesian(f, save_cartesian)
1187 CALL trexio_error(rc)
1188
1189 rc = trexio_read_ao_num(f, ao_num)
1190 CALL trexio_error(rc)
1191
1192 rc = trexio_read_basis_shell_num(f, shell_num)
1193 CALL trexio_error(rc)
1194 END IF
1195
1196 CALL para_env%bcast(save_cartesian, para_env%source)
1197 CALL para_env%bcast(ao_num, para_env%source)
1198 CALL para_env%bcast(shell_num, para_env%source)
1199
1200 IF (save_cartesian == 1) THEN
1201 cpabort('Reading Cartesian AOs is not yet supported.')
1202 END IF
1203
1204 ! check that the number of AOs and shells is the same
1205 CALL get_qs_env(qs_env, qs_kind_set=kind_set)
1206 CALL get_qs_kind_set(kind_set, nsgf=nsgf, nshell=nshell)
1207 cpassert(ao_num == nsgf)
1208 cpassert(shell_num == nshell)
1209
1210 ALLOCATE (shell_ang_mom(shell_num))
1211 shell_ang_mom(:) = 0
1212
1213 IF (ionode) THEN
1214 IF (myprint > medium_print_level) THEN
1215 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading shell angular momenta...'
1216 END IF
1217 rc = trexio_read_basis_shell_ang_mom(f, shell_ang_mom)
1218 CALL trexio_error(rc)
1219 END IF
1220
1221 CALL para_env%bcast(shell_ang_mom, para_env%source)
1222
1223 ! AO order map from TREXIO to CP2K convention
1224 ! from m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
1225 ! to m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
1226 ALLOCATE (trexio_to_cp2k_ang_mom(nsgf))
1227 i = 0
1228 DO ishell = 1, shell_num
1229 l = shell_ang_mom(ishell)
1230 DO k = 1, 2*l + 1
1231 m = (-1)**k*floor(real(k, kind=dp)/2.0_dp)
1232 trexio_to_cp2k_ang_mom(i + l + 1 + m) = i + k
1233 END DO
1234 i = i + 2*l + 1
1235 END DO
1236 cpassert(i == nsgf)
1237
1238 ! check whether we want to read MOs
1239 IF (PRESENT(mo_set_trexio)) THEN
1240 IF (output_unit > 1) THEN
1241 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading molecular orbitals...'
1242 END IF
1243
1244 ! at the moment, we assume that the basis set is the same
1245 ! first we read all arrays lengths we need from the trexio file
1246 IF (ionode) THEN
1247 rc = trexio_read_mo_num(f, mo_num)
1248 CALL trexio_error(rc)
1249
1250 rc = trexio_read_electron_up_num(f, electron_num(1))
1251 CALL trexio_error(rc)
1252
1253 rc = trexio_read_electron_dn_num(f, electron_num(2))
1254 CALL trexio_error(rc)
1255 END IF
1256
1257 ! broadcast information to all processors and allocate arrays
1258 CALL para_env%bcast(mo_num, para_env%source)
1259 CALL para_env%bcast(electron_num, para_env%source)
1260
1261 ! check that the number of MOs is the same
1262 CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1263 nspins = dft_control%nspins
1264 nmo_spin(:) = 0
1265 DO ispin = 1, nspins
1266 CALL get_mo_set(mos(ispin), nmo=nmo)
1267 nmo_spin(ispin) = nmo
1268 END DO
1269 cpassert(mo_num == sum(nmo_spin))
1270
1271 ALLOCATE (mo_coefficient(ao_num, mo_num))
1272 ALLOCATE (mo_energy(mo_num))
1273 ALLOCATE (mo_occupation(mo_num))
1274 ALLOCATE (mo_spin(mo_num))
1275
1276 mo_coefficient(:, :) = 0.0_dp
1277 mo_energy(:) = 0.0_dp
1278 mo_occupation(:) = 0.0_dp
1279 mo_spin(:) = 0
1280
1281 ! read the MOs info
1282 IF (ionode) THEN
1283 IF (myprint > medium_print_level) THEN
1284 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO coefficients...'
1285 END IF
1286 rc = trexio_read_mo_coefficient(f, mo_coefficient)
1287 CALL trexio_error(rc)
1288
1289 IF (myprint > medium_print_level) THEN
1290 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO energies...'
1291 END IF
1292 rc = trexio_read_mo_energy(f, mo_energy)
1293 CALL trexio_error(rc)
1294
1295 IF (myprint > medium_print_level) THEN
1296 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO occupations...'
1297 END IF
1298 rc = trexio_read_mo_occupation(f, mo_occupation)
1299 CALL trexio_error(rc)
1300
1301 IF (myprint > medium_print_level) THEN
1302 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO spins...'
1303 END IF
1304 rc = trexio_read_mo_spin(f, mo_spin)
1305 CALL trexio_error(rc)
1306 END IF
1307
1308 ! broadcast the data to all processors
1309 CALL para_env%bcast(mo_coefficient, para_env%source)
1310 CALL para_env%bcast(mo_energy, para_env%source)
1311 CALL para_env%bcast(mo_occupation, para_env%source)
1312 CALL para_env%bcast(mo_spin, para_env%source)
1313
1314 ! assume nspins and nmo_spin match the ones in the trexio file
1315 ! reorder magnetic quantum number
1316 DO ispin = 1, nspins
1317 ! global MOs index
1318 imo = (ispin - 1)*nmo_spin(1)
1319 ! get number of MOs for this spin
1320 nmo = nmo_spin(ispin)
1321 ! allocate local temp array to read MOs
1322 ALLOCATE (mos_sgf(nsgf, nmo))
1323 mos_sgf(:, :) = 0.0_dp
1324
1325 ! we need to reorder the MOs according to CP2K convention
1326 DO i = 1, nsgf
1327 mos_sgf(i, :) = mo_coefficient(trexio_to_cp2k_ang_mom(i), imo + 1:imo + nmo)
1328 END DO
1329
1330 IF (nspins == 1) THEN
1331 maxocc = 2.0_dp
1332 nelectron = electron_num(1) + electron_num(2)
1333 ELSE
1334 maxocc = 1.0_dp
1335 nelectron = electron_num(ispin)
1336 END IF
1337 ! the right number of active electrons per spin channel is initialized further down
1338 CALL allocate_mo_set(mo_set_trexio(ispin), nsgf, nmo, nelectron, 0.0_dp, maxocc, 0.0_dp)
1339
1340 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff_ref)
1341 CALL init_mo_set(mo_set_trexio(ispin), fm_ref=mo_coeff_ref, name="TREXIO MOs")
1342
1343 CALL get_mo_set(mo_set_trexio(ispin), mo_coeff=mo_coeff_target)
1344 DO j = 1, nmo
1345 ! make sure I copy the right spin channel
1346 cpassert(mo_spin(j) == ispin - 1)
1347 mo_set_trexio(ispin)%eigenvalues(j) = mo_energy(imo + j)
1348 mo_set_trexio(ispin)%occupation_numbers(j) = mo_occupation(imo + j)
1349 DO i = 1, nsgf
1350 CALL cp_fm_set_element(mo_coeff_target, i, j, mos_sgf(i, j))
1351 END DO
1352 END DO
1353
1354 DEALLOCATE (mos_sgf)
1355 END DO
1356
1357 DEALLOCATE (mo_coefficient)
1358 DEALLOCATE (mo_energy)
1359 DEALLOCATE (mo_occupation)
1360 DEALLOCATE (mo_spin)
1361
1362 END IF ! if MOs should be read
1363
1364 ! check whether we want to read derivatives
1365 IF (PRESENT(energy_derivative)) THEN
1366 IF (output_unit > 1) THEN
1367 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading energy derivatives...'
1368 END IF
1369
1370 ! Temporary solution: allocate here the energy derivatives matrix here
1371 ! assuming that nsgf is the same as the number read from the dEdP file
1372 ! TODO: once available in TREXIO, first read size and then allocate
1373 ! in the same way done for the MOs
1374 ALLOCATE (temp(nsgf, nsgf))
1375 temp(:, :) = 0.0_dp
1376
1377 ! check if file exists and open it
1378 IF (ionode) THEN
1379 IF (file_exists(filename_de)) THEN
1380 CALL open_file(file_name=filename_de, file_status="OLD", unit_number=unit_de)
1381 ELSE
1382 cpabort("Energy derivatives file "//trim(filename_de)//" not found")
1383 END IF
1384
1385 ! read the header and check everything is fine
1386 IF (myprint > medium_print_level) THEN
1387 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading header information...'
1388 END IF
1389 READ (unit_de, *) nrows, ncols
1390 IF (myprint > medium_print_level) THEN
1391 WRITE (output_unit, "((T2,A))") 'TREXIO| Check size of dEdP matrix...'
1392 END IF
1393 cpassert(nrows == nsgf)
1394 cpassert(ncols == nsgf)
1395
1396 ! read the data
1397 IF (myprint > medium_print_level) THEN
1398 WRITE (output_unit, "((T2,A))") 'TREXIO| Reading dEdP matrix...'
1399 END IF
1400 ! Read the data matrix
1401 DO i = 1, nrows
1402 READ (unit_de, *) (temp(i, j), j=1, ncols)
1403 END DO
1404
1405 CALL close_file(unit_number=unit_de)
1406 END IF
1407
1408 ! send data to all processes
1409 CALL para_env%bcast(temp, para_env%source)
1410
1411 ! Reshuffle
1412 ALLOCATE (dedp(nsgf, nsgf))
1413 dedp(:, :) = 0.0_dp
1414
1415 ! Reorder rows and columns according to trexio_to_cp2k_ang_mom mapping
1416 DO j = 1, nsgf
1417 DO i = 1, nsgf
1418 ! either this
1419 dedp(i, j) = temp(trexio_to_cp2k_ang_mom(i), trexio_to_cp2k_ang_mom(j))
1420 ! or this
1421 ! dEdP(cp2k_to_trexio_ang_mom(i), cp2k_to_trexio_ang_mom(j)) = temp(i, j)
1422 END DO
1423 END DO
1424
1425 DEALLOCATE (temp)
1426
1427 CALL get_qs_env(qs_env, matrix_s=matrix_s)
1428 DO ispin = 1, nspins
1429 ALLOCATE (energy_derivative(ispin)%matrix)
1430
1431 ! we use the overlap matrix as a template, copying it but removing the sparsity
1432 CALL dbcsr_copy(energy_derivative(ispin)%matrix, matrix_s(1)%matrix, &
1433 name='Energy Derivative', keep_sparsity=.false.)
1434 CALL dbcsr_set(energy_derivative(ispin)%matrix, 0.0_dp)
1435
1436 CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
1437 DO WHILE (dbcsr_iterator_blocks_left(iter))
1438 CALL dbcsr_iterator_next_block(iter, row, col, data_block, &
1439 row_size=row_size, col_size=col_size, &
1440 row_offset=row_offset, col_offset=col_offset)
1441
1442 ! Copy data from array to block
1443 DO i = 1, row_size
1444 DO j = 1, col_size
1445 data_block(i, j) = dedp(row_offset + i - 1, col_offset + j - 1)
1446 END DO
1447 END DO
1448 END DO
1449 CALL dbcsr_iterator_stop(iter)
1450 END DO
1451
1452 DEALLOCATE (dedp)
1453 END IF ! finished reading energy derivatives
1454
1455 ! Clean up
1456 IF (ALLOCATED(shell_ang_mom)) DEALLOCATE (shell_ang_mom)
1457 IF (ALLOCATED(trexio_to_cp2k_ang_mom)) DEALLOCATE (trexio_to_cp2k_ang_mom)
1458
1459 ! Close the TREXIO file
1460 IF (ionode) THEN
1461 WRITE (output_unit, "((T2,A,A))") 'TREXIO| Closing file named ', trim(filename)
1462 rc = trexio_close(f)
1463 CALL trexio_error(rc)
1464 END IF
1465
1466#else
1467 mark_used(qs_env)
1468 mark_used(trexio_filename)
1469 mark_used(mo_set_trexio)
1470 mark_used(energy_derivative)
1471 cpwarn('TREXIO support has not been enabled in this build.')
1472 cpabort('TREXIO Not Available')
1473#endif
1474 CALL timestop(handle)
1475
1476 END SUBROUTINE read_trexio
1477
1478#ifdef __TREXIO
1479! **************************************************************************************************
1480!> \brief Handles TREXIO errors
1481!> \param rc the TREXIO return code
1482! **************************************************************************************************
1483 SUBROUTINE trexio_error(rc)
1484 INTEGER(trexio_exit_code), INTENT(IN) :: rc
1485
1486 CHARACTER(LEN=128) :: err_msg
1487
1488 IF (rc /= trexio_success) THEN
1489 CALL trexio_string_of_error(rc, err_msg)
1490 cpabort('TREXIO Error: '//trim(err_msg))
1491 END IF
1492
1493 END SUBROUTINE trexio_error
1494
1495! **************************************************************************************************
1496!> \brief Computes the nuclear repulsion energy of a molecular system
1497!> \param particle_set the set of particles in the system
1498!> \param kind_set the set of qs_kinds in the system
1499!> \param e_nn the nuclear repulsion energy
1500! **************************************************************************************************
1501 SUBROUTINE nuclear_repulsion_energy(particle_set, kind_set, e_nn)
1502 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1503 POINTER :: particle_set
1504 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1505 POINTER :: kind_set
1506 REAL(kind=dp), INTENT(OUT) :: e_nn
1507
1508 INTEGER :: i, ikind, j, jkind, natoms
1509 REAL(kind=dp) :: r_ij, zeff_i, zeff_j
1510
1511 natoms = SIZE(particle_set)
1512 e_nn = 0.0_dp
1513 DO i = 1, natoms
1514 CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
1515 CALL get_qs_kind(kind_set(ikind), zeff=zeff_i)
1516 DO j = i + 1, natoms
1517 r_ij = norm2(particle_set(i)%r - particle_set(j)%r)
1518
1519 CALL get_atomic_kind(particle_set(j)%atomic_kind, kind_number=jkind)
1520 CALL get_qs_kind(kind_set(jkind), zeff=zeff_j)
1521
1522 e_nn = e_nn + zeff_i*zeff_j/r_ij
1523 END DO
1524 END DO
1525
1526 END SUBROUTINE nuclear_repulsion_energy
1527
1528! **************************************************************************************************
1529!> \brief Computes a spherical to cartesian MO transformation (solid harmonics in reality)
1530!> \param mos_sgf the MO coefficients in spherical AO basis
1531!> \param particle_set the set of particles in the system
1532!> \param qs_kind_set the set of qs_kinds in the system
1533!> \param mos_cgf the transformed MO coefficients in Cartesian AO basis
1534! **************************************************************************************************
1535 SUBROUTINE spherical_to_cartesian_mo(mos_sgf, particle_set, qs_kind_set, mos_cgf)
1536 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: mos_sgf
1537 TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1538 POINTER :: particle_set
1539 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1540 POINTER :: qs_kind_set
1541 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: mos_cgf
1542
1543 INTEGER :: iatom, icgf, ikind, iset, isgf, ishell, &
1544 lshell, ncgf, nmo, nset, nsgf
1545 INTEGER, DIMENSION(:), POINTER :: nshell
1546 INTEGER, DIMENSION(:, :), POINTER :: l
1547 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1548
1549 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1550
1551 mos_cgf = 0.0_dp
1552 nmo = SIZE(mos_sgf, 2)
1553
1554 ! Transform spherical MOs to Cartesian MOs
1555 icgf = 1
1556 isgf = 1
1557 DO iatom = 1, SIZE(particle_set)
1558 NULLIFY (orb_basis_set)
1559 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1560 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1561
1562 IF (ASSOCIATED(orb_basis_set)) THEN
1563 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1564 nset=nset, &
1565 nshell=nshell, &
1566 l=l)
1567 DO iset = 1, nset
1568 DO ishell = 1, nshell(iset)
1569 lshell = l(ishell, iset)
1570 CALL dgemm("T", "N", nco(lshell), nmo, nso(lshell), 1.0_dp, &
1571 orbtramat(lshell)%c2s, nso(lshell), &
1572 mos_sgf(isgf, 1), nsgf, 0.0_dp, &
1573 mos_cgf(icgf, 1), ncgf)
1574 icgf = icgf + nco(lshell)
1575 isgf = isgf + nso(lshell)
1576 END DO
1577 END DO
1578 ELSE
1579 ! assume atom without basis set
1580 cpabort("Unknown basis set type")
1581 END IF
1582 END DO ! iatom
1583
1584 END SUBROUTINE spherical_to_cartesian_mo
1585#endif
1586
1587END 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:206
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
external_potential_types::get_potential
Definition external_potential_types.F:276
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:138
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)
...
Definition basis_set_types.F:624
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp2k_info
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:16
cp2k_info::cp2k_version
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:42
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
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_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)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:847
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
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:295
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)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
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)
...
Definition cp_dbcsr_output.F:163
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:501
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:351
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:1009
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:2002
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:528
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:894
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:693
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
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:515
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:73
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:704
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:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
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:779
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)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
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
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:1197
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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, init_u_ramping_each_scf, reltmat, ghost, 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:454
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, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
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:206
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)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:245
trexio_utils
The module to read/write TREX IO files for interfacing CP2K with other programs.
Definition trexio_utils.F:14
trexio_utils::write_trexio
subroutine, public write_trexio(qs_env, trexio_section, energy_derivative)
Write a trexio file.
Definition trexio_utils.F:123
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:1057
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:55
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:570
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
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:113
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:172
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:85
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
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:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46