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