Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 : USE cp_fm_struct, ONLY: cp_fm_struct_create, &
27 : cp_fm_struct_release, &
28 : cp_fm_struct_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger, &
30 : cp_logger_get_default_io_unit, &
31 : cp_logger_type
32 : USE external_potential_types, ONLY: sgp_potential_type, get_potential
33 : USE input_section_types, ONLY: section_vals_type, section_vals_get, &
34 : section_vals_val_get
35 : USE kinds, ONLY: default_path_length, dp
36 : USE kpoint_types, ONLY: kpoint_type, kpoint_env_p_type, &
37 : get_kpoint_info, get_kpoint_env
38 : USE mathconstants, ONLY: fourpi, pi
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE orbital_pointers, ONLY: nco, nso
41 : USE orbital_transformation_matrices, ONLY: orbtramat
42 : USE particle_types, ONLY: particle_type
43 : USE qs_energy_types, ONLY: qs_energy_type
44 : USE qs_environment_types, ONLY: get_qs_env, &
45 : qs_environment_type
46 : USE qs_kind_types, ONLY: get_qs_kind, get_qs_kind_set, &
47 : qs_kind_type
48 : USE qs_mo_types, ONLY: mo_set_type, get_mo_set
49 : #ifdef __TREXIO
50 : USE trexio, ONLY: trexio_open, trexio_close, &
51 : TREXIO_HDF5, TREXIO_SUCCESS, &
52 : trexio_string_of_error, trexio_t, trexio_exit_code, &
53 : trexio_write_metadata_code, trexio_write_metadata_code_num, &
54 : trexio_write_nucleus_coord, trexio_write_nucleus_num, &
55 : trexio_write_nucleus_charge, trexio_write_nucleus_label, &
56 : trexio_write_nucleus_repulsion, &
57 : trexio_write_cell_a, trexio_write_cell_b, trexio_write_cell_c, &
58 : trexio_write_cell_g_a, trexio_write_cell_g_b, &
59 : trexio_write_cell_g_c, trexio_write_cell_two_pi, &
60 : trexio_write_pbc_periodic, trexio_write_pbc_k_point_num, &
61 : trexio_write_pbc_k_point, trexio_write_pbc_k_point_weight, &
62 : trexio_write_electron_num, trexio_write_electron_up_num, &
63 : trexio_write_electron_dn_num, &
64 : trexio_write_state_num, trexio_write_state_id, &
65 : trexio_write_state_energy, &
66 : trexio_write_basis_type, trexio_write_basis_prim_num, &
67 : trexio_write_basis_shell_num, trexio_write_basis_nucleus_index, &
68 : trexio_write_basis_shell_ang_mom, trexio_write_basis_shell_factor, &
69 : trexio_write_basis_r_power, trexio_write_basis_shell_index, &
70 : trexio_write_basis_exponent, trexio_write_basis_coefficient, &
71 : trexio_write_basis_prim_factor, &
72 : trexio_write_ecp_z_core, trexio_write_ecp_max_ang_mom_plus_1, &
73 : trexio_write_ecp_num, trexio_write_ecp_ang_mom, &
74 : trexio_write_ecp_nucleus_index, trexio_write_ecp_exponent, &
75 : trexio_write_ecp_coefficient, trexio_write_ecp_power, &
76 : trexio_write_ao_cartesian, trexio_write_ao_num, &
77 : trexio_write_ao_shell, trexio_write_ao_normalization, &
78 : trexio_write_mo_num, trexio_write_mo_energy, &
79 : trexio_write_mo_occupation, trexio_write_mo_spin, &
80 : trexio_write_mo_class, trexio_write_mo_coefficient, &
81 : trexio_write_mo_coefficient_im, trexio_write_mo_k_point, &
82 : trexio_write_mo_type
83 : #endif
84 : #include "./base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'trexio_utils'
91 :
92 : PUBLIC :: write_trexio
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief Write a trexio file
98 : !> \param qs_env the qs environment with all the info of the computation
99 : !> \param trexio_section the section with the trexio info
100 : ! **************************************************************************************************
101 8 : SUBROUTINE write_trexio(qs_env, trexio_section)
102 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
103 : TYPE(section_vals_type), INTENT(IN), POINTER :: trexio_section
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_trexio'
106 :
107 : INTEGER :: handle
108 :
109 : #ifdef __TREXIO
110 : INTEGER :: output_unit, unit_trexio
111 : CHARACTER(len=default_path_length) :: filename
112 : INTEGER(trexio_t) :: f ! The TREXIO file handle
113 : INTEGER(trexio_exit_code) :: rc ! TREXIO return code
114 : LOGICAL :: explicit, do_kpoints, ecp_semi_local, &
115 : ecp_local, sgp_potential_present, ionode, &
116 : use_real_wfn, save_cartesian
117 : REAL(KIND=dp) :: e_nn, zeff, expzet, prefac, zeta, gcca
118 : TYPE(cell_type), POINTER :: cell => Null()
119 : TYPE(cp_logger_type), POINTER :: logger => Null()
120 : TYPE(dft_control_type), POINTER :: dft_control => Null()
121 : TYPE(gto_basis_set_type), POINTER :: basis_set
122 : TYPE(kpoint_type), POINTER :: kpoints => Null()
123 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => Null()
124 : TYPE(qs_energy_type), POINTER :: energy => Null()
125 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set => Null()
126 : TYPE(sgp_potential_type), POINTER :: sgp_potential => Null()
127 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos => Null()
128 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp => Null()
129 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env => Null()
130 : TYPE(mp_para_env_type), POINTER :: para_env => Null(), para_env_inter_kp => Null()
131 : TYPE(cp_blacs_env_type), POINTER :: blacs_env => Null()
132 : TYPE(cp_fm_struct_type), POINTER :: fm_struct => Null()
133 : TYPE(cp_fm_type) :: fm_mo_coeff, fm_dummy, fm_mo_coeff_im
134 :
135 : CHARACTER(LEN=2) :: element_symbol
136 8 : CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: label
137 : INTEGER :: iatom, natoms, periodic, nkp, nel_tot, &
138 : nspins, ikind, ishell_loc, ishell, &
139 : shell_num, prim_num, nset, iset, ipgf, z, &
140 : sl_lmax, ecp_num, nloc, nsemiloc, sl_l, iecp, &
141 : igf, icgf, ncgf, ngf_shell, lshell, ao_num, nmo, &
142 : mo_num, ispin, ikp, imo, ikp_loc, nsgf, &
143 : i, k, l, m
144 : INTEGER, DIMENSION(2) :: nel_spin, kp_range, nmo_spin
145 : INTEGER, DIMENSION(3) :: nkp_grid
146 : INTEGER, DIMENSION(0:10) :: npot
147 8 : INTEGER, DIMENSION(:), ALLOCATABLE :: nucleus_index, shell_ang_mom, r_power, &
148 8 : shell_index, z_core, max_ang_mom_plus_1, &
149 8 : ang_mom, powers, ao_shell, mo_spin, mo_kpoint
150 : INTEGER, DIMENSION(:), POINTER :: nshell => Null(), npgf => Null()
151 : INTEGER, DIMENSION(:, :), POINTER :: l_shell_set => Null()
152 8 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: charge, shell_factor, exponents, coefficients, &
153 8 : prim_factor, ao_normalization, mo_energy, &
154 8 : mo_occupation
155 : REAL(KIND=dp), DIMENSION(:), POINTER :: wkp => Null(), norm_cgf => Null()
156 8 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: coord, mo_coefficient, mo_coefficient_im, &
157 8 : mos_sgf, diag_nsgf, diag_ncgf, temp
158 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zetas => Null()
159 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc => Null()
160 : #endif
161 :
162 8 : CALL timeset(routineN, handle)
163 :
164 : #ifdef __TREXIO
165 8 : logger => cp_get_default_logger()
166 8 : output_unit = cp_logger_get_default_io_unit(logger)
167 :
168 8 : CPASSERT(ASSOCIATED(qs_env))
169 :
170 : ! get filename
171 8 : CALL section_vals_val_get(trexio_section, "FILENAME", c_val=filename, explicit=explicit)
172 8 : IF (.NOT. explicit) THEN
173 6 : filename = TRIM(logger%iter_info%project_name)//'-TREXIO.h5'
174 : ELSE
175 2 : filename = TRIM(filename)//'.h5'
176 : END IF
177 :
178 8 : CALL get_qs_env(qs_env, para_env=para_env)
179 8 : ionode = para_env%is_source()
180 :
181 : ! inquire whether a file with the same name already exists, if yes, delete it
182 8 : IF (ionode) THEN
183 4 : IF (file_exists(filename)) THEN
184 0 : CALL open_file(filename, unit_number=unit_trexio)
185 0 : CALL close_file(unit_number=unit_trexio, file_status="DELETE")
186 : END IF
187 :
188 : !========================================================================================!
189 : ! Open the TREXIO file
190 : !========================================================================================!
191 4 : f = trexio_open(filename, 'w', TREXIO_HDF5, rc)
192 4 : CALL trexio_error(rc)
193 :
194 : !========================================================================================!
195 : ! Metadata group
196 : !========================================================================================!
197 4 : rc = trexio_write_metadata_code_num(f, 1)
198 4 : CALL trexio_error(rc)
199 :
200 4 : rc = trexio_write_metadata_code(f, cp2k_version, LEN_TRIM(cp2k_version) + 1)
201 4 : CALL trexio_error(rc)
202 :
203 : !========================================================================================!
204 : ! Nucleus group
205 : !========================================================================================!
206 4 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=kind_set, natom=natoms)
207 :
208 4 : rc = trexio_write_nucleus_num(f, natoms)
209 4 : CALL trexio_error(rc)
210 :
211 12 : ALLOCATE (coord(3, natoms))
212 8 : ALLOCATE (label(natoms))
213 12 : ALLOCATE (charge(natoms))
214 17 : DO iatom = 1, natoms
215 : ! store the coordinates
216 52 : coord(:, iatom) = particle_set(iatom)%r(1:3)
217 : ! figure out the element symbol and to which kind_set entry this atomic_kind corresponds to
218 13 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=element_symbol, kind_number=ikind)
219 : ! store the element symbol
220 13 : label(iatom) = element_symbol
221 : ! get and store the effective nuclear charge of this kind_type (ikind)
222 13 : CALL get_qs_kind(kind_set(ikind), zeff=zeff)
223 17 : charge(iatom) = zeff
224 : END DO
225 :
226 4 : rc = trexio_write_nucleus_coord(f, coord)
227 4 : CALL trexio_error(rc)
228 4 : DEALLOCATE (coord)
229 :
230 4 : rc = trexio_write_nucleus_charge(f, charge)
231 4 : CALL trexio_error(rc)
232 4 : DEALLOCATE (charge)
233 :
234 4 : rc = trexio_write_nucleus_label(f, label, 3)
235 4 : CALL trexio_error(rc)
236 4 : DEALLOCATE (label)
237 :
238 : ! nuclear repulsion energy well-defined for molecules only
239 16 : IF (SUM(cell%perd) == 0) THEN
240 2 : CALL nuclear_repulsion_energy(particle_set, kind_set, e_nn)
241 2 : rc = trexio_write_nucleus_repulsion(f, e_nn)
242 2 : CALL trexio_error(rc)
243 : END IF
244 :
245 : !========================================================================================!
246 : ! Cell group
247 : !========================================================================================!
248 4 : rc = trexio_write_cell_a(f, cell%hmat(:, 1))
249 4 : CALL trexio_error(rc)
250 :
251 4 : rc = trexio_write_cell_b(f, cell%hmat(:, 2))
252 4 : CALL trexio_error(rc)
253 :
254 4 : rc = trexio_write_cell_c(f, cell%hmat(:, 3))
255 4 : CALL trexio_error(rc)
256 :
257 4 : rc = trexio_write_cell_g_a(f, cell%h_inv(:, 1))
258 4 : CALL trexio_error(rc)
259 :
260 4 : rc = trexio_write_cell_g_b(f, cell%h_inv(:, 2))
261 4 : CALL trexio_error(rc)
262 :
263 4 : rc = trexio_write_cell_g_c(f, cell%h_inv(:, 3))
264 4 : CALL trexio_error(rc)
265 :
266 4 : rc = trexio_write_cell_two_pi(f, 0)
267 4 : CALL trexio_error(rc)
268 :
269 : !========================================================================================!
270 : ! PBC group
271 : !========================================================================================!
272 4 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints)
273 :
274 4 : periodic = 0
275 16 : IF (SUM(cell%perd) /= 0) periodic = 1
276 4 : rc = trexio_write_pbc_periodic(f, periodic)
277 4 : CALL trexio_error(rc)
278 :
279 4 : IF (do_kpoints) THEN
280 1 : CALL get_kpoint_info(kpoints, nkp=nkp, nkp_grid=nkp_grid, wkp=wkp)
281 :
282 1 : rc = trexio_write_pbc_k_point_num(f, nkp)
283 1 : CALL trexio_error(rc)
284 :
285 4 : rc = trexio_write_pbc_k_point(f, REAL(nkp_grid, KIND=dp))
286 1 : CALL trexio_error(rc)
287 :
288 1 : rc = trexio_write_pbc_k_point_weight(f, wkp)
289 1 : CALL trexio_error(rc)
290 : END IF
291 :
292 : !========================================================================================!
293 : ! Electron group
294 : !========================================================================================!
295 4 : CALL get_qs_env(qs_env, dft_control=dft_control, nelectron_total=nel_tot)
296 :
297 4 : rc = trexio_write_electron_num(f, nel_tot)
298 4 : CALL trexio_error(rc)
299 :
300 4 : nspins = dft_control%nspins
301 4 : IF (nspins == 1) THEN
302 : ! it is a spin-restricted calculation and we need to split the electrons manually,
303 : ! because in CP2K they are all otherwise weirdly stored in nelectron_spin(1)
304 3 : nel_spin(1) = nel_tot/2
305 3 : nel_spin(2) = nel_tot/2
306 : ELSE
307 : ! for UKS/ROKS, the two spin channels are populated correctly and according to
308 : ! the multiplicity
309 1 : CALL get_qs_env(qs_env, nelectron_spin=nel_spin)
310 : END IF
311 4 : rc = trexio_write_electron_up_num(f, nel_spin(1))
312 4 : CALL trexio_error(rc)
313 4 : rc = trexio_write_electron_dn_num(f, nel_spin(2))
314 4 : CALL trexio_error(rc)
315 :
316 : !========================================================================================!
317 : ! State group
318 : !========================================================================================!
319 4 : CALL get_qs_env(qs_env, energy=energy)
320 :
321 4 : rc = trexio_write_state_num(f, 1)
322 4 : CALL trexio_error(rc)
323 :
324 4 : rc = trexio_write_state_id(f, 1)
325 4 : CALL trexio_error(rc)
326 :
327 4 : rc = trexio_write_state_energy(f, energy%total)
328 4 : CALL trexio_error(rc)
329 :
330 : END IF ! ionode
331 :
332 : !========================================================================================!
333 : ! Basis group
334 : !========================================================================================!
335 8 : CALL get_qs_env(qs_env, qs_kind_set=kind_set, natom=natoms, particle_set=particle_set)
336 8 : CALL get_qs_kind_set(kind_set, nshell=shell_num, npgf_seg=prim_num)
337 :
338 8 : IF (ionode) THEN
339 4 : rc = trexio_write_basis_type(f, 'Gaussian', LEN_TRIM('Gaussian') + 1)
340 4 : CALL trexio_error(rc)
341 :
342 4 : rc = trexio_write_basis_shell_num(f, shell_num)
343 4 : CALL trexio_error(rc)
344 :
345 4 : rc = trexio_write_basis_prim_num(f, prim_num)
346 4 : CALL trexio_error(rc)
347 : END IF ! ionode
348 :
349 : ! one-to-one mapping between shells and ...
350 24 : ALLOCATE (nucleus_index(shell_num)) ! ...atomic indices
351 16 : ALLOCATE (shell_ang_mom(shell_num)) ! ...angular momenta
352 24 : ALLOCATE (shell_index(prim_num)) ! ...indices of primitive functions
353 24 : ALLOCATE (exponents(prim_num)) ! ...primitive exponents
354 16 : ALLOCATE (coefficients(prim_num)) ! ...contraction coefficients
355 16 : ALLOCATE (prim_factor(prim_num)) ! ...primitive normalization factors
356 :
357 8 : ishell = 0
358 8 : ipgf = 0
359 34 : DO iatom = 1, natoms
360 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
361 26 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
362 : ! get the primary (orbital) basis set associated to this qs_kind
363 26 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
364 : ! get the info from the basis set
365 : CALL get_gto_basis_set(basis_set, &
366 : nset=nset, &
367 : nshell=nshell, &
368 : npgf=npgf, &
369 : zet=zetas, &
370 : gcc=gcc, &
371 26 : l=l_shell_set)
372 :
373 148 : DO iset = 1, nset
374 252 : DO ishell_loc = 1, nshell(iset)
375 138 : ishell = ishell + 1
376 :
377 : ! nucleus_index array
378 138 : nucleus_index(ishell) = iatom
379 :
380 : ! shell_ang_mom array
381 138 : l = l_shell_set(ishell_loc, iset)
382 138 : shell_ang_mom(ishell) = l
383 :
384 : ! shell_index array
385 512 : shell_index(ipgf + 1:ipgf + npgf(iset)) = ishell
386 :
387 : ! exponents array
388 512 : exponents(ipgf + 1:ipgf + npgf(iset)) = zetas(1:npgf(iset), iset)
389 :
390 : ! compute on the fly the normalization factor as in normalise_gcc_orb
391 : ! and recover the original contraction coefficients to store them separately
392 138 : expzet = 0.25_dp*REAL(2*l + 3, dp)
393 138 : prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
394 512 : DO i = 1, npgf(iset)
395 374 : gcca = gcc(i, ishell_loc, iset)
396 374 : zeta = zetas(i, iset)
397 :
398 : ! primitives normalization factors array
399 374 : prim_factor(i + ipgf) = prefac*zeta**expzet
400 :
401 : ! contractio coefficients array
402 512 : coefficients(i + ipgf) = gcca/prim_factor(i + ipgf)
403 : END DO
404 :
405 226 : ipgf = ipgf + npgf(iset)
406 : END DO
407 : END DO
408 : END DO
409 : ! just a failsafe check
410 8 : CPASSERT(ishell == shell_num)
411 8 : CPASSERT(ipgf == prim_num)
412 :
413 8 : IF (ionode) THEN
414 4 : rc = trexio_write_basis_nucleus_index(f, nucleus_index)
415 4 : CALL trexio_error(rc)
416 :
417 4 : rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
418 4 : CALL trexio_error(rc)
419 :
420 : ! Normalization factors are shoved in the AO group
421 12 : ALLOCATE (shell_factor(shell_num)) ! 1-to-1 map bw shells and normalization factors
422 73 : shell_factor(:) = 1.0_dp
423 4 : rc = trexio_write_basis_shell_factor(f, shell_factor)
424 4 : CALL trexio_error(rc)
425 4 : DEALLOCATE (shell_factor)
426 :
427 : ! This is always 0 for Gaussian basis sets
428 12 : ALLOCATE (r_power(shell_num)) ! 1-to-1 map bw shells radial function powers
429 73 : r_power(:) = 0
430 4 : rc = trexio_write_basis_r_power(f, r_power)
431 4 : CALL trexio_error(rc)
432 4 : DEALLOCATE (r_power)
433 :
434 4 : rc = trexio_write_basis_shell_index(f, shell_index)
435 4 : CALL trexio_error(rc)
436 :
437 4 : rc = trexio_write_basis_exponent(f, exponents)
438 4 : CALL trexio_error(rc)
439 :
440 4 : rc = trexio_write_basis_coefficient(f, coefficients)
441 4 : CALL trexio_error(rc)
442 :
443 : ! Normalization factors are shoved in the AO group
444 4 : rc = trexio_write_basis_prim_factor(f, prim_factor)
445 4 : CALL trexio_error(rc)
446 : END IF
447 :
448 8 : DEALLOCATE (nucleus_index)
449 8 : DEALLOCATE (shell_index)
450 8 : DEALLOCATE (exponents)
451 8 : DEALLOCATE (coefficients)
452 8 : DEALLOCATE (prim_factor)
453 : ! shell_ang_mom is needed in the MO group, so will be deallocated there
454 :
455 : !========================================================================================!
456 : ! ECP group
457 : !========================================================================================!
458 8 : IF (ionode) THEN
459 4 : CALL get_qs_kind_set(kind_set, sgp_potential_present=sgp_potential_present)
460 :
461 : ! figure out whether we actually have ECP potentials
462 4 : ecp_num = 0
463 4 : IF (sgp_potential_present) THEN
464 4 : DO iatom = 1, natoms
465 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
466 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
467 : ! get the the sgp_potential associated to this qs_kind
468 2 : CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential)
469 :
470 : ! get the info on the potential
471 6 : IF (ASSOCIATED(sgp_potential)) THEN
472 2 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
473 2 : IF (ecp_local) THEN
474 : ! get number of local terms
475 2 : CALL get_potential(potential=sgp_potential, nloc=nloc)
476 2 : ecp_num = ecp_num + nloc
477 : END IF
478 2 : IF (ecp_semi_local) THEN
479 : ! get number of semilocal terms
480 2 : CALL get_potential(potential=sgp_potential, npot=npot)
481 24 : ecp_num = ecp_num + SUM(npot)
482 : END IF
483 : END IF
484 : END DO
485 : END IF
486 :
487 : ! if we have ECP potentials, populate the ECP group
488 2 : IF (ecp_num > 0) THEN
489 6 : ALLOCATE (z_core(natoms))
490 4 : ALLOCATE (max_ang_mom_plus_1(natoms))
491 4 : max_ang_mom_plus_1(:) = 0
492 :
493 6 : ALLOCATE (ang_mom(ecp_num))
494 4 : ALLOCATE (nucleus_index(ecp_num))
495 6 : ALLOCATE (exponents(ecp_num))
496 4 : ALLOCATE (coefficients(ecp_num))
497 4 : ALLOCATE (powers(ecp_num))
498 :
499 4 : iecp = 0
500 4 : DO iatom = 1, natoms
501 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
502 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind, z=z)
503 : ! get the the sgp_potential associated to this qs_kind
504 2 : CALL get_qs_kind(kind_set(ikind), sgp_potential=sgp_potential, zeff=zeff)
505 :
506 : ! number of core electrons removed by the ECP
507 2 : z_core(iatom) = z - INT(zeff)
508 :
509 : ! get the info on the potential
510 4 : IF (ASSOCIATED(sgp_potential)) THEN
511 2 : CALL get_potential(potential=sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
512 :
513 : ! deal with the local part
514 2 : IF (ecp_local) THEN
515 2 : CALL get_potential(potential=sgp_potential, nloc=nloc, sl_lmax=sl_lmax)
516 4 : ang_mom(iecp + 1:iecp + nloc) = sl_lmax + 1
517 4 : nucleus_index(iecp + 1:iecp + nloc) = iatom
518 4 : exponents(iecp + 1:iecp + nloc) = sgp_potential%bloc(1:nloc)
519 4 : coefficients(iecp + 1:iecp + nloc) = sgp_potential%aloc(1:nloc)
520 4 : powers(iecp + 1:iecp + nloc) = sgp_potential%nrloc(1:nloc) - 2
521 : iecp = iecp + nloc
522 : END IF
523 :
524 : ! deal with the semilocal part
525 2 : IF (ecp_semi_local) THEN
526 2 : CALL get_potential(potential=sgp_potential, npot=npot, sl_lmax=sl_lmax)
527 2 : max_ang_mom_plus_1(iatom) = sl_lmax + 1
528 :
529 8 : DO sl_l = 0, sl_lmax
530 6 : nsemiloc = npot(sl_l)
531 16 : ang_mom(iecp + 1:iecp + nsemiloc) = sl_l
532 16 : nucleus_index(iecp + 1:iecp + nsemiloc) = iatom
533 16 : exponents(iecp + 1:iecp + nsemiloc) = sgp_potential%bpot(1:nsemiloc, sl_l)
534 16 : coefficients(iecp + 1:iecp + nsemiloc) = sgp_potential%apot(1:nsemiloc, sl_l)
535 16 : powers(iecp + 1:iecp + nsemiloc) = sgp_potential%nrpot(1:nsemiloc, sl_l) - 2
536 8 : iecp = iecp + nsemiloc
537 : END DO
538 : END IF
539 : END IF
540 : END DO
541 :
542 : ! fail-safe check
543 2 : CPASSERT(iecp == ecp_num)
544 :
545 2 : rc = trexio_write_ecp_num(f, ecp_num)
546 2 : CALL trexio_error(rc)
547 :
548 2 : rc = trexio_write_ecp_z_core(f, z_core)
549 2 : CALL trexio_error(rc)
550 2 : DEALLOCATE (z_core)
551 :
552 2 : rc = trexio_write_ecp_max_ang_mom_plus_1(f, max_ang_mom_plus_1)
553 2 : CALL trexio_error(rc)
554 2 : DEALLOCATE (max_ang_mom_plus_1)
555 :
556 2 : rc = trexio_write_ecp_ang_mom(f, ang_mom)
557 2 : CALL trexio_error(rc)
558 2 : DEALLOCATE (ang_mom)
559 :
560 2 : rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
561 2 : CALL trexio_error(rc)
562 2 : DEALLOCATE (nucleus_index)
563 :
564 2 : rc = trexio_write_ecp_exponent(f, exponents)
565 2 : CALL trexio_error(rc)
566 2 : DEALLOCATE (exponents)
567 :
568 2 : rc = trexio_write_ecp_coefficient(f, coefficients)
569 2 : CALL trexio_error(rc)
570 2 : DEALLOCATE (coefficients)
571 :
572 2 : rc = trexio_write_ecp_power(f, powers)
573 2 : CALL trexio_error(rc)
574 2 : DEALLOCATE (powers)
575 : END IF
576 :
577 : END IF ! ionode
578 :
579 : !========================================================================================!
580 : ! Grid group
581 : !========================================================================================!
582 : ! TODO
583 :
584 : !========================================================================================!
585 : ! AO group
586 : !========================================================================================!
587 8 : CALL get_qs_env(qs_env, qs_kind_set=kind_set)
588 8 : CALL get_qs_kind_set(kind_set, ncgf=ncgf, nsgf=nsgf)
589 :
590 8 : CALL section_vals_val_get(trexio_section, "CARTESIAN", l_val=save_cartesian)
591 8 : IF (save_cartesian) THEN
592 4 : ao_num = ncgf
593 : ELSE
594 4 : ao_num = nsgf
595 : END IF
596 :
597 8 : IF (ionode) THEN
598 4 : IF (save_cartesian) THEN
599 2 : rc = trexio_write_ao_cartesian(f, 1)
600 : ELSE
601 2 : rc = trexio_write_ao_cartesian(f, 0)
602 : END IF
603 4 : CALL trexio_error(rc)
604 :
605 4 : rc = trexio_write_ao_num(f, ao_num)
606 4 : CALL trexio_error(rc)
607 : END IF
608 :
609 : ! one-to-one mapping between AOs and ...
610 24 : ALLOCATE (ao_shell(ao_num)) ! ..shells
611 24 : ALLOCATE (ao_normalization(ao_num)) ! ..normalization factors
612 :
613 : ! we need to be consistent with the basis group on the shell indices
614 8 : ishell = 0 ! global shell index
615 8 : igf = 0 ! global AO index
616 34 : DO iatom = 1, natoms
617 : ! get the qs_kind (index position in kind_set) for this atom (atomic_kind)
618 26 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
619 : ! get the primary (orbital) basis set associated to this qs_kind
620 26 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_set, basis_type="ORB")
621 : ! get the info from the basis set
622 : CALL get_gto_basis_set(basis_set, &
623 : nset=nset, &
624 : nshell=nshell, &
625 : norm_cgf=norm_cgf, &
626 : ncgf=ncgf, &
627 : nsgf=nsgf, &
628 26 : l=l_shell_set)
629 :
630 26 : icgf = 0
631 114 : DO iset = 1, nset
632 252 : DO ishell_loc = 1, nshell(iset)
633 : ! global shell index
634 138 : ishell = ishell + 1
635 : ! angular momentum l of this shell
636 138 : lshell = l_shell_set(ishell_loc, iset)
637 :
638 : ! number of AOs in this shell
639 138 : IF (save_cartesian) THEN
640 34 : ngf_shell = nco(lshell)
641 : ELSE
642 104 : ngf_shell = nso(lshell)
643 : END IF
644 :
645 : ! one-to-one mapping between AOs and shells
646 524 : ao_shell(igf + 1:igf + ngf_shell) = ishell
647 :
648 : ! one-to-one mapping between AOs and normalization factors
649 138 : IF (save_cartesian) THEN
650 136 : ao_normalization(igf + 1:igf + ngf_shell) = norm_cgf(icgf + 1:icgf + ngf_shell)
651 : ELSE
652 : ! allocate some temporary arrays
653 416 : ALLOCATE (diag_ncgf(nco(lshell), nco(lshell)))
654 416 : ALLOCATE (diag_nsgf(nso(lshell), nso(lshell)))
655 416 : ALLOCATE (temp(nso(lshell), nco(lshell)))
656 1808 : diag_ncgf = 0.0_dp
657 1436 : diag_nsgf = 0.0_dp
658 1616 : temp = 0.0_dp
659 :
660 416 : DO i = 1, nco(lshell)
661 416 : diag_ncgf(i, i) = norm_cgf(icgf + i)
662 : END DO
663 :
664 : ! transform the normalization factors from Cartesian to solid harmonics
665 42200 : temp(:, :) = MATMUL(orbtramat(lshell)%c2s, diag_ncgf)
666 24520 : diag_nsgf(:, :) = MATMUL(temp, TRANSPOSE(orbtramat(lshell)%s2c))
667 388 : DO i = 1, nso(lshell)
668 388 : ao_normalization(igf + i) = diag_nsgf(i, i)
669 : END DO
670 :
671 104 : DEALLOCATE (diag_ncgf)
672 104 : DEALLOCATE (diag_nsgf)
673 104 : DEALLOCATE (temp)
674 : END IF
675 :
676 138 : igf = igf + ngf_shell
677 226 : icgf = icgf + nco(lshell)
678 : END DO
679 : END DO
680 : ! just a failsafe check
681 86 : CPASSERT(icgf == ncgf)
682 : END DO
683 :
684 8 : IF (ionode) THEN
685 4 : rc = trexio_write_ao_shell(f, ao_shell)
686 4 : CALL trexio_error(rc)
687 :
688 4 : rc = trexio_write_ao_normalization(f, ao_normalization)
689 4 : CALL trexio_error(rc)
690 : END IF
691 :
692 8 : DEALLOCATE (ao_shell)
693 8 : DEALLOCATE (ao_normalization)
694 :
695 : !========================================================================================!
696 : ! MO group
697 : !========================================================================================!
698 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, kpoints=kpoints, dft_control=dft_control, &
699 8 : particle_set=particle_set, qs_kind_set=kind_set, blacs_env=blacs_env)
700 8 : nspins = dft_control%nspins
701 8 : CALL get_qs_kind_set(kind_set, nsgf=nsgf)
702 8 : nmo_spin = 0
703 :
704 : ! figure out that total number of MOs
705 8 : mo_num = 0
706 8 : IF (do_kpoints) THEN
707 2 : CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, use_real_wfn=use_real_wfn)
708 2 : CALL get_kpoint_env(kp_env(1)%kpoint_env, mos=mos_kp)
709 4 : DO ispin = 1, nspins
710 2 : CALL get_mo_set(mos_kp(1, ispin), nmo=nmo)
711 4 : nmo_spin(ispin) = nmo
712 : END DO
713 6 : mo_num = nkp*SUM(nmo_spin)
714 :
715 : ! we create a distributed fm matrix to gather the MOs from everywhere (in sph basis)
716 : CALL cp_fm_struct_create(fm_struct, para_env=para_env, context=blacs_env, &
717 2 : nrow_global=nsgf, ncol_global=mo_num)
718 2 : CALL cp_fm_create(fm_mo_coeff, fm_struct)
719 2 : CALL cp_fm_set_all(fm_mo_coeff, 0.0_dp)
720 2 : IF (.NOT. use_real_wfn) THEN
721 2 : CALL cp_fm_create(fm_mo_coeff_im, fm_struct)
722 2 : CALL cp_fm_set_all(fm_mo_coeff_im, 0.0_dp)
723 : END IF
724 2 : CALL cp_fm_struct_release(fm_struct)
725 : ELSE
726 6 : CALL get_qs_env(qs_env, mos=mos)
727 14 : DO ispin = 1, nspins
728 8 : CALL get_mo_set(mos(ispin), nmo=nmo)
729 14 : nmo_spin(ispin) = nmo
730 : END DO
731 18 : mo_num = SUM(nmo_spin)
732 : END IF
733 :
734 : ! allocate all the arrays
735 32 : ALLOCATE (mo_coefficient(ao_num, mo_num))
736 33432 : mo_coefficient(:, :) = 0.0_dp
737 24 : ALLOCATE (mo_energy(mo_num))
738 436 : mo_energy(:) = 0.0_dp
739 16 : ALLOCATE (mo_occupation(mo_num))
740 436 : mo_occupation(:) = 0.0_dp
741 24 : ALLOCATE (mo_spin(mo_num))
742 436 : mo_spin(:) = 0
743 : ! extra arrays for kpoints
744 8 : IF (do_kpoints) THEN
745 6 : ALLOCATE (mo_coefficient_im(ao_num, mo_num))
746 26882 : mo_coefficient_im(:, :) = 0.0_dp
747 4 : ALLOCATE (mo_kpoint(mo_num))
748 258 : mo_kpoint(:) = 0
749 : END IF
750 :
751 : ! in case of kpoints, we do this in 2 steps:
752 : ! 1. we gather the MOs of each kpt and pipe them into a single large distributed fm matrix;
753 : ! 2. we possibly transform the MOs of each kpt to Cartesian AOs and write them in the single large local array;
754 8 : IF (do_kpoints) THEN
755 2 : CALL get_kpoint_info(kpoints, kp_env=kp_env, nkp=nkp, kp_range=kp_range)
756 :
757 4 : DO ispin = 1, nspins
758 20 : DO ikp = 1, nkp
759 16 : nmo = nmo_spin(ispin)
760 : ! global index to store the MOs
761 16 : imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
762 :
763 : ! do I have this kpoint on this rank?
764 18 : IF (ikp >= kp_range(1) .AND. ikp <= kp_range(2)) THEN
765 16 : ikp_loc = ikp - kp_range(1) + 1
766 : ! get the mo set for this kpoint
767 16 : CALL get_kpoint_env(kp_env(ikp_loc)%kpoint_env, mos=mos_kp)
768 :
769 : ! if MOs are stored with dbcsr, copy them to fm
770 16 : IF (mos_kp(1, ispin)%use_mo_coeff_b) THEN
771 0 : CALL copy_dbcsr_to_fm(mos_kp(1, ispin)%mo_coeff_b, mos_kp(1, ispin)%mo_coeff)
772 : END IF
773 : ! copy real part of MO coefficients to large distributed fm matrix
774 : CALL cp_fm_to_fm_submat_general(mos_kp(1, ispin)%mo_coeff, fm_mo_coeff, &
775 16 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
776 :
777 : ! copy MO energies to local arrays
778 272 : mo_energy(imo + 1:imo + nmo) = mos_kp(1, ispin)%eigenvalues(1:nmo)
779 :
780 : ! copy MO occupations to local arrays
781 272 : mo_occupation(imo + 1:imo + nmo) = mos_kp(1, ispin)%occupation_numbers(1:nmo)
782 :
783 : ! same for the imaginary part of MO coefficients
784 16 : IF (.NOT. use_real_wfn) THEN
785 16 : IF (mos_kp(2, ispin)%use_mo_coeff_b) THEN
786 0 : CALL copy_dbcsr_to_fm(mos_kp(2, ispin)%mo_coeff_b, mos_kp(2, ispin)%mo_coeff)
787 : END IF
788 : CALL cp_fm_to_fm_submat_general(mos_kp(2, ispin)%mo_coeff, fm_mo_coeff_im, &
789 16 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
790 : END IF
791 : ELSE
792 : ! call with a dummy fm for receiving the data
793 : CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff, &
794 0 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
795 0 : IF (.NOT. use_real_wfn) THEN
796 : CALL cp_fm_to_fm_submat_general(fm_dummy, fm_mo_coeff_im, &
797 0 : nsgf, nmo, 1, 1, 1, imo + 1, blacs_env)
798 : END IF
799 : END IF
800 : END DO
801 : END DO
802 : END IF
803 :
804 : ! reduce MO energies and occupations to the master node
805 8 : IF (do_kpoints) THEN
806 2 : CALL get_kpoint_info(kpoints, para_env_inter_kp=para_env_inter_kp)
807 2 : CALL para_env_inter_kp%sum(mo_energy)
808 2 : CALL para_env_inter_kp%sum(mo_occupation)
809 : END IF
810 :
811 : ! second step: here we actually put everything in the local arrays for writing to trexio
812 18 : DO ispin = 1, nspins
813 : ! get number of MOs for this spin
814 10 : nmo = nmo_spin(ispin)
815 : ! allocate local temp array to transform the MOs of each kpoint/spin
816 40 : ALLOCATE (mos_sgf(nsgf, nmo))
817 :
818 10 : IF (do_kpoints) THEN
819 18 : DO ikp = 1, nkp
820 : ! global index to store the MOs
821 16 : imo = (ikp - 1)*nmo + (ispin - 1)*nmo_spin(1)*nkp
822 :
823 : ! store kpoint index
824 272 : mo_kpoint(imo + 1:imo + nmo) = ikp
825 : ! store the MO spins
826 272 : mo_spin(imo + 1:imo + nmo) = ispin - 1
827 :
828 : ! transform and store the MO coefficients
829 16 : CALL cp_fm_get_submatrix(fm_mo_coeff, mos_sgf, 1, imo + 1, nsgf, nmo)
830 16 : IF (save_cartesian) THEN
831 0 : CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
832 : ELSE
833 : ! we have to reorder the MOs since CP2K and TREXIO have different conventions
834 : ! from m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
835 : ! to m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
836 16 : i = 0
837 656 : DO ishell = 1, shell_num
838 640 : l = shell_ang_mom(ishell)
839 2304 : DO k = 1, 2*l + 1
840 : ! map from running index k: 1...2l+1 to magnetic quantum number m in TREXIO convention
841 1664 : m = (-1)**k*FLOOR(REAL(k, KIND=dp)/2.0_dp)
842 28928 : mo_coefficient(i + k, imo + 1:imo + nmo) = mos_sgf(i + l + m + 1, :)
843 : END DO
844 656 : i = i + 2*l + 1
845 : END DO
846 16 : CPASSERT(i == nsgf)
847 : END IF
848 :
849 : ! we have to do it for the imaginary part as well
850 18 : IF (.NOT. use_real_wfn) THEN
851 16 : CALL cp_fm_get_submatrix(fm_mo_coeff_im, mos_sgf, 1, imo + 1, nsgf, nmo)
852 16 : IF (save_cartesian) THEN
853 0 : CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient_im(:, imo + 1:imo + nmo))
854 : ELSE
855 : ! we have to reorder the MOs since CP2K and TREXIO have different conventions
856 : ! from m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
857 : ! to m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
858 16 : i = 0
859 656 : DO ishell = 1, shell_num
860 640 : l = shell_ang_mom(ishell)
861 2304 : DO k = 1, 2*l + 1
862 : ! map from running index k: 1...2l+1 to magnetic quantum number m in TREXIO convention
863 1664 : m = (-1)**k*FLOOR(REAL(k, KIND=dp)/2.0_dp)
864 28928 : mo_coefficient_im(i + k, imo + 1:imo + nmo) = mos_sgf(i + l + m + 1, :)
865 : END DO
866 656 : i = i + 2*l + 1
867 : END DO
868 16 : CPASSERT(i == nsgf)
869 : END IF
870 : END IF
871 : END DO
872 : ELSE ! no k-points
873 : ! global index to store the MOs
874 8 : imo = (ispin - 1)*nmo_spin(1)
875 : ! store the MO energies
876 180 : mo_energy(imo + 1:imo + nmo) = mos(ispin)%eigenvalues
877 : ! store the MO occupations
878 180 : mo_occupation(imo + 1:imo + nmo) = mos(ispin)%occupation_numbers
879 : ! store the MO spins
880 180 : mo_spin(imo + 1:imo + nmo) = ispin - 1
881 :
882 : ! check if we are using the dbcsr mo_coeff and copy them to fm if needed
883 8 : IF (mos(ispin)%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, mos(ispin)%mo_coeff)
884 :
885 : ! allocate a normal fortran array to store the spherical MO coefficients
886 8 : CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, mos_sgf)
887 :
888 8 : IF (save_cartesian) THEN
889 6 : CALL spherical_to_cartesian_mo(mos_sgf, particle_set, kind_set, mo_coefficient(:, imo + 1:imo + nmo))
890 : ELSE
891 : ! we have to reorder the MOs since CP2K and TREXIO have different conventions
892 : ! from m = -l, -l+1, ..., 0, ..., l-1, l of CP2K
893 : ! to m = 0, +1, -1, +2, -2, ..., +l, -l of TREXIO
894 2 : i = 0
895 26 : DO ishell = 1, shell_num
896 24 : l = shell_ang_mom(ishell)
897 100 : DO k = 1, 2*l + 1
898 : ! map from running index k: 1...2l+1 to magnetic quantum number m in TREXIO convention
899 76 : m = (-1)**k*FLOOR(REAL(k, KIND=dp)/2.0_dp)
900 2988 : mo_coefficient(i + k, imo + 1:imo + nmo) = mos_sgf(i + l + m + 1, :)
901 : END DO
902 26 : i = i + 2*l + 1
903 : END DO
904 2 : CPASSERT(i == nsgf)
905 : END IF
906 : END IF
907 :
908 18 : DEALLOCATE (mos_sgf)
909 : END DO
910 :
911 8 : IF (ionode) THEN
912 4 : rc = trexio_write_mo_type(f, 'Canonical', LEN_TRIM('Canonical') + 1)
913 4 : CALL trexio_error(rc)
914 :
915 4 : rc = trexio_write_mo_num(f, mo_num)
916 4 : CALL trexio_error(rc)
917 :
918 4 : rc = trexio_write_mo_coefficient(f, mo_coefficient)
919 4 : CALL trexio_error(rc)
920 :
921 4 : rc = trexio_write_mo_energy(f, mo_energy)
922 4 : CALL trexio_error(rc)
923 :
924 4 : rc = trexio_write_mo_occupation(f, mo_occupation)
925 4 : CALL trexio_error(rc)
926 :
927 4 : rc = trexio_write_mo_spin(f, mo_spin)
928 4 : CALL trexio_error(rc)
929 :
930 4 : IF (do_kpoints) THEN
931 1 : rc = trexio_write_mo_coefficient_im(f, mo_coefficient_im)
932 1 : CALL trexio_error(rc)
933 :
934 1 : rc = trexio_write_mo_k_point(f, mo_kpoint)
935 1 : CALL trexio_error(rc)
936 : END IF
937 : END IF
938 :
939 8 : DEALLOCATE (shell_ang_mom)
940 8 : DEALLOCATE (mo_coefficient)
941 8 : DEALLOCATE (mo_energy)
942 8 : DEALLOCATE (mo_occupation)
943 8 : DEALLOCATE (mo_spin)
944 8 : IF (do_kpoints) THEN
945 2 : DEALLOCATE (mo_coefficient_im)
946 2 : DEALLOCATE (mo_kpoint)
947 2 : CALL cp_fm_release(fm_mo_coeff)
948 2 : CALL cp_fm_release(fm_mo_coeff_im)
949 : END IF
950 :
951 : !========================================================================================!
952 : ! RDM group
953 : !========================================================================================!
954 : !TODO
955 :
956 : !========================================================================================!
957 : ! Close the TREXIO file
958 : !========================================================================================!
959 8 : IF (ionode) THEN
960 4 : rc = trexio_close(f)
961 4 : CALL trexio_error(rc)
962 : END IF
963 : #else
964 : MARK_USED(qs_env)
965 : MARK_USED(trexio_section)
966 : CPWARN('TREXIO support has not been enabled in this build.')
967 : #endif
968 8 : CALL timestop(handle)
969 :
970 48 : END SUBROUTINE write_trexio
971 :
972 : #ifdef __TREXIO
973 : ! **************************************************************************************************
974 : !> \brief Handles TREXIO errors
975 : !> \param rc the TREXIO return code
976 : ! **************************************************************************************************
977 195 : SUBROUTINE trexio_error(rc)
978 : INTEGER(trexio_exit_code), INTENT(IN) :: rc
979 :
980 : CHARACTER(LEN=128) :: err_msg
981 :
982 195 : IF (rc /= TREXIO_SUCCESS) THEN
983 0 : CALL trexio_string_of_error(rc, err_msg)
984 0 : CPABORT('TREXIO Error: '//TRIM(err_msg))
985 : END IF
986 :
987 195 : END SUBROUTINE trexio_error
988 :
989 : ! **************************************************************************************************
990 : !> \brief Computes the nuclear repulsion energy of a molecular system
991 : !> \param particle_set the set of particles in the system
992 : !> \param kind_set the set of qs_kinds in the system
993 : !> \param e_nn the nuclear repulsion energy
994 : ! **************************************************************************************************
995 2 : SUBROUTINE nuclear_repulsion_energy(particle_set, kind_set, e_nn)
996 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
997 : POINTER :: particle_set
998 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
999 : POINTER :: kind_set
1000 : REAL(KIND=dp), INTENT(OUT) :: e_nn
1001 :
1002 : INTEGER :: i, ikind, j, jkind, natoms
1003 : REAL(KIND=dp) :: r_ij, zeff_i, zeff_j
1004 :
1005 2 : natoms = SIZE(particle_set)
1006 2 : e_nn = 0.0_dp
1007 4 : DO i = 1, natoms
1008 2 : CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
1009 2 : CALL get_qs_kind(kind_set(ikind), zeff=zeff_i)
1010 4 : DO j = i + 1, natoms
1011 0 : r_ij = NORM2(particle_set(i)%r - particle_set(j)%r)
1012 :
1013 0 : CALL get_atomic_kind(particle_set(j)%atomic_kind, kind_number=jkind)
1014 0 : CALL get_qs_kind(kind_set(jkind), zeff=zeff_j)
1015 :
1016 2 : e_nn = e_nn + zeff_i*zeff_j/r_ij
1017 : END DO
1018 : END DO
1019 :
1020 2 : END SUBROUTINE nuclear_repulsion_energy
1021 :
1022 : ! **************************************************************************************************
1023 : !> \brief Computes a spherical to cartesian MO transformation (solid harmonics in reality)
1024 : !> \param mos_sgf the MO coefficients in spherical AO basis
1025 : !> \param particle_set the set of particles in the system
1026 : !> \param qs_kind_set the set of qs_kinds in the system
1027 : !> \param mos_cgf the transformed MO coefficients in Cartesian AO basis
1028 : ! **************************************************************************************************
1029 6 : SUBROUTINE spherical_to_cartesian_mo(mos_sgf, particle_set, qs_kind_set, mos_cgf)
1030 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: mos_sgf
1031 : TYPE(particle_type), DIMENSION(:), INTENT(IN), &
1032 : POINTER :: particle_set
1033 : TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
1034 : POINTER :: qs_kind_set
1035 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: mos_cgf
1036 :
1037 : INTEGER :: iatom, icgf, ikind, iset, isgf, ishell, &
1038 : lshell, ncgf, nmo, nset, nsgf
1039 6 : INTEGER, DIMENSION(:), POINTER :: nshell
1040 6 : INTEGER, DIMENSION(:, :), POINTER :: l
1041 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1042 :
1043 6 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
1044 :
1045 3586 : mos_cgf = 0.0_dp
1046 6 : nmo = SIZE(mos_sgf, 2)
1047 :
1048 : ! Transform spherical MOs to Cartesian MOs
1049 6 : icgf = 1
1050 6 : isgf = 1
1051 20 : DO iatom = 1, SIZE(particle_set)
1052 14 : NULLIFY (orb_basis_set)
1053 14 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1054 14 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1055 :
1056 34 : IF (ASSOCIATED(orb_basis_set)) THEN
1057 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1058 : nset=nset, &
1059 : nshell=nshell, &
1060 14 : l=l)
1061 54 : DO iset = 1, nset
1062 98 : DO ishell = 1, nshell(iset)
1063 44 : lshell = l(ishell, iset)
1064 : CALL dgemm("T", "N", nco(lshell), nmo, nso(lshell), 1.0_dp, &
1065 : orbtramat(lshell)%c2s, nso(lshell), &
1066 : mos_sgf(isgf, 1), nsgf, 0.0_dp, &
1067 44 : mos_cgf(icgf, 1), ncgf)
1068 44 : icgf = icgf + nco(lshell)
1069 84 : isgf = isgf + nso(lshell)
1070 : END DO
1071 : END DO
1072 : ELSE
1073 : ! assume atom without basis set
1074 0 : CPABORT("Unknown basis set type")
1075 : END IF
1076 : END DO ! iatom
1077 :
1078 6 : END SUBROUTINE spherical_to_cartesian_mo
1079 : #endif
1080 :
1081 : END MODULE trexio_utils
|
