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