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