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