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 : !> \author Stefano Battaglia
13 : ! **************************************************************************************************
14 : MODULE trexio_utils
15 :
16 : USE ai_onecenter, ONLY: sg_overlap
17 : USE atomic_kind_types, ONLY: get_atomic_kind
18 : USE basis_set_types, ONLY: gto_basis_set_type, get_gto_basis_set
19 : USE cell_types, ONLY: cell_type
20 : USE cp2k_info, ONLY: cp2k_version
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
24 : USE cp_files, ONLY: close_file, file_exists, open_file
25 : USE cp_fm_types, ONLY: cp_fm_get_info, cp_fm_type, cp_fm_create, cp_fm_set_all, &
26 : cp_fm_get_submatrix, cp_fm_to_fm_submat_general, cp_fm_release, &
27 : cp_fm_set_element
28 : USE cp_fm_struct, ONLY: cp_fm_struct_create, &
29 : cp_fm_struct_release, &
30 : cp_fm_struct_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger, &
32 : cp_logger_get_default_io_unit, &
33 : cp_logger_type
34 : USE cp_dbcsr_api, ONLY: dbcsr_p_type, dbcsr_iterator_type, dbcsr_iterator_start, &
35 : dbcsr_iterator_stop, dbcsr_iterator_blocks_left, &
36 : dbcsr_iterator_next_block, dbcsr_copy, dbcsr_set, &
37 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
38 : dbcsr_type_symmetric, dbcsr_get_matrix_type
39 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
40 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
41 : USE cp_output_handling, ONLY: medium_print_level
42 : USE external_potential_types, ONLY: sgp_potential_type, get_potential
43 : USE input_section_types, ONLY: section_vals_type, section_vals_get, &
44 : section_vals_val_get
45 : USE kinds, ONLY: default_path_length, dp
46 : USE kpoint_types, ONLY: kpoint_type, kpoint_env_p_type, &
47 : get_kpoint_info, get_kpoint_env
48 : USE mathconstants, ONLY: fourpi, pi, fac
49 : USE mathlib, ONLY: symmetrize_matrix
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE orbital_pointers, ONLY: nco, nso
52 : USE orbital_transformation_matrices, ONLY: orbtramat
53 : USE particle_types, ONLY: particle_type
54 : USE qs_energy_types, ONLY: qs_energy_type
55 : USE qs_environment_types, ONLY: get_qs_env, &
56 : qs_environment_type
57 : USE qs_kind_types, ONLY: get_qs_kind, get_qs_kind_set, &
58 : qs_kind_type
59 : USE qs_mo_types, ONLY: mo_set_type, get_mo_set, init_mo_set, allocate_mo_set
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 :
115 : CONTAINS
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 8 : 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 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
146 : TYPE(qs_energy_type), POINTER :: energy
147 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
148 : TYPE(sgp_potential_type), POINTER :: sgp_potential
149 8 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
150 8 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
151 8 : 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 8 : 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 8 : INTEGER, DIMENSION(:), ALLOCATABLE :: nucleus_index, shell_ang_mom, r_power, &
172 8 : shell_index, z_core, max_ang_mom_plus_1, &
173 8 : ang_mom, powers, ao_shell, mo_spin, mo_kpoint, &
174 8 : cp2k_to_trexio_ang_mom
175 8 : INTEGER, DIMENSION(:), POINTER :: nshell, npgf
176 8 : INTEGER, DIMENSION(:, :), POINTER :: l_shell_set
177 8 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: charge, shell_factor, exponents, coefficients, &
178 8 : prim_factor, ao_normalization, mo_energy, &
179 8 : mo_occupation, Sgcc, ecp_coefficients, &
180 8 : sgf_coefficients
181 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp, norm_cgf
182 8 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: coord, mo_coefficient, mo_coefficient_im, &
183 8 : mos_sgf, dEdP, Sloc
184 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zetas, data_block, xkp
185 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
186 :
187 8 : CALL timeset(routineN, handle)
188 :
189 8 : NULLIFY (cell, logger, dft_control, basis_set, kpoints, particle_set, energy, kind_set)
190 8 : NULLIFY (sgp_potential, mos, mos_kp, kp_env, para_env, para_env_inter_kp, blacs_env)
191 8 : NULLIFY (fm_struct, nshell, npgf, l_shell_set, wkp, norm_cgf, zetas, data_block, gcc)
192 :
193 8 : logger => cp_get_default_logger()
194 8 : output_unit = cp_logger_get_default_io_unit(logger)
195 :
196 8 : CPASSERT(ASSOCIATED(qs_env))
197 :
198 : ! get filename
199 8 : CALL section_vals_val_get(trexio_section, "FILENAME", c_val=filename, explicit=explicit)
200 8 : IF (.NOT. explicit) THEN
201 6 : filename = TRIM(logger%iter_info%project_name)//'-TREXIO.h5'
202 : ELSE
203 2 : filename = TRIM(filename)//'.h5'
204 : END IF
205 :
206 8 : CALL get_qs_env(qs_env, para_env=para_env)
207 8 : ionode = para_env%is_source()
208 :
209 : ! inquire whether a file with the same name already exists, if yes, delete it
210 8 : IF (ionode) THEN
211 4 : IF (file_exists(filename)) THEN
212 0 : CALL open_file(filename, unit_number=unit_trexio)
213 0 : CALL close_file(unit_number=unit_trexio, file_status="DELETE")
214 : END IF
215 :
216 : !========================================================================================!
217 : ! Open the TREXIO file
218 : !========================================================================================!
219 4 : WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing trexio file ', TRIM(filename)
220 4 : f = trexio_open(filename, 'w', TREXIO_HDF5, rc)
221 4 : CALL trexio_error(rc)
222 :
223 : !========================================================================================!
224 : ! Metadata group
225 : !========================================================================================!
226 4 : rc = trexio_write_metadata_code_num(f, 1)
227 4 : CALL trexio_error(rc)
228 :
229 4 : rc = trexio_write_metadata_code(f, cp2k_version, LEN_TRIM(cp2k_version) + 1)
230 4 : CALL trexio_error(rc)
231 :
232 : !========================================================================================!
233 : ! Nucleus group
234 : !========================================================================================!
235 4 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
236 :
237 4 : rc = trexio_write_nucleus_num(f, natoms)
238 4 : CALL trexio_error(rc)
239 :
240 12 : ALLOCATE (coord(3, natoms))
241 8 : ALLOCATE (label(natoms))
242 12 : ALLOCATE (charge(natoms))
243 17 : DO iatom = 1, natoms
244 : ! store the coordinates
245 52 : 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 13 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
248 : ! store the element symbol
249 13 : label(iatom) = element_symbol
250 : ! get and store the effective nuclear charge of this kind_type (ikind)
251 13 : CALL get_qs_kind(kind_set(ikind), zeff=zeff)
252 17 : charge(iatom) = zeff
253 : END DO
254 :
255 4 : rc = trexio_write_nucleus_coord(f, coord)
256 4 : CALL trexio_error(rc)
257 4 : DEALLOCATE (coord)
258 :
259 4 : rc = trexio_write_nucleus_charge(f, charge)
260 4 : CALL trexio_error(rc)
261 4 : DEALLOCATE (charge)
262 :
263 4 : rc = trexio_write_nucleus_label(f, label, 3)
264 4 : CALL trexio_error(rc)
265 4 : DEALLOCATE (label)
266 :
267 : ! nuclear repulsion energy well-defined for molecules only
268 16 : IF (SUM(cell%perd) == 0) THEN
269 2 : CALL nuclear_repulsion_energy(particle_set, kind_set, e_nn)
270 2 : rc = trexio_write_nucleus_repulsion(f, e_nn)
271 2 : CALL trexio_error(rc)
272 : END IF
273 :
274 : !========================================================================================!
275 : ! Cell group
276 : !========================================================================================!
277 4 : rc = trexio_write_cell_a(f, cell%hmat(:, 1))
278 4 : CALL trexio_error(rc)
279 :
280 4 : rc = trexio_write_cell_b(f, cell%hmat(:, 2))
281 4 : CALL trexio_error(rc)
282 :
283 4 : rc = trexio_write_cell_c(f, cell%hmat(:, 3))
284 4 : CALL trexio_error(rc)
285 :
286 4 : rc = trexio_write_cell_g_a(f, cell%h_inv(:, 1))
287 4 : CALL trexio_error(rc)
288 :
289 4 : rc = trexio_write_cell_g_b(f, cell%h_inv(:, 2))
290 4 : CALL trexio_error(rc)
291 :
292 4 : rc = trexio_write_cell_g_c(f, cell%h_inv(:, 3))
293 4 : CALL trexio_error(rc)
294 :
295 4 : rc = trexio_write_cell_two_pi(f, 0)
296 4 : CALL trexio_error(rc)
297 :
298 : !========================================================================================!
299 : ! PBC group
300 : !========================================================================================!
301 4 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints)
302 :
303 4 : periodic = 0
304 16 : IF (SUM(cell%perd) /= 0) periodic = 1
305 4 : rc = trexio_write_pbc_periodic(f, periodic)
306 4 : CALL trexio_error(rc)
307 :
308 4 : IF (do_kpoints) THEN
309 1 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp)
310 :
311 1 : rc = trexio_write_pbc_k_point_num(f, nkp)
312 1 : CALL trexio_error(rc)
313 :
314 1 : rc = trexio_write_pbc_k_point(f, xkp)
315 1 : CALL trexio_error(rc)
316 :
317 1 : rc = trexio_write_pbc_k_point_weight(f, wkp)
318 1 : CALL trexio_error(rc)
319 : END IF
320 :
321 : !========================================================================================!
322 : ! Electron group
323 : !========================================================================================!
324 4 : CALL get_qs_env(qs_env, dft_control=dft_control, nelectron_total=nel_tot)
325 :
326 4 : rc = trexio_write_electron_num(f, nel_tot)
327 4 : CALL trexio_error(rc)
328 :
329 4 : nspins = dft_control%nspins
330 4 : 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 3 : nel_spin(1) = nel_tot/2
334 3 : 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 1 : CALL get_qs_env(qs_env, nelectron_spin=nel_spin)
339 : END IF
340 4 : rc = trexio_write_electron_up_num(f, nel_spin(1))
341 4 : CALL trexio_error(rc)
342 4 : rc = trexio_write_electron_dn_num(f, nel_spin(2))
343 4 : CALL trexio_error(rc)
344 :
345 : !========================================================================================!
346 : ! State group
347 : !========================================================================================!
348 4 : CALL get_qs_env(qs_env, energy=energy)
349 :
350 4 : rc = trexio_write_state_num(f, 1)
351 4 : CALL trexio_error(rc)
352 :
353 4 : rc = trexio_write_state_id(f, 1)
354 4 : CALL trexio_error(rc)
355 :
356 : ! rc = trexio_write_state_energy(f, energy%total)
357 4 : CALL trexio_error(rc)
358 :
359 : END IF ! ionode
360 :
361 : !========================================================================================!
362 : ! Basis group
363 : !========================================================================================!
364 8 : CALL get_qs_env(qs_env, qs_kind_set=kind_set, natom=natoms, particle_set=particle_set)
365 8 : CALL get_qs_kind_set(kind_set, nshell=shell_num, npgf_seg=prim_num)
366 :
367 8 : CALL section_vals_val_get(trexio_section, "CARTESIAN", l_val=save_cartesian)
368 :
369 8 : IF (ionode) THEN
370 4 : rc = trexio_write_basis_type(f, 'Gaussian', LEN_TRIM('Gaussian') + 1)
371 4 : CALL trexio_error(rc)
372 :
373 4 : rc = trexio_write_basis_shell_num(f, shell_num)
374 4 : CALL trexio_error(rc)
375 :
376 4 : rc = trexio_write_basis_prim_num(f, prim_num)
377 4 : CALL trexio_error(rc)
378 : END IF ! ionode
379 :
380 : ! one-to-one mapping between shells and ...
381 24 : ALLOCATE (nucleus_index(shell_num)) ! ...atomic indices
382 16 : ALLOCATE (shell_ang_mom(shell_num)) ! ...angular momenta
383 24 : ALLOCATE (shell_index(prim_num)) ! ...indices of primitive functions
384 24 : ALLOCATE (exponents(prim_num)) ! ...primitive exponents
385 16 : ALLOCATE (coefficients(prim_num)) ! ...contraction coefficients
386 16 : ALLOCATE (prim_factor(prim_num)) ! ...primitive normalization factors
387 :
388 : ! needed in AO group
389 8 : IF (.NOT. save_cartesian) THEN
390 8 : ALLOCATE (sgf_coefficients(prim_num)) ! ...contraction coefficients
391 : END IF
392 :
393 8 : ishell = 0 ! global shell index
394 8 : ipgf = 0 ! global primitives index
395 34 : DO iatom = 1, natoms
396 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
397 26 : 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 26 : 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 26 : l=l_shell_set)
408 :
409 148 : DO iset = 1, nset
410 252 : DO ishell_loc = 1, nshell(iset)
411 138 : ishell = ishell + 1
412 :
413 : ! nucleus_index array
414 138 : nucleus_index(ishell) = iatom
415 :
416 : ! shell_ang_mom array
417 138 : l = l_shell_set(ishell_loc, iset)
418 138 : shell_ang_mom(ishell) = l
419 :
420 : ! shell_index array
421 512 : shell_index(ipgf + 1:ipgf + npgf(iset)) = ishell
422 :
423 : ! exponents array
424 512 : 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 138 : expzet = 0.25_dp*REAL(2*l + 3, dp)
429 138 : prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
430 512 : DO i = 1, npgf(iset)
431 374 : gcca = gcc(i, ishell_loc, iset)
432 374 : zeta = zetas(i, iset)
433 374 : prim_cart_fac = prefac*zeta**expzet
434 :
435 : ! contraction coefficients array
436 374 : coefficients(ipgf + i) = gcca/prim_cart_fac
437 :
438 512 : IF (save_cartesian) THEN
439 : ! primitives normalization factors array
440 66 : prim_factor(ipgf + i) = prim_cart_fac
441 : ELSE
442 : ! for spherical harmonics we have a different factor
443 308 : prim_factor(ipgf + i) = sgf_norm(l, exponents(ipgf + i))
444 : ! we need these later in the AO group
445 308 : sgf_coefficients(ipgf + i) = coefficients(ipgf + i)*prim_factor(ipgf + i)
446 : END IF
447 :
448 : END DO
449 :
450 226 : ipgf = ipgf + npgf(iset)
451 : END DO
452 : END DO
453 : END DO
454 : ! just a failsafe check
455 8 : CPASSERT(ishell == shell_num)
456 8 : CPASSERT(ipgf == prim_num)
457 :
458 8 : IF (ionode) THEN
459 4 : rc = trexio_write_basis_nucleus_index(f, nucleus_index)
460 4 : CALL trexio_error(rc)
461 :
462 4 : rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
463 4 : CALL trexio_error(rc)
464 :
465 : ! Normalization factors are shoved in the AO group
466 12 : ALLOCATE (shell_factor(shell_num)) ! 1-to-1 map bw shells and normalization factors
467 73 : shell_factor(:) = 1.0_dp
468 4 : rc = trexio_write_basis_shell_factor(f, shell_factor)
469 4 : CALL trexio_error(rc)
470 4 : DEALLOCATE (shell_factor)
471 :
472 : ! This is always 0 for Gaussian basis sets
473 12 : ALLOCATE (r_power(shell_num)) ! 1-to-1 map bw shells radial function powers
474 73 : r_power(:) = 0
475 4 : rc = trexio_write_basis_r_power(f, r_power)
476 4 : CALL trexio_error(rc)
477 4 : DEALLOCATE (r_power)
478 :
479 4 : rc = trexio_write_basis_shell_index(f, shell_index)
480 4 : CALL trexio_error(rc)
481 :
482 4 : rc = trexio_write_basis_exponent(f, exponents)
483 4 : CALL trexio_error(rc)
484 :
485 4 : rc = trexio_write_basis_coefficient(f, coefficients)
486 4 : CALL trexio_error(rc)
487 :
488 : ! Normalization factors are shoved in the AO group
489 4 : rc = trexio_write_basis_prim_factor(f, prim_factor)
490 4 : CALL trexio_error(rc)
491 : END IF
492 :
493 8 : DEALLOCATE (nucleus_index)
494 8 : DEALLOCATE (shell_index)
495 8 : DEALLOCATE (exponents)
496 8 : DEALLOCATE (coefficients)
497 8 : 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 8 : IF (ionode) THEN
504 4 : CALL get_qs_kind_set(kind_set, sgp_potential_present=sgp_potential_present)
505 :
506 : ! figure out whether we actually have ECP potentials
507 4 : ecp_num = 0
508 4 : IF (sgp_potential_present) THEN
509 4 : DO iatom = 1, natoms
510 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
511 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
512 : ! get the the sgp_potential associated to this qs_kind
513 2 : CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential)
514 :
515 : ! get the info on the potential
516 6 : IF (ASSOCIATED(sgp_potential)) THEN
517 2 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
518 2 : IF (ecp_local) THEN
519 : ! get number of local terms
520 2 : CALL get_potential(potential=sgp_potential, nloc=nloc)
521 2 : ecp_num = ecp_num + nloc
522 : END IF
523 2 : IF (ecp_semi_local) THEN
524 : ! get number of semilocal terms
525 2 : CALL get_potential(potential=sgp_potential, npot=npot)
526 24 : 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 2 : IF (ecp_num > 0) THEN
534 6 : ALLOCATE (z_core(natoms))
535 4 : ALLOCATE (max_ang_mom_plus_1(natoms))
536 4 : max_ang_mom_plus_1(:) = 0
537 :
538 6 : ALLOCATE (ang_mom(ecp_num))
539 4 : ALLOCATE (nucleus_index(ecp_num))
540 6 : ALLOCATE (exponents(ecp_num))
541 4 : ALLOCATE (ecp_coefficients(ecp_num))
542 4 : ALLOCATE (powers(ecp_num))
543 :
544 4 : iecp = 0
545 4 : DO iatom = 1, natoms
546 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
547 2 : 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 2 : CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential, zeff=zeff)
550 :
551 : ! number of core electrons removed by the ECP
552 2 : z_core(iatom) = z - INT(zeff)
553 :
554 : ! get the info on the potential
555 4 : IF (ASSOCIATED(sgp_potential)) THEN
556 2 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
557 :
558 : ! deal with the local part
559 2 : IF (ecp_local) THEN
560 2 : CALL get_potential(potential=sgp_potential, nloc=nloc, sl_lmax=sl_lmax)
561 4 : ang_mom(iecp + 1:iecp + nloc) = sl_lmax + 1
562 4 : nucleus_index(iecp + 1:iecp + nloc) = iatom
563 4 : exponents(iecp + 1:iecp + nloc) = sgp_potential%bloc(1:nloc)
564 4 : ecp_coefficients(iecp + 1:iecp + nloc) = sgp_potential%aloc(1:nloc)
565 4 : 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 2 : IF (ecp_semi_local) THEN
571 2 : CALL get_potential(potential=sgp_potential, npot=npot, sl_lmax=sl_lmax)
572 2 : max_ang_mom_plus_1(iatom) = sl_lmax + 1
573 :
574 8 : DO sl_l = 0, sl_lmax
575 6 : nsemiloc = npot(sl_l)
576 16 : ang_mom(iecp + 1:iecp + nsemiloc) = sl_l
577 16 : nucleus_index(iecp + 1:iecp + nsemiloc) = iatom
578 16 : exponents(iecp + 1:iecp + nsemiloc) = sgp_potential%bpot(1:nsemiloc, sl_l)
579 16 : ecp_coefficients(iecp + 1:iecp + nsemiloc) = sgp_potential%apot(1:nsemiloc, sl_l)
580 16 : powers(iecp + 1:iecp + nsemiloc) = sgp_potential%nrpot(1:nsemiloc, sl_l) - 2
581 8 : iecp = iecp + nsemiloc
582 : END DO
583 : END IF
584 : END IF
585 : END DO
586 :
587 : ! fail-safe check
588 2 : CPASSERT(iecp == ecp_num)
589 :
590 2 : rc = trexio_write_ecp_num(f, ecp_num)
591 2 : CALL trexio_error(rc)
592 :
593 2 : rc = trexio_write_ecp_z_core(f, z_core)
594 2 : CALL trexio_error(rc)
595 2 : DEALLOCATE (z_core)
596 :
597 2 : rc = trexio_write_ecp_max_ang_mom_plus_1(f, max_ang_mom_plus_1)
598 2 : CALL trexio_error(rc)
599 2 : DEALLOCATE (max_ang_mom_plus_1)
600 :
601 2 : rc = trexio_write_ecp_ang_mom(f, ang_mom)
602 2 : CALL trexio_error(rc)
603 2 : DEALLOCATE (ang_mom)
604 :
605 2 : rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
606 2 : CALL trexio_error(rc)
607 2 : DEALLOCATE (nucleus_index)
608 :
609 2 : rc = trexio_write_ecp_exponent(f, exponents)
610 2 : CALL trexio_error(rc)
611 2 : DEALLOCATE (exponents)
612 :
613 2 : rc = trexio_write_ecp_coefficient(f, ecp_coefficients)
614 2 : CALL trexio_error(rc)
615 2 : DEALLOCATE (ecp_coefficients)
616 :
617 2 : rc = trexio_write_ecp_power(f, powers)
618 2 : CALL trexio_error(rc)
619 2 : 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 8 : CALL get_qs_env(qs_env, qs_kind_set=kind_set)
633 8 : CALL get_qs_kind_set(kind_set, ncgf=ncgf, nsgf=nsgf)
634 :
635 8 : IF (save_cartesian) THEN
636 4 : ao_num = ncgf
637 : ELSE
638 4 : ao_num = nsgf
639 : END IF
640 :
641 8 : IF (ionode) THEN
642 4 : IF (save_cartesian) THEN
643 2 : rc = trexio_write_ao_cartesian(f, 1)
644 : ELSE
645 2 : rc = trexio_write_ao_cartesian(f, 0)
646 : END IF
647 4 : CALL trexio_error(rc)
648 :
649 4 : rc = trexio_write_ao_num(f, ao_num)
650 4 : CALL trexio_error(rc)
651 : END IF
652 :
653 : ! one-to-one mapping between AOs and ...
654 24 : ALLOCATE (ao_shell(ao_num)) ! ..shells
655 24 : ALLOCATE (ao_normalization(ao_num)) ! ..normalization factors
656 :
657 8 : 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 8 : ALLOCATE (cp2k_to_trexio_ang_mom(ao_num))
662 4 : i = 0
663 108 : DO ishell = 1, shell_num
664 104 : l = shell_ang_mom(ishell)
665 388 : DO k = 1, 2*l + 1
666 284 : m = (-1)**k*FLOOR(REAL(k, KIND=dp)/2.0_dp)
667 388 : cp2k_to_trexio_ang_mom(i + k) = i + l + 1 + m
668 : END DO
669 108 : i = i + 2*l + 1
670 : END DO
671 4 : CPASSERT(i == ao_num)
672 : END IF
673 :
674 : ! we need to be consistent with the basis group on the shell indices
675 8 : ishell = 0 ! global shell index
676 8 : iao = 0 ! global AO index
677 8 : ipgf = 0 ! global primitives index
678 34 : DO iatom = 1, natoms
679 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
680 26 : 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 26 : 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 26 : l=l_shell_set)
692 :
693 26 : icgf_atom = 0
694 114 : DO iset = 1, nset
695 252 : DO ishell_loc = 1, nshell(iset)
696 : ! global shell index
697 138 : ishell = ishell + 1
698 : ! angular momentum l of this shell
699 138 : l = l_shell_set(ishell_loc, iset)
700 :
701 : ! number of AOs in this shell
702 138 : IF (save_cartesian) THEN
703 34 : nao_shell = nco(l)
704 : ELSE
705 104 : nao_shell = nso(l)
706 : END IF
707 :
708 : ! one-to-one mapping between AOs and shells
709 524 : ao_shell(iao + 1:iao + nao_shell) = ishell
710 :
711 : ! one-to-one mapping between AOs and normalization factors
712 138 : IF (save_cartesian) THEN
713 136 : 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 416 : ALLOCATE (Sloc(npgf(iset), npgf(iset)))
717 312 : ALLOCATE (Sgcc(npgf(iset)))
718 104 : 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 2968 : Sgcc(:) = MATMUL(Sloc, sgf_coefficients(ipgf + 1:ipgf + npgf(iset)))
722 412 : Nsgto = 1.0_dp/SQRT(DOT_PRODUCT(sgf_coefficients(ipgf + 1:ipgf + npgf(iset)), Sgcc))
723 :
724 104 : DEALLOCATE (Sloc)
725 104 : 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 388 : ao_normalization(iao + 1:iao + nao_shell) = Nsgto*SQRT((2*l + 1)/(4*pi))
731 : END IF
732 :
733 138 : ipgf = ipgf + npgf(iset)
734 138 : iao = iao + nao_shell
735 226 : icgf_atom = icgf_atom + nco(l)
736 : END DO
737 : END DO
738 : ! just a failsafe check
739 86 : CPASSERT(icgf_atom == ncgf_atom)
740 : END DO
741 :
742 8 : IF (ionode) THEN
743 4 : rc = trexio_write_ao_shell(f, ao_shell)
744 4 : CALL trexio_error(rc)
745 :
746 4 : rc = trexio_write_ao_normalization(f, ao_normalization)
747 4 : CALL trexio_error(rc)
748 : END IF
749 :
750 8 : DEALLOCATE (ao_shell)
751 8 : DEALLOCATE (ao_normalization)
752 8 : 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 8 : particle_set=particle_set, qs_kind_set=kind_set, blacs_env=blacs_env)
759 8 : nspins = dft_control%nspins
760 8 : CALL get_qs_kind_set(kind_set, nsgf=nsgf, ncgf=ncgf)
761 8 : nmo_spin = 0
762 :
763 : ! figure out that total number of MOs
764 8 : mo_num = 0
765 8 : IF (do_kpoints) THEN
766 2 : CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, use_real_wfn=use_real_wfn)
767 2 : CALL get_kpoint_env(kp_env(1)%kpoint_env, mos=mos_kp)
768 4 : DO ispin = 1, nspins
769 2 : CALL get_mo_set(mos_kp(1, ispin), nmo=nmo)
770 4 : nmo_spin(ispin) = nmo
771 : END DO
772 6 : 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 2 : nrow_global=nsgf, ncol_global=mo_num)
777 2 : CALL cp_fm_create(fm_mo_coeff, fm_struct)
778 2 : CALL cp_fm_set_all(fm_mo_coeff, 0.0_dp)
779 2 : IF (.NOT. use_real_wfn) THEN
780 2 : CALL cp_fm_create(fm_mo_coeff_im, fm_struct)
781 2 : CALL cp_fm_set_all(fm_mo_coeff_im, 0.0_dp)
782 : END IF
783 2 : CALL cp_fm_struct_release(fm_struct)
784 : ELSE
785 6 : CALL get_qs_env(qs_env, mos=mos)
786 14 : DO ispin = 1, nspins
787 8 : CALL get_mo_set(mos(ispin), nmo=nmo)
788 14 : nmo_spin(ispin) = nmo
789 : END DO
790 18 : mo_num = SUM(nmo_spin)
791 : END IF
792 :
793 : ! allocate all the arrays
794 32 : ALLOCATE (mo_coefficient(ao_num, mo_num))
795 33432 : mo_coefficient(:, :) = 0.0_dp
796 24 : ALLOCATE (mo_energy(mo_num))
797 436 : mo_energy(:) = 0.0_dp
798 16 : ALLOCATE (mo_occupation(mo_num))
799 436 : mo_occupation(:) = 0.0_dp
800 24 : ALLOCATE (mo_spin(mo_num))
801 436 : mo_spin(:) = 0
802 : ! extra arrays for kpoints
803 8 : IF (do_kpoints) THEN
804 6 : ALLOCATE (mo_coefficient_im(ao_num, mo_num))
805 26882 : mo_coefficient_im(:, :) = 0.0_dp
806 4 : ALLOCATE (mo_kpoint(mo_num))
807 258 : 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 8 : IF (do_kpoints) THEN
814 2 : CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, kp_range=kp_range)
815 :
816 4 : DO ispin = 1, nspins
817 20 : DO ikp = 1, nkp
818 16 : nmo = nmo_spin(ispin)
819 : ! global index to store the MOs
820 16 : imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
821 :
822 : ! do I have this kpoint on this rank?
823 18 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
824 16 : ikp_loc = ikp - kp_range(1) + 1
825 : ! get the mo set for this kpoint
826 16 : 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 16 : IF (mos_kp(1, ispin)%use_mo_coeff_b) THEN
830 0 : 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 16 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
835 :
836 : ! copy MO energies to local arrays
837 272 : mo_energy(imo + 1:imo + nmo) = mos_kp(1, ispin)%eigenvalues(1:nmo)
838 :
839 : ! copy MO occupations to local arrays
840 272 : 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 16 : IF (.NOT. use_real_wfn) THEN
844 16 : IF (mos_kp(2, ispin)%use_mo_coeff_b) THEN
845 0 : 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 16 : 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 0 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
854 0 : IF (.NOT. use_real_wfn) THEN
855 : CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff_im, &
856 0 : 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 8 : IF (do_kpoints) THEN
865 2 : CALL get_kpoint_info(kpoints, para_env_inter_kp=para_env_inter_kp)
866 2 : CALL para_env_inter_kp%sum(mo_energy)
867 2 : 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 18 : DO ispin = 1, nspins
872 : ! get number of MOs for this spin
873 10 : nmo = nmo_spin(ispin)
874 : ! allocate local temp array to transform the MOs of each kpoint/spin
875 40 : ALLOCATE (mos_sgf(nsgf, nmo))
876 9458 : mos_sgf(:, :) = 0.0_dp
877 :
878 10 : IF (do_kpoints) THEN
879 18 : DO ikp = 1, nkp
880 : ! global index to store the MOs
881 16 : imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
882 :
883 : ! store kpoint index
884 272 : mo_kpoint(imo + 1:imo + nmo) = ikp
885 : ! store the MO spins
886 272 : mo_spin(imo + 1:imo + nmo) = ispin - 1
887 :
888 : ! transform and store the MO coefficients
889 16 : CALL cp_fm_get_submatrix(fm_mo_coeff, mos_sgf, 1, imo + 1, nsgf, nmo)
890 16 : IF (save_cartesian) THEN
891 0 : 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 1680 : DO i = 1, nsgf
895 28304 : 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 18 : IF (.NOT. use_real_wfn) THEN
901 16 : CALL cp_fm_get_submatrix(fm_mo_coeff_im, mos_sgf, 1, imo + 1, nsgf, nmo)
902 16 : IF (save_cartesian) THEN
903 0 : 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 1680 : DO i = 1, nsgf
907 28304 : 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 8 : imo = (ispin - 1)*nmo_spin(1)
915 : ! store the MO energies
916 180 : mo_energy(imo + 1:imo + nmo) = mos(ispin)%eigenvalues
917 : ! store the MO occupations
918 180 : mo_occupation(imo + 1:imo + nmo) = mos(ispin)%occupation_numbers
919 : ! store the MO spins
920 180 : 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 8 : 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 8 : CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, mos_sgf)
927 :
928 8 : IF (save_cartesian) THEN
929 6 : 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 78 : DO i = 1, nsgf
933 2966 : 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 18 : DEALLOCATE (mos_sgf)
939 : END DO
940 :
941 8 : IF (ionode) THEN
942 4 : rc = trexio_write_mo_type(f, 'Canonical', LEN_TRIM('Canonical') + 1)
943 4 : CALL trexio_error(rc)
944 :
945 4 : rc = trexio_write_mo_num(f, mo_num)
946 4 : CALL trexio_error(rc)
947 :
948 4 : rc = trexio_write_mo_coefficient(f, mo_coefficient)
949 4 : CALL trexio_error(rc)
950 :
951 4 : rc = trexio_write_mo_energy(f, mo_energy)
952 4 : CALL trexio_error(rc)
953 :
954 4 : rc = trexio_write_mo_occupation(f, mo_occupation)
955 4 : CALL trexio_error(rc)
956 :
957 4 : rc = trexio_write_mo_spin(f, mo_spin)
958 4 : CALL trexio_error(rc)
959 :
960 4 : IF (do_kpoints) THEN
961 1 : rc = trexio_write_mo_coefficient_im(f, mo_coefficient_im)
962 1 : CALL trexio_error(rc)
963 :
964 1 : rc = trexio_write_mo_k_point(f, mo_kpoint)
965 1 : CALL trexio_error(rc)
966 : END IF
967 : END IF
968 :
969 8 : DEALLOCATE (mo_coefficient)
970 8 : DEALLOCATE (mo_energy)
971 8 : DEALLOCATE (mo_occupation)
972 8 : DEALLOCATE (mo_spin)
973 8 : IF (do_kpoints) THEN
974 2 : DEALLOCATE (mo_coefficient_im)
975 2 : DEALLOCATE (mo_kpoint)
976 2 : CALL cp_fm_release(fm_mo_coeff)
977 2 : 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 8 : IF (PRESENT(energy_derivative)) THEN
989 0 : filename_dE = TRIM(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
990 :
991 0 : ALLOCATE (dEdP(nsgf, nsgf))
992 0 : dEdP(:, :) = 0.0_dp
993 :
994 0 : DO ispin = 1, nspins
995 0 : CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
996 0 : 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 0 : row_offset=row_offset, col_offset=col_offset)
1001 :
1002 : ! Copy data from block to array
1003 0 : DO i = 1, row_size
1004 0 : DO j = 1, col_size
1005 0 : dEdP(row_offset + i - 1, col_offset + j - 1) = data_block(i, j)
1006 : END DO
1007 : END DO
1008 : END DO
1009 0 : CALL dbcsr_iterator_stop(iter)
1010 :
1011 : ! symmetrize the matrix if needed
1012 0 : SELECT CASE (dbcsr_get_matrix_type(energy_derivative(ispin)%matrix))
1013 : CASE (dbcsr_type_symmetric)
1014 0 : CALL symmetrize_matrix(dEdP, "upper_to_lower")
1015 : CASE (dbcsr_type_antisymmetric)
1016 0 : CALL symmetrize_matrix(dEdP, "anti_upper_to_lower")
1017 : CASE (dbcsr_type_no_symmetry)
1018 : CASE DEFAULT
1019 0 : CPABORT("Unknown matrix type for energy derivative")
1020 : END SELECT
1021 : END DO
1022 :
1023 : ! reduce the dEdP matrix to the master node
1024 0 : CALL para_env%sum(dEdP)
1025 :
1026 : ! print the dEdP matrix to a file
1027 0 : IF (ionode) THEN
1028 0 : WRITE (output_unit, "((T2,A,A))") 'TREXIO| Writing derivative file ', TRIM(filename_dE)
1029 :
1030 0 : unit_dE = 10
1031 : CALL open_file(file_name=filename_dE, &
1032 : file_action="WRITE", &
1033 : file_status="UNKNOWN", &
1034 0 : unit_number=unit_dE)
1035 0 : WRITE (unit_dE, '(I0, 1X, I0)') nsgf, nsgf
1036 0 : DO i = 1, nsgf
1037 : WRITE (unit_dE, '(*(1X, F15.8))') (dEdP(cp2k_to_trexio_ang_mom(i), &
1038 0 : cp2k_to_trexio_ang_mom(j)), j=1, nsgf)
1039 : END DO
1040 0 : CALL close_file(unit_number=unit_dE)
1041 : END IF
1042 :
1043 0 : DEALLOCATE (dEdP)
1044 : END IF
1045 :
1046 : ! Deallocate arrays used throughout the subroutine
1047 8 : IF (ALLOCATED(shell_ang_mom)) DEALLOCATE (shell_ang_mom)
1048 8 : IF (ALLOCATED(cp2k_to_trexio_ang_mom)) DEALLOCATE (cp2k_to_trexio_ang_mom)
1049 :
1050 : !========================================================================================!
1051 : ! Close the TREXIO file
1052 : !========================================================================================!
1053 8 : IF (ionode) THEN
1054 4 : rc = trexio_close(f)
1055 4 : CALL trexio_error(rc)
1056 : END IF
1057 :
1058 8 : 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 48 : 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 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:), ALLOCATABLE :: mo_spin, shell_ang_mom, trexio_to_cp2k_ang_mom
1103 :
1104 : REAL(KIND=dp) :: zeff, maxocc
1105 0 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: mo_energy, mo_occupation, charge
1106 0 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: mo_coefficient, mos_sgf, coord, dEdP, temp
1107 0 : 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 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1114 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1115 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1116 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1117 : TYPE(dbcsr_iterator_type) :: iter
1118 :
1119 0 : CALL timeset(routineN, handle)
1120 :
1121 0 : NULLIFY (logger, mo_coeff_ref, mo_coeff_target, para_env, dft_control, matrix_s, kind_set, mos, particle_set)
1122 :
1123 0 : logger => cp_get_default_logger()
1124 0 : output_unit = cp_logger_get_default_io_unit(logger)
1125 0 : myprint = logger%iter_info%print_level
1126 :
1127 0 : CPASSERT(ASSOCIATED(qs_env))
1128 :
1129 : ! get filename
1130 0 : IF (.NOT. PRESENT(trexio_filename)) THEN
1131 0 : filename = TRIM(logger%iter_info%project_name)//'-TREXIO.h5'
1132 0 : filename_dE = TRIM(logger%iter_info%project_name)//'-TREXIO.dEdP.dat'
1133 : ELSE
1134 0 : filename = TRIM(trexio_filename)//'.h5'
1135 0 : filename_dE = TRIM(trexio_filename)//'.dEdP.dat'
1136 : END IF
1137 :
1138 0 : CALL get_qs_env(qs_env, para_env=para_env)
1139 0 : ionode = para_env%is_source()
1140 :
1141 : ! Open the TREXIO file and check that we have the same molecule as in qs_env
1142 0 : IF (ionode) THEN
1143 0 : WRITE (output_unit, "((T2,A,A))") 'TREXIO| Opening file named ', TRIM(filename)
1144 0 : f = trexio_open(filename, 'r', TREXIO_HDF5, rc)
1145 0 : CALL trexio_error(rc)
1146 :
1147 0 : IF (myprint > medium_print_level) THEN
1148 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading molecule information...'
1149 : END IF
1150 0 : rc = trexio_read_nucleus_num(f, nucleus_num)
1151 0 : CALL trexio_error(rc)
1152 :
1153 0 : IF (myprint > medium_print_level) THEN
1154 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear coordinates...'
1155 : END IF
1156 0 : ALLOCATE (coord(3, nucleus_num))
1157 0 : rc = trexio_read_nucleus_coord(f, coord)
1158 0 : CALL trexio_error(rc)
1159 :
1160 0 : IF (myprint > medium_print_level) THEN
1161 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear labels...'
1162 : END IF
1163 0 : ALLOCATE (label(nucleus_num))
1164 0 : rc = trexio_read_nucleus_label(f, label, 3)
1165 0 : CALL trexio_error(rc)
1166 :
1167 0 : IF (myprint > medium_print_level) THEN
1168 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading nuclear charges...'
1169 : END IF
1170 0 : ALLOCATE (charge(nucleus_num))
1171 0 : rc = trexio_read_nucleus_charge(f, charge)
1172 0 : CALL trexio_error(rc)
1173 :
1174 : ! get the same info from qs_env
1175 0 : 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 0 : CPASSERT(nucleus_num == natoms)
1179 :
1180 0 : DO iatom = 1, natoms
1181 : ! compare the coordinates within a certain tolerance
1182 0 : DO i = 1, 3
1183 0 : 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 0 : 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 0 : CPASSERT(TRIM(element_symbol) == TRIM(label(iatom)))
1190 :
1191 : ! get the effective nuclear charge for this kind
1192 0 : CALL get_qs_kind(kind_set(ikind), zeff=zeff)
1193 : ! check that the nuclear charge is also the same
1194 0 : CPASSERT(charge(iatom) == zeff)
1195 : END DO
1196 :
1197 0 : 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 0 : DEALLOCATE (coord)
1200 0 : DEALLOCATE (label)
1201 0 : DEALLOCATE (charge)
1202 :
1203 : ! get info from trexio to map cp2k and trexio AOs
1204 0 : rc = trexio_read_ao_cartesian(f, save_cartesian)
1205 0 : CALL trexio_error(rc)
1206 :
1207 0 : rc = trexio_read_ao_num(f, ao_num)
1208 0 : CALL trexio_error(rc)
1209 :
1210 0 : rc = trexio_read_basis_shell_num(f, shell_num)
1211 0 : CALL trexio_error(rc)
1212 : END IF
1213 :
1214 0 : CALL para_env%bcast(save_cartesian, para_env%source)
1215 0 : CALL para_env%bcast(ao_num, para_env%source)
1216 0 : CALL para_env%bcast(shell_num, para_env%source)
1217 :
1218 0 : IF (save_cartesian == 1) THEN
1219 0 : 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 0 : CALL get_qs_env(qs_env, qs_kind_set=kind_set)
1224 0 : CALL get_qs_kind_set(kind_set, nsgf=nsgf, nshell=nshell)
1225 0 : CPASSERT(ao_num == nsgf)
1226 0 : CPASSERT(shell_num == nshell)
1227 :
1228 0 : ALLOCATE (shell_ang_mom(shell_num))
1229 0 : shell_ang_mom(:) = 0
1230 :
1231 0 : IF (ionode) THEN
1232 0 : IF (myprint > medium_print_level) THEN
1233 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading shell angular momenta...'
1234 : END IF
1235 0 : rc = trexio_read_basis_shell_ang_mom(f, shell_ang_mom)
1236 0 : CALL trexio_error(rc)
1237 : END IF
1238 :
1239 0 : 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 0 : ALLOCATE (trexio_to_cp2k_ang_mom(nsgf))
1245 0 : i = 0
1246 0 : DO ishell = 1, shell_num
1247 0 : l = shell_ang_mom(ishell)
1248 0 : DO k = 1, 2*l + 1
1249 0 : m = (-1)**k*FLOOR(REAL(k, KIND=dp)/2.0_dp)
1250 0 : trexio_to_cp2k_ang_mom(i + l + 1 + m) = i + k
1251 : END DO
1252 0 : i = i + 2*l + 1
1253 : END DO
1254 0 : CPASSERT(i == nsgf)
1255 :
1256 : ! check whether we want to read MOs
1257 0 : IF (PRESENT(mo_set_trexio)) THEN
1258 0 : IF (output_unit > 1) THEN
1259 0 : 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 0 : IF (ionode) THEN
1265 0 : rc = trexio_read_mo_num(f, mo_num)
1266 0 : CALL trexio_error(rc)
1267 :
1268 0 : rc = trexio_read_electron_up_num(f, electron_num(1))
1269 0 : CALL trexio_error(rc)
1270 :
1271 0 : rc = trexio_read_electron_dn_num(f, electron_num(2))
1272 0 : CALL trexio_error(rc)
1273 : END IF
1274 :
1275 : ! broadcast information to all processors and allocate arrays
1276 0 : CALL para_env%bcast(mo_num, para_env%source)
1277 0 : CALL para_env%bcast(electron_num, para_env%source)
1278 :
1279 : ! check that the number of MOs is the same
1280 0 : CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1281 0 : nspins = dft_control%nspins
1282 0 : nmo_spin(:) = 0
1283 0 : DO ispin = 1, nspins
1284 0 : CALL get_mo_set(mos(ispin), nmo=nmo)
1285 0 : nmo_spin(ispin) = nmo
1286 : END DO
1287 0 : CPASSERT(mo_num == SUM(nmo_spin))
1288 :
1289 0 : ALLOCATE (mo_coefficient(ao_num, mo_num))
1290 0 : ALLOCATE (mo_energy(mo_num))
1291 0 : ALLOCATE (mo_occupation(mo_num))
1292 0 : ALLOCATE (mo_spin(mo_num))
1293 :
1294 0 : mo_coefficient(:, :) = 0.0_dp
1295 0 : mo_energy(:) = 0.0_dp
1296 0 : mo_occupation(:) = 0.0_dp
1297 0 : mo_spin(:) = 0
1298 :
1299 : ! read the MOs info
1300 0 : IF (ionode) THEN
1301 0 : IF (myprint > medium_print_level) THEN
1302 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO coefficients...'
1303 : END IF
1304 0 : rc = trexio_read_mo_coefficient(f, mo_coefficient)
1305 0 : CALL trexio_error(rc)
1306 :
1307 0 : IF (myprint > medium_print_level) THEN
1308 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO energies...'
1309 : END IF
1310 0 : rc = trexio_read_mo_energy(f, mo_energy)
1311 0 : CALL trexio_error(rc)
1312 :
1313 0 : IF (myprint > medium_print_level) THEN
1314 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO occupations...'
1315 : END IF
1316 0 : rc = trexio_read_mo_occupation(f, mo_occupation)
1317 0 : CALL trexio_error(rc)
1318 :
1319 0 : IF (myprint > medium_print_level) THEN
1320 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading MO spins...'
1321 : END IF
1322 0 : rc = trexio_read_mo_spin(f, mo_spin)
1323 0 : CALL trexio_error(rc)
1324 : END IF
1325 :
1326 : ! broadcast the data to all processors
1327 0 : CALL para_env%bcast(mo_coefficient, para_env%source)
1328 0 : CALL para_env%bcast(mo_energy, para_env%source)
1329 0 : CALL para_env%bcast(mo_occupation, para_env%source)
1330 0 : 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 0 : DO ispin = 1, nspins
1335 : ! global MOs index
1336 0 : imo = (ispin - 1)*nmo_spin(1)
1337 : ! get number of MOs for this spin
1338 0 : nmo = nmo_spin(ispin)
1339 : ! allocate local temp array to read MOs
1340 0 : ALLOCATE (mos_sgf(nsgf, nmo))
1341 0 : mos_sgf(:, :) = 0.0_dp
1342 :
1343 : ! we need to reorder the MOs according to CP2K convention
1344 0 : DO i = 1, nsgf
1345 0 : mos_sgf(i, :) = mo_coefficient(trexio_to_cp2k_ang_mom(i), imo + 1:imo + nmo)
1346 : END DO
1347 :
1348 0 : IF (nspins == 1) THEN
1349 0 : maxocc = 2.0_dp
1350 0 : nelectron = electron_num(1) + electron_num(2)
1351 : ELSE
1352 0 : maxocc = 1.0_dp
1353 0 : nelectron = electron_num(ispin)
1354 : END IF
1355 : ! the right number of active electrons per spin channel is initialized further down
1356 0 : CALL allocate_mo_set(mo_set_trexio(ispin), nsgf, nmo, nelectron, 0.0_dp, maxocc, 0.0_dp)
1357 :
1358 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff_ref)
1359 0 : CALL init_mo_set(mo_set_trexio(ispin), fm_ref=mo_coeff_ref, name="TREXIO MOs")
1360 :
1361 0 : CALL get_mo_set(mo_set_trexio(ispin), mo_coeff=mo_coeff_target)
1362 0 : DO j = 1, nmo
1363 : ! make sure I copy the right spin channel
1364 0 : CPASSERT(mo_spin(j) == ispin - 1)
1365 0 : mo_set_trexio(ispin)%eigenvalues(j) = mo_energy(imo + j)
1366 0 : mo_set_trexio(ispin)%occupation_numbers(j) = mo_occupation(imo + j)
1367 0 : DO i = 1, nsgf
1368 0 : CALL cp_fm_set_element(mo_coeff_target, i, j, mos_sgf(i, j))
1369 : END DO
1370 : END DO
1371 :
1372 0 : DEALLOCATE (mos_sgf)
1373 : END DO
1374 :
1375 0 : DEALLOCATE (mo_coefficient)
1376 0 : DEALLOCATE (mo_energy)
1377 0 : DEALLOCATE (mo_occupation)
1378 0 : DEALLOCATE (mo_spin)
1379 :
1380 : END IF ! if MOs should be read
1381 :
1382 : ! check whether we want to read derivatives
1383 0 : IF (PRESENT(energy_derivative)) THEN
1384 0 : IF (output_unit > 1) THEN
1385 0 : 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 0 : ALLOCATE (temp(nsgf, nsgf))
1393 0 : temp(:, :) = 0.0_dp
1394 :
1395 : ! check if file exists and open it
1396 0 : IF (ionode) THEN
1397 0 : IF (file_exists(filename_dE)) THEN
1398 0 : CALL open_file(file_name=filename_dE, file_status="OLD", unit_number=unit_dE)
1399 : ELSE
1400 0 : CPABORT("Energy derivatives file "//TRIM(filename_dE)//" not found")
1401 : END IF
1402 :
1403 : ! read the header and check everything is fine
1404 0 : IF (myprint > medium_print_level) THEN
1405 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading header information...'
1406 : END IF
1407 0 : READ (unit_dE, *) nrows, ncols
1408 0 : IF (myprint > medium_print_level) THEN
1409 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Check size of dEdP matrix...'
1410 : END IF
1411 0 : CPASSERT(nrows == nsgf)
1412 0 : CPASSERT(ncols == nsgf)
1413 :
1414 : ! read the data
1415 0 : IF (myprint > medium_print_level) THEN
1416 0 : WRITE (output_unit, "((T2,A))") 'TREXIO| Reading dEdP matrix...'
1417 : END IF
1418 : ! Read the data matrix
1419 0 : DO i = 1, nrows
1420 0 : READ (unit_dE, *) (temp(i, j), j=1, ncols)
1421 : END DO
1422 :
1423 0 : CALL close_file(unit_number=unit_dE)
1424 : END IF
1425 :
1426 : ! send data to all processes
1427 0 : CALL para_env%bcast(temp, para_env%source)
1428 :
1429 : ! Reshuffle
1430 0 : ALLOCATE (dEdP(nsgf, nsgf))
1431 0 : dEdP(:, :) = 0.0_dp
1432 :
1433 : ! Reorder rows and columns according to trexio_to_cp2k_ang_mom mapping
1434 0 : DO j = 1, nsgf
1435 0 : DO i = 1, nsgf
1436 : ! either this
1437 0 : 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 0 : DEALLOCATE (temp)
1444 :
1445 0 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1446 0 : DO ispin = 1, nspins
1447 0 : 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 0 : name='Energy Derivative', keep_sparsity=.FALSE.)
1452 0 : CALL dbcsr_set(energy_derivative(ispin)%matrix, 0.0_dp)
1453 :
1454 0 : CALL dbcsr_iterator_start(iter, energy_derivative(ispin)%matrix)
1455 0 : 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 0 : row_offset=row_offset, col_offset=col_offset)
1459 :
1460 : ! Copy data from array to block
1461 0 : DO i = 1, row_size
1462 0 : DO j = 1, col_size
1463 0 : data_block(i, j) = dEdP(row_offset + i - 1, col_offset + j - 1)
1464 : END DO
1465 : END DO
1466 : END DO
1467 0 : CALL dbcsr_iterator_stop(iter)
1468 : END DO
1469 :
1470 0 : DEALLOCATE (dEdP)
1471 : END IF ! finished reading energy derivatives
1472 :
1473 : ! Clean up
1474 0 : IF (ALLOCATED(shell_ang_mom)) DEALLOCATE (shell_ang_mom)
1475 0 : IF (ALLOCATED(trexio_to_cp2k_ang_mom)) DEALLOCATE (trexio_to_cp2k_ang_mom)
1476 :
1477 : ! Close the TREXIO file
1478 0 : IF (ionode) THEN
1479 0 : WRITE (output_unit, "((T2,A,A))") 'TREXIO| Closing file named ', TRIM(filename)
1480 0 : rc = trexio_close(f)
1481 0 : CALL trexio_error(rc)
1482 : END IF
1483 :
1484 0 : 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 0 : END SUBROUTINE read_trexio
1496 :
1497 : #ifdef __TREXIO
1498 : ! **************************************************************************************************
1499 : !> \brief Handles TREXIO errors
1500 : !> \param rc the TREXIO return code
1501 : ! **************************************************************************************************
1502 195 : SUBROUTINE trexio_error(rc)
1503 : INTEGER(trexio_exit_code), INTENT(IN) :: rc
1504 :
1505 : CHARACTER(LEN=128) :: err_msg
1506 :
1507 195 : IF (rc /= TREXIO_SUCCESS) THEN
1508 0 : CALL trexio_string_of_error(rc, err_msg)
1509 0 : CPABORT('TREXIO Error: '//TRIM(err_msg))
1510 : END IF
1511 :
1512 195 : 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 2 : 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 2 : natoms = SIZE(particle_set)
1531 2 : e_nn = 0.0_dp
1532 4 : DO i = 1, natoms
1533 2 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
1534 2 : CALL get_qs_kind(kind_set(ikind), zeff=zeff_i)
1535 4 : DO j = i + 1, natoms
1536 0 : r_ij = NORM2(particle_set(i)%r - particle_set(j)%r)
1537 :
1538 0 : CALL get_atomic_kind(particle_set(j)%atomic_kind, kind_number=jkind)
1539 0 : CALL get_qs_kind(kind_set(jkind), zeff=zeff_j)
1540 :
1541 2 : e_nn = e_nn + zeff_i*zeff_j/r_ij
1542 : END DO
1543 : END DO
1544 :
1545 2 : 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 308 : 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 308 : IF (l >= 0) THEN
1559 308 : norm = SQRT(2**(2*l + 3)*fac(l + 1)*(2*expnt)**(l + 1.5)/(fac(2*l + 2)*SQRT(pi)))
1560 : ELSE
1561 0 : CPABORT("The angular momentum should be >= 0!")
1562 : END IF
1563 :
1564 308 : 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 6 : 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 6 : INTEGER, DIMENSION(:), POINTER :: nshell
1584 6 : INTEGER, DIMENSION(:, :), POINTER :: l
1585 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1586 :
1587 6 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1588 :
1589 3586 : mos_cgf(:, :) = 0.0_dp
1590 6 : nmo = SIZE(mos_sgf, 2)
1591 :
1592 : ! Transform spherical MOs to Cartesian MOs
1593 6 : icgf = 1
1594 6 : isgf = 1
1595 20 : DO iatom = 1, SIZE(particle_set)
1596 14 : NULLIFY (orb_basis_set)
1597 14 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1598 14 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1599 :
1600 34 : 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 14 : l=l)
1605 54 : DO iset = 1, nset
1606 98 : DO ishell = 1, nshell(iset)
1607 44 : 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 44 : mos_cgf(icgf, 1), ncgf)
1612 44 : icgf = icgf + nco(lshell)
1613 84 : isgf = isgf + nso(lshell)
1614 : END DO
1615 : END DO
1616 : ELSE
1617 : ! assume atom without basis set
1618 0 : CPABORT("Unknown basis set type")
1619 : END IF
1620 : END DO ! iatom
1621 :
1622 6 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: nshell
1642 0 : INTEGER, DIMENSION(:, :), POINTER :: l
1643 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1644 :
1645 0 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1646 :
1647 0 : mos_sgf(:, :) = 0.0_dp
1648 0 : nmo = SIZE(mos_cgf, 2)
1649 :
1650 : ! Transform Cartesian MOs to spherical MOs
1651 0 : icgf = 1
1652 0 : isgf = 1
1653 0 : DO iatom = 1, SIZE(particle_set)
1654 0 : NULLIFY (orb_basis_set)
1655 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1656 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1657 :
1658 0 : 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 0 : l=l)
1663 0 : DO iset = 1, nset
1664 0 : DO ishell = 1, nshell(iset)
1665 0 : 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 0 : mos_sgf(isgf, 1), nsgf)
1670 0 : icgf = icgf + nco(lshell)
1671 0 : isgf = isgf + nso(lshell)
1672 : END DO
1673 : END DO
1674 : ELSE
1675 : ! assume atom without basis set
1676 0 : CPABORT("Unknown basis set type")
1677 : END IF
1678 : END DO ! iatom
1679 :
1680 0 : END SUBROUTINE cartesian_to_spherical_mo
1681 : #endif
1682 :
1683 : END MODULE trexio_utils
|
