Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief The module to read/write QCSchema HDF5 files for interfacing CP2K with other programs
10 : !> \par History
11 : !> 10.2022 created [SB]
12 : !> \author Stefano Battaglia
13 : ! **************************************************************************************************
14 : MODULE qcschema
15 :
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE basis_set_types, ONLY: gto_basis_set_type
18 : USE cp2k_info, ONLY: cp2k_version
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger, &
21 : cp_logger_get_default_io_unit, &
22 : cp_logger_type
23 : #ifdef __HDF5
24 : USE hdf5_wrapper, ONLY: &
25 : h5aread_double_scalar, h5awrite_boolean, h5awrite_double_scalar, h5awrite_double_simple, &
26 : h5awrite_fixlen_string, h5awrite_integer_scalar, h5awrite_integer_simple, &
27 : h5awrite_string_simple, h5close, h5dread_double_simple, h5dwrite_double_simple, h5fclose, &
28 : h5fcreate, h5fopen, h5gclose, h5gcreate, h5gopen, h5open, hdf5_id
29 : #endif
30 : USE input_section_types, ONLY: section_vals_get, &
31 : section_vals_get_subs_vals, &
32 : section_vals_type
33 : USE kinds, ONLY: default_path_length, &
34 : default_string_length, &
35 : dp, &
36 : int_8
37 : USE mp2_types, ONLY: mp2_type
38 : USE particle_types, ONLY: particle_type
39 : USE periodic_table, ONLY: get_ptable_info
40 : USE qs_active_space_types, ONLY: active_space_type
41 : USE qs_active_space_utils, ONLY: eri_to_array, &
42 : subspace_matrix_to_array
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_force_types, ONLY: qs_force_type
47 : USE qs_kind_types, ONLY: get_qs_kind, &
48 : qs_kind_type
49 : USE qs_ks_types, ONLY: qs_ks_env_type
50 : USE qs_scf_types, ONLY: qs_scf_env_type
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qcschema'
58 :
59 : PUBLIC :: qcschema_type
60 : PUBLIC :: qcschema_env_create, qcschema_env_release, qcschema_to_hdf5
61 :
62 : ! **************************************************************************************************
63 : !> \brief A derived type to store the program information that generated the QCSchema file.
64 : !> For more information refer to:
65 : !> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#provenance
66 : ! **************************************************************************************************
67 : TYPE qcschema_provenance
68 : CHARACTER(LEN=default_string_length) :: creator = "" ! The name of the creator of this object
69 : CHARACTER(LEN=default_string_length) :: version = "" ! The version of the creator of this object
70 : CHARACTER(LEN=default_string_length) :: routine = "" ! The routine that was used to create this object
71 : END TYPE qcschema_provenance
72 :
73 : ! **************************************************************************************************
74 : !> \brief A derived type to store the topological information of the physical system.
75 : !> For more information refer to:
76 : !> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#topology
77 : ! **************************************************************************************************
78 : TYPE qcschema_topology
79 : CHARACTER(LEN=default_string_length) :: name = "" ! of the molecule
80 : CHARACTER(LEN=2), DIMENSION(:), ALLOCATABLE :: symbols ! of the atoms
81 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: geometry ! row major, in bohr
82 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: masses
83 : INTEGER, DIMENSION(:), ALLOCATABLE :: atomic_numbers
84 : INTEGER :: molecular_charge = 0
85 : INTEGER :: molecular_multiplicity = 1
86 : CHARACTER(LEN=default_string_length) :: schema_name = ""
87 : INTEGER :: schema_version = 0
88 : TYPE(qcschema_provenance) :: provenance = qcschema_provenance()
89 : END TYPE qcschema_topology
90 :
91 : ! **************************************************************************************************
92 : !> \brief A derived type to store the information of a single electron shell in a basis set.
93 : !> For more information refer to:
94 : !> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L43
95 : ! **************************************************************************************************
96 : TYPE qcschema_electron_shell
97 : ! The angular momenta of this electron shell as a list of integers
98 : INTEGER, DIMENSION(:), POINTER :: angular_momentum => NULL()
99 : ! The type of this shell: spherical or cartesian
100 : CHARACTER(LEN=9) :: harmonic_type = ""
101 : ! The exponents of this contracted shell. The official spec stores these values as strings
102 : REAL(KIND=dp), DIMENSION(:), POINTER :: exponents => NULL()
103 : ! The general contraction coefficients of this contracted shell
104 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: coefficients => NULL()
105 : END TYPE qcschema_electron_shell
106 :
107 : ! **************************************************************************************************
108 : !> \brief A derived type to store the information of an ECP in a basis set.
109 : !> For more information refer to:
110 : !> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L90
111 : ! **************************************************************************************************
112 : TYPE qcschema_ecp
113 : ! The type of this potential
114 : CHARACTER(LEN=default_string_length) :: ecp_type = ""
115 : ! The angular momenta of this potential as a list of integers
116 : INTEGER, DIMENSION(:), POINTER :: angular_momentum => NULL()
117 : ! The exponents of the r terms
118 : INTEGER, DIMENSION(:), POINTER :: r_exponents => NULL()
119 : ! The exponents of the Gaussian terms
120 : REAL(KIND=dp), DIMENSION(:), POINTER :: gaussian_exponents => NULL()
121 : ! The general contraction coefficients of this potential
122 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: coefficients => NULL()
123 : END TYPE qcschema_ecp
124 :
125 : ! **************************************************************************************************
126 : !> \brief A derived type to store the information of a single atom/center in the basis.
127 : !> For more information refer to:
128 : !> https://github.com/MolSSI/QCSchema/blob/1d5ff3baa5/qcschema/dev/definitions.py#L146
129 : ! **************************************************************************************************
130 : TYPE qcschema_center_basis
131 : ! The list of electronic shells for this element
132 : TYPE(qcschema_electron_shell), DIMENSION(:), POINTER :: electron_shells => NULL()
133 : ! The list of effective core potentials for this element
134 : TYPE(qcschema_ecp), DIMENSION(:), POINTER :: ecp_potentials => NULL()
135 : ! The number of electrons replaced by an ECP
136 : INTEGER :: ecp_electrons = 0
137 : END TYPE qcschema_center_basis
138 :
139 : ! **************************************************************************************************
140 : !> \brief A derived type to store the information of the basis set used in the calculation.
141 : !> For more information refer to:
142 : !> https://molssi-qc-schema.readthedocs.io/en/latest/auto_basis.html#basis-set-schema
143 : ! **************************************************************************************************
144 : TYPE qcschema_basis_set
145 : ! The name of the basis set
146 : CHARACTER(LEN=default_string_length) :: name = ""
147 : ! A dictionary mapping the keys provided by `atom_map` to their basis center data
148 : TYPE(qcschema_center_basis), DIMENSION(:), POINTER :: center_data => NULL()
149 : ! The list of atomic kinds, indicating the keys used to store the basis in `center_data`
150 : ! Not clear if this will be of the length of the basis set size, or rather just one
151 : ! entry for atomic kind. E.g. only one entry for hydrogen even though there might be
152 : ! many hydrogen atoms in the molecule. If this is the case, then we really need a
153 : ! hash table for `center_data`
154 : CHARACTER(LEN=2), DIMENSION(:), POINTER :: atom_map => NULL()
155 : ! The version of this specific schema
156 : INTEGER :: schema_version = -1
157 : ! The name of this schema. This value is expected to be `qcschema_basis`
158 : CHARACTER(LEN=default_string_length) :: schema_name = ""
159 : ! A description of this basis set
160 : CHARACTER(LEN=default_string_length) :: description = ""
161 : END TYPE qcschema_basis_set
162 :
163 : ! **************************************************************************************************
164 : !> \brief A derived type to store any additional computed wavefunction properties.
165 : !> Matrix quantities are stored as flat, column-major arrays.
166 : !> For more information refer to:
167 : !> https://molssi-qc-schema.readthedocs.io/en/latest/auto_wf.html#wavefunction-schema
168 : ! **************************************************************************************************
169 : TYPE qcschema_wavefunction
170 :
171 : ! The name of the method used to obtain the wf
172 : CHARACTER(LEN=default_string_length) :: method = ""
173 :
174 : ! The basis set used during the computation
175 : TYPE(qcschema_basis_set) :: basis_set = qcschema_basis_set()
176 :
177 : ! SCF quantities in AO or MO basis
178 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_orbitals_a
179 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_orbitals_b
180 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eigenvalues_a
181 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eigenvalues_b
182 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_occupations_a
183 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_occupations_b
184 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_density_mo_a
185 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_density_mo_b
186 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_fock_mo_a
187 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_fock_mo_b
188 :
189 : ! Electron repulsion integrals
190 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri
191 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_aa
192 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_ab
193 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: scf_eri_mo_bb
194 :
195 : ! Quantities with localized orbitals. All `nmo` orbitals are included,
196 : ! even if only a subset were localized
197 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_orbitals_a
198 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_orbitals_b
199 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_fock_a
200 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: localized_fock_b
201 :
202 : ! Whether the computation used restricted spin orbitals
203 : LOGICAL :: restricted = .FALSE.
204 :
205 : END TYPE qcschema_wavefunction
206 :
207 : ! **************************************************************************************************
208 : !> \brief A derived type to store the computed properties of the original calculation.
209 : !> For more information refer to:
210 : !> https://molssi-qc-schema.readthedocs.io/en/latest/auto_props.html#properties-schema
211 : ! **************************************************************************************************
212 : TYPE qcschema_properties
213 :
214 : REAL(KIND=dp) :: return_energy = 0.0_dp
215 :
216 : INTEGER :: calcinfo_nbasis = 0 ! AO basis size
217 : INTEGER :: calcinfo_nmo = 0 ! MO basis size
218 : INTEGER :: calcinfo_nalpha = 0 ! # of alpha electrons
219 : INTEGER :: calcinfo_nbeta = 0 ! # of beta electrons
220 : INTEGER :: calcinfo_natom = 0
221 :
222 : ! SCF results
223 : INTEGER :: scf_iterations = 0
224 : REAL(KIND=dp) :: scf_one_electron_energy = 0.0_dp
225 : REAL(KIND=dp) :: scf_two_electron_energy = 0.0_dp
226 : REAL(KIND=dp) :: nuclear_repulsion_energy = 0.0_dp
227 : REAL(KIND=dp) :: scf_vv10_energy = 0.0_dp
228 : REAL(KIND=dp) :: scf_xc_energy = 0.0_dp
229 : REAL(KIND=dp) :: scf_dispersion_correction_energy = 0.0_dp
230 : REAL(KIND=dp) :: scf_total_energy = 0.0_dp
231 : ! the dipole moment is calculated on the fly and not stored
232 : REAL(KIND=dp), DIMENSION(3) :: scf_dipole_moment = 0.0_dp
233 :
234 : ! MP2 results
235 : REAL(KIND=dp) :: mp2_same_spin_correlation_energy = 0.0_dp
236 : REAL(KIND=dp) :: mp2_opposite_spin_correlation_energy = 0.0_dp
237 : REAL(KIND=dp) :: mp2_singles_energy = 0.0_dp
238 : REAL(KIND=dp) :: mp2_doubles_energy = 0.0_dp
239 : ! these are the only two that are saved
240 : REAL(KIND=dp) :: mp2_correlation_energy = 0.0_dp
241 : REAL(KIND=dp) :: mp2_total_energy = 0.0_dp
242 :
243 : ! internal flags to know the type of calculation
244 : LOGICAL :: mp2 = .FALSE.
245 :
246 : END TYPE qcschema_properties
247 :
248 : ! **************************************************************************************************
249 : !> \brief The full QCSchema output type.
250 : !> For more information refer to:
251 : !> https://molssi-qc-schema.readthedocs.io/en/latest/spec_components.html#output-components
252 : ! **************************************************************************************************
253 : TYPE qcschema_type
254 : TYPE(qcschema_topology) :: topology = qcschema_topology()
255 : TYPE(qcschema_provenance) :: provenance = qcschema_provenance()
256 : TYPE(qcschema_properties) :: properties = qcschema_properties()
257 : TYPE(qcschema_wavefunction) :: wavefunction = qcschema_wavefunction()
258 : TYPE(qcschema_basis_set) :: basis = qcschema_basis_set()
259 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: return_result
260 : CHARACTER(LEN=default_string_length) :: driver = ""
261 : LOGICAL :: success = .FALSE.
262 : END TYPE qcschema_type
263 :
264 : CONTAINS
265 :
266 : ! **************************************************************************************************
267 : !> \brief Create and initialize a qcschema object from a quickstep environment
268 : !> \param qcschema_env the qcschema environment to populate
269 : !> \param qs_env the qs environment with all the info of the computation
270 : ! **************************************************************************************************
271 4 : SUBROUTINE qcschema_env_create(qcschema_env, qs_env)
272 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
273 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
274 :
275 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qcschema_env_create'
276 :
277 : CHARACTER(LEN=2) :: atomic_symbol
278 : CHARACTER(LEN=default_string_length) :: basis_set_name, method
279 : INTEGER :: atomic_number, handle, i, i_glb, iatom, &
280 : ikind, nalpha, nao, natoms, nbeta, &
281 : nel, nmo, nspins, output_unit
282 : LOGICAL :: do_hfx
283 : REAL(KIND=dp) :: dispersion, mass, one_el_en, two_el_en
284 : TYPE(active_space_type), POINTER :: active_space_env
285 : TYPE(cp_logger_type), POINTER :: logger
286 : TYPE(dft_control_type), POINTER :: dft_control
287 : TYPE(gto_basis_set_type), POINTER :: basis_set
288 : TYPE(mp2_type), POINTER :: mp2_env
289 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
290 : TYPE(qs_energy_type), POINTER :: energy
291 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
292 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
293 : TYPE(qs_ks_env_type), POINTER :: ks_env
294 : TYPE(qs_scf_env_type), POINTER :: scf_env
295 : TYPE(section_vals_type), POINTER :: hfx_sections, input
296 :
297 4 : CALL timeset(routineN, handle)
298 :
299 4 : logger => cp_get_default_logger()
300 4 : output_unit = cp_logger_get_default_io_unit(logger)
301 :
302 : ! reset everything
303 4 : CALL qcschema_env_release(qcschema_env)
304 :
305 : ! collect environment info
306 4 : IF (ASSOCIATED(qs_env)) THEN
307 : CALL get_qs_env(qs_env, ks_env=ks_env, energy=energy, &
308 : dft_control=dft_control, force=force, &
309 : particle_set=particle_set, &
310 : scf_env=scf_env, mp2_env=mp2_env, &
311 : input=input, qs_kind_set=kind_set, &
312 4 : active_space=active_space_env)
313 : ELSE
314 0 : CPABORT("QS environment not associated, QCSchema interface quitting")
315 : END IF
316 :
317 : ! we need the AS environemnt to get all the SCF data
318 4 : IF (.NOT. ASSOCIATED(active_space_env)) THEN
319 0 : CPABORT("Active space environment not associated, QCSchema interface quitting")
320 : END IF
321 :
322 : !========================================================================================!
323 : ! *** QCSchema provenance ***
324 : !========================================================================================!
325 :
326 4 : qcschema_env%provenance%creator = 'CP2K'
327 4 : qcschema_env%provenance%version = cp2k_version
328 4 : qcschema_env%provenance%routine = routineN
329 :
330 : !========================================================================================!
331 : ! *** QCSchema topology ***
332 : !========================================================================================!
333 :
334 4 : qcschema_env%topology%schema_name = 'qcschema'
335 4 : qcschema_env%topology%schema_version = 3
336 :
337 4 : natoms = SIZE(particle_set)
338 :
339 12 : ALLOCATE (qcschema_env%topology%geometry(3*natoms))
340 8 : ALLOCATE (qcschema_env%topology%symbols(natoms))
341 12 : ALLOCATE (qcschema_env%topology%atomic_numbers(natoms))
342 12 : ALLOCATE (qcschema_env%topology%masses(natoms))
343 :
344 12 : DO iatom = 1, natoms
345 : ! set the geometry as a flat array
346 32 : qcschema_env%topology%geometry((iatom - 1)*3 + 1:(iatom)*3) = particle_set(iatom)%r(1:3)
347 :
348 : ! set the atomic symbols
349 8 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=atomic_symbol)
350 8 : qcschema_env%topology%symbols(iatom) = atomic_symbol
351 :
352 : ! set the atomic numbers and masses
353 8 : CALL get_ptable_info(atomic_symbol, number=atomic_number, amass=mass)
354 8 : qcschema_env%topology%atomic_numbers(iatom) = atomic_number
355 12 : qcschema_env%topology%masses(iatom) = mass
356 : END DO
357 :
358 4 : qcschema_env%topology%molecular_charge = dft_control%charge
359 4 : qcschema_env%topology%molecular_multiplicity = dft_control%multiplicity
360 :
361 : !========================================================================================!
362 : ! *** QCSchema properties ***
363 : !========================================================================================!
364 :
365 4 : nspins = active_space_env%nspins
366 :
367 4 : nao = active_space_env%mos_active(1)%nao
368 4 : nmo = active_space_env%nmo_active
369 4 : nel = active_space_env%nelec_active
370 :
371 4 : IF (nspins == 1) THEN
372 4 : nalpha = active_space_env%nelec_active/2
373 4 : nbeta = nalpha
374 : ELSE
375 0 : nalpha = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2
376 0 : nbeta = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2
377 : END IF
378 :
379 4 : qcschema_env%properties%calcinfo_natom = natoms
380 4 : qcschema_env%properties%calcinfo_nbasis = nao
381 4 : qcschema_env%properties%calcinfo_nmo = nmo
382 4 : qcschema_env%properties%calcinfo_nalpha = nalpha
383 4 : qcschema_env%properties%calcinfo_nbeta = nbeta
384 :
385 : ! energy results
386 4 : qcschema_env%properties%return_energy = energy%total
387 4 : qcschema_env%properties%scf_total_energy = energy%total
388 : ! here we abuse the nuclear repulsion energy to store the inactive energy
389 4 : qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive
390 : ! SCF info
391 4 : qcschema_env%properties%scf_iterations = scf_env%iter_count
392 : ! one-electron energy is the sum of all core terms
393 4 : one_el_en = energy%core_overlap + energy%core_self + energy%core
394 4 : qcschema_env%properties%scf_two_electron_energy = one_el_en
395 : ! two-electron energy is the sum of hartree and exact exchange (if there)
396 4 : two_el_en = energy%hartree + energy%ex + energy%hartree_1c
397 4 : qcschema_env%properties%scf_one_electron_energy = two_el_en
398 : ! xc energy
399 : qcschema_env%properties%scf_xc_energy = &
400 4 : energy%exc + energy%exc_aux_fit + energy%exc1 + energy%exc1_aux_fit
401 : ! dispersion energy
402 4 : dispersion = energy%dispersion + energy%gcp
403 4 : qcschema_env%properties%scf_dispersion_correction_energy = dispersion
404 :
405 : ! Some methods of CP2K are not supported by QCSchema, let's warn the user
406 4 : IF (dft_control%smear) CPABORT('WARNING: smearing not supported in QCSchema')
407 4 : IF (dft_control%dft_plus_u) CPABORT('WARNING: DFT+U not supported in QCSchema')
408 4 : IF (dft_control%do_sccs) CPABORT('WARNING: SCCS not supported in QCSchema')
409 4 : IF (qs_env%qmmm) CPABORT('WARNING: QM/MM not supported in QCSchema')
410 4 : IF (dft_control%qs_control%mulliken_restraint) THEN
411 0 : CPABORT('WARNING: Mulliken restrains not supported in QCSchema')
412 : END IF
413 4 : IF (dft_control%qs_control%semi_empirical) THEN
414 0 : CPABORT('WARNING: semi_empirical methods not supported in QCSchema')
415 : END IF
416 4 : IF (dft_control%qs_control%dftb) CPABORT('WARNING: DFTB not supported in QCSchema')
417 4 : IF (dft_control%qs_control%xtb) CPABORT('WARNING: xTB not supported in QCSchema')
418 :
419 : ! MP2 info
420 4 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
421 0 : qcschema_env%properties%mp2 = .TRUE.
422 : ! this info is computed on the fly, but not stored!
423 : ! qcschema_env%properties%mp2_same_spin_correlation_energy
424 : ! qcschema_env%properties%mp2_opposite_spin_correlation_energy
425 :
426 0 : qcschema_env%properties%mp2_correlation_energy = energy%mp2
427 0 : qcschema_env%properties%mp2_total_energy = energy%total
428 :
429 : ! update the scf energy
430 0 : qcschema_env%properties%scf_total_energy = energy%total - energy%mp2
431 : END IF
432 :
433 : !========================================================================================!
434 : ! *** QCSchema wavefunction ***
435 : !========================================================================================!
436 :
437 4 : IF (nspins == 1) THEN
438 4 : qcschema_env%wavefunction%restricted = .TRUE.
439 : ELSE
440 0 : qcschema_env%wavefunction%restricted = .FALSE.
441 : END IF
442 :
443 : ! alpha MO energies
444 12 : ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a(nmo))
445 22 : DO i = 1, nmo
446 18 : i_glb = active_space_env%active_orbitals(i, 1)
447 : qcschema_env%wavefunction%scf_eigenvalues_a(i) = &
448 22 : active_space_env%mos_active(1)%eigenvalues(i_glb)
449 : END DO
450 :
451 : ! alpha MO occupations
452 8 : ALLOCATE (qcschema_env%wavefunction%scf_occupations_a(nmo))
453 22 : DO i = 1, nmo
454 18 : i_glb = active_space_env%active_orbitals(i, 1)
455 : qcschema_env%wavefunction%scf_occupations_a(i) = &
456 22 : active_space_env%mos_active(1)%occupation_numbers(i_glb)
457 : END DO
458 :
459 : ! alpha Fock matrix
460 12 : ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a(nmo*nmo))
461 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), &
462 : qcschema_env%wavefunction%scf_fock_mo_a, &
463 : active_space_env%active_orbitals(:, 1), &
464 4 : active_space_env%active_orbitals(:, 1))
465 :
466 : ! alpha density matrix
467 8 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo))
468 : CALL subspace_matrix_to_array(active_space_env%p_active(1), &
469 : qcschema_env%wavefunction%scf_density_mo_a, &
470 : active_space_env%active_orbitals(:, 1), &
471 4 : active_space_env%active_orbitals(:, 1))
472 :
473 : ! alpha MOs coefficients
474 12 : ALLOCATE (qcschema_env%wavefunction%scf_orbitals_a(nao*nmo))
475 : CALL subspace_matrix_to_array(active_space_env%mos_active(1)%mo_coeff, &
476 : qcschema_env%wavefunction%scf_orbitals_a, &
477 84 : [(i, i=1, nao)], active_space_env%active_orbitals(:, 1))
478 :
479 4 : IF (nspins == 2) THEN
480 : ! beta MO energies
481 0 : ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b(nmo))
482 0 : DO i = 1, nmo
483 0 : i_glb = active_space_env%active_orbitals(i, 2)
484 : qcschema_env%wavefunction%scf_eigenvalues_b(i) = &
485 0 : active_space_env%mos_active(2)%eigenvalues(i_glb)
486 : END DO
487 :
488 : ! beta MO occupations
489 0 : ALLOCATE (qcschema_env%wavefunction%scf_occupations_b(nmo))
490 0 : DO i = 1, nmo
491 0 : i_glb = active_space_env%active_orbitals(i, 2)
492 : qcschema_env%wavefunction%scf_occupations_b(i) = &
493 0 : active_space_env%mos_active(2)%occupation_numbers(i_glb)
494 : END DO
495 :
496 : ! beta Fock matrix
497 0 : ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b(nmo*nmo))
498 : CALL subspace_matrix_to_array(active_space_env%fock_sub(2), &
499 : qcschema_env%wavefunction%scf_fock_mo_b, &
500 : active_space_env%active_orbitals(:, 2), &
501 0 : active_space_env%active_orbitals(:, 2))
502 :
503 : ! beta density matrix
504 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo))
505 : CALL subspace_matrix_to_array(active_space_env%p_active(2), &
506 : qcschema_env%wavefunction%scf_density_mo_b, &
507 : active_space_env%active_orbitals(:, 2), &
508 0 : active_space_env%active_orbitals(:, 2))
509 :
510 : ! beta MOs coefficients
511 0 : ALLOCATE (qcschema_env%wavefunction%scf_orbitals_b(nao*nmo))
512 : CALL subspace_matrix_to_array(active_space_env%mos_active(2)%mo_coeff, &
513 : qcschema_env%wavefunction%scf_orbitals_b, &
514 0 : [(i, i=1, nao)], active_space_env%active_orbitals(:, 2))
515 : END IF
516 :
517 : ! get the alpha-alpha eri
518 12 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa(nmo**4))
519 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_aa, &
520 4 : active_space_env%active_orbitals, 1, 1)
521 :
522 4 : IF (nspins == 2) THEN
523 : ! get the alpha-beta eri
524 0 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab(nmo**4))
525 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_ab, &
526 0 : active_space_env%active_orbitals, 1, 2)
527 :
528 : ! get the beta-beta eri
529 0 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb(nmo**4))
530 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_bb, &
531 0 : active_space_env%active_orbitals, 2, 2)
532 : END IF
533 :
534 : !========================================================================================!
535 : ! *** QCSchema model ***
536 : !========================================================================================!
537 :
538 12 : DO iatom = 1, natoms
539 8 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
540 8 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_set)
541 :
542 8 : basis_set_name = basis_set%name
543 :
544 : ! make sure that we do not run a mixed basis set
545 20 : IF (iatom > 1) THEN
546 4 : CPASSERT(basis_set_name == basis_set%name)
547 : END IF
548 : END DO
549 4 : qcschema_env%wavefunction%basis_set%name = basis_set_name
550 :
551 : ! figure out which method was used for the calculation
552 4 : IF (dft_control%uks) THEN
553 0 : method = 'U'
554 4 : ELSE IF (dft_control%roks) THEN
555 0 : method = 'RO'
556 : ELSE
557 4 : method = 'R'
558 : END IF
559 :
560 4 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
561 4 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
562 :
563 4 : IF (do_hfx) THEN
564 2 : method = TRIM(method)//'HF'
565 2 : ELSE IF (qcschema_env%properties%mp2) THEN
566 0 : method = TRIM(method)//'MP2'
567 : ELSE
568 2 : method = TRIM(method)//'KS'
569 : END IF
570 :
571 4 : qcschema_env%wavefunction%method = TRIM(method)
572 :
573 : !========================================================================================!
574 : ! *** QCSchema root ***
575 : !========================================================================================!
576 :
577 : ! driver
578 4 : IF (ASSOCIATED(force)) THEN
579 0 : qcschema_env%driver = 'gradient'
580 : ELSE
581 4 : qcschema_env%driver = 'energy'
582 : END IF
583 :
584 : ! success
585 : ! TODO: how to check if the calculation was succesful?
586 4 : qcschema_env%success = .TRUE.
587 :
588 : ! return result
589 : IF (qcschema_env%success) THEN
590 4 : IF (qcschema_env%driver == 'energy') THEN
591 4 : ALLOCATE (qcschema_env%return_result(1))
592 4 : qcschema_env%return_result(1) = energy%total
593 : ELSE
594 0 : ALLOCATE (qcschema_env%return_result(3*SIZE(particle_set)))
595 : ! TODO: populate with forces!!
596 0 : qcschema_env%return_result = 0.0_dp
597 : END IF
598 : ELSE
599 : CPABORT("The calculation to build the AS is unsuccessful")
600 : END IF
601 :
602 4 : CALL timestop(handle)
603 :
604 8 : END SUBROUTINE qcschema_env_create
605 :
606 : ! **************************************************************************************************
607 : !> \brief Releases the allocated memory of a qcschema environment
608 : !> \param qcschema_env the qcschema environment to release
609 : ! **************************************************************************************************
610 8 : SUBROUTINE qcschema_env_release(qcschema_env)
611 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
612 :
613 8 : IF (ALLOCATED(qcschema_env%return_result)) THEN
614 4 : DEALLOCATE (qcschema_env%return_result)
615 : END IF
616 :
617 8 : IF (ALLOCATED(qcschema_env%topology%atomic_numbers)) THEN
618 4 : DEALLOCATE (qcschema_env%topology%atomic_numbers)
619 : END IF
620 :
621 8 : IF (ALLOCATED(qcschema_env%topology%masses)) THEN
622 4 : DEALLOCATE (qcschema_env%topology%masses)
623 : END IF
624 :
625 8 : IF (ALLOCATED(qcschema_env%topology%geometry)) THEN
626 4 : DEALLOCATE (qcschema_env%topology%geometry)
627 : END IF
628 :
629 8 : IF (ALLOCATED(qcschema_env%topology%symbols)) THEN
630 4 : DEALLOCATE (qcschema_env%topology%symbols)
631 : END IF
632 :
633 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN
634 4 : DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_a)
635 : END IF
636 :
637 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN
638 0 : DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_b)
639 : END IF
640 :
641 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_a)) THEN
642 4 : DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a)
643 : END IF
644 :
645 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_b)) THEN
646 0 : DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b)
647 : END IF
648 :
649 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_a)) THEN
650 4 : DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_a)
651 : END IF
652 :
653 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_b)) THEN
654 0 : DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_b)
655 : END IF
656 :
657 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_a)) THEN
658 4 : DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a)
659 : END IF
660 :
661 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_b)) THEN
662 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b)
663 : END IF
664 :
665 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_a)) THEN
666 4 : DEALLOCATE (qcschema_env%wavefunction%scf_occupations_a)
667 : END IF
668 :
669 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_b)) THEN
670 0 : DEALLOCATE (qcschema_env%wavefunction%scf_occupations_b)
671 : END IF
672 :
673 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri)) THEN
674 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri)
675 : END IF
676 :
677 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_aa)) THEN
678 4 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa)
679 : END IF
680 :
681 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_bb)) THEN
682 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb)
683 : END IF
684 :
685 8 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_ab)) THEN
686 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab)
687 : END IF
688 :
689 8 : IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_a)) THEN
690 0 : DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_a)
691 : END IF
692 :
693 8 : IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_b)) THEN
694 0 : DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_b)
695 : END IF
696 :
697 8 : IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_a)) THEN
698 0 : DEALLOCATE (qcschema_env%wavefunction%localized_fock_a)
699 : END IF
700 :
701 8 : IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_b)) THEN
702 0 : DEALLOCATE (qcschema_env%wavefunction%localized_fock_b)
703 : END IF
704 :
705 8 : END SUBROUTINE qcschema_env_release
706 :
707 : ! **************************************************************************************************
708 : !> \brief Updates the Fock matrix and the inactive energy in a qcschema object
709 : !> \param qcschema_env the qcschema environment
710 : !> \param active_space_env the active space environment with the updated data
711 : ! **************************************************************************************************
712 0 : SUBROUTINE qcschema_update_fock(qcschema_env, active_space_env)
713 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
714 : TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env
715 :
716 : ! alpha Fock matrix
717 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), &
718 : qcschema_env%wavefunction%scf_fock_mo_a, &
719 : active_space_env%active_orbitals(:, 1), &
720 0 : active_space_env%active_orbitals(:, 1))
721 :
722 : ! beta Fock matrix
723 0 : IF (active_space_env%nspins == 2) THEN
724 : CALL subspace_matrix_to_array(active_space_env%fock_sub(2), &
725 : qcschema_env%wavefunction%scf_fock_mo_b, &
726 : active_space_env%active_orbitals(:, 2), &
727 0 : active_space_env%active_orbitals(:, 2))
728 : END IF
729 :
730 : ! update inactive energy
731 0 : qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive
732 :
733 0 : END SUBROUTINE qcschema_update_fock
734 :
735 : ! **************************************************************************************************
736 : !> \brief Writes a qcschema object to an hdf5 file
737 : !> \param qcschema_env the qcschema environment to write to file
738 : !> \param filename ...
739 : ! **************************************************************************************************
740 28 : SUBROUTINE qcschema_to_hdf5(qcschema_env, filename)
741 : TYPE(qcschema_type), INTENT(IN) :: qcschema_env
742 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
743 : #ifndef __HDF5
744 : CPABORT("CP2K was compiled without the HDF5 library")
745 : MARK_USED(filename)
746 : MARK_USED(qcschema_env)
747 : #else
748 : INTEGER :: output_unit
749 : INTEGER(KIND=hdf5_id) :: file_id, group_id
750 : INTEGER(KIND=int_8) :: nresult
751 : TYPE(cp_logger_type), POINTER :: logger
752 :
753 4 : logger => cp_get_default_logger()
754 4 : output_unit = cp_logger_get_default_io_unit(logger)
755 :
756 : ! initialize HDF5 Fortran API
757 4 : CALL h5open()
758 :
759 : ! create qcschema hdf5 file
760 : ! filename = TRIM(logger%iter_info%project_name) // 'hdf5'
761 4 : CALL h5fcreate(TRIM(filename), file_id)
762 :
763 : ! !===========================================================================!
764 : ! *** Root group ***
765 : ! !===========================================================================!
766 : ! driver
767 4 : CALL h5awrite_fixlen_string(file_id, 'driver', TRIM(qcschema_env%driver))
768 : ! return result
769 4 : nresult = SIZE(qcschema_env%return_result)
770 4 : IF (SIZE(qcschema_env%return_result) == 1) THEN
771 4 : CALL h5awrite_double_scalar(file_id, 'return_result', qcschema_env%return_result(1))
772 : ELSE
773 0 : CALL h5awrite_double_simple(file_id, 'return_result', qcschema_env%return_result)
774 : END IF
775 : ! schema name
776 4 : CALL h5awrite_fixlen_string(file_id, 'schema_name', TRIM(qcschema_env%topology%schema_name))
777 : ! schema version
778 4 : CALL h5awrite_integer_scalar(file_id, 'schema_version', qcschema_env%topology%schema_version)
779 : ! success
780 4 : CALL h5awrite_boolean(file_id, 'success', qcschema_env%success)
781 :
782 : !========================================================================================!
783 : ! *** QCSchema provenance ***
784 : !========================================================================================!
785 : ! create the provenance group
786 4 : CALL h5gcreate(file_id, 'provenance', group_id)
787 : ! populate provenance
788 4 : CALL h5awrite_fixlen_string(group_id, 'creator', TRIM(qcschema_env%provenance%creator))
789 4 : CALL h5awrite_fixlen_string(group_id, 'routine', TRIM(qcschema_env%provenance%routine))
790 4 : CALL h5awrite_fixlen_string(group_id, 'version', TRIM(qcschema_env%provenance%version))
791 : ! close provenance group
792 4 : CALL h5gclose(group_id)
793 :
794 : !========================================================================================!
795 : ! *** QCSchema molecule ***
796 : !========================================================================================!
797 : ! create the molecule group
798 4 : CALL h5gcreate(file_id, 'molecule', group_id)
799 : ! populate molecule
800 4 : CALL h5awrite_double_simple(group_id, 'geometry', qcschema_env%topology%geometry)
801 4 : CALL h5awrite_integer_simple(group_id, 'atomic_numbers', qcschema_env%topology%atomic_numbers)
802 4 : CALL h5awrite_double_simple(group_id, 'masses', qcschema_env%topology%masses)
803 4 : CALL h5awrite_integer_scalar(group_id, 'molecular_charge', qcschema_env%topology%molecular_charge)
804 4 : CALL h5awrite_integer_scalar(group_id, 'molecular_multiplicity', qcschema_env%topology%molecular_multiplicity)
805 4 : CALL h5awrite_string_simple(group_id, 'symbols', qcschema_env%topology%symbols)
806 :
807 4 : CALL h5awrite_fixlen_string(group_id, 'schema_name', 'qcschema_molecule')
808 4 : CALL h5awrite_integer_scalar(group_id, 'schema_version', 2)
809 : ! close molecule group
810 4 : CALL h5gclose(group_id)
811 :
812 : !========================================================================================!
813 : ! *** QCSchema properties ***
814 : !========================================================================================!
815 : ! create the properties group
816 4 : CALL h5gcreate(file_id, 'properties', group_id)
817 : ! populate properties
818 4 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_natom', qcschema_env%properties%calcinfo_natom)
819 4 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbasis', qcschema_env%properties%calcinfo_nbasis)
820 4 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nmo', qcschema_env%properties%calcinfo_nmo)
821 4 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nalpha', qcschema_env%properties%calcinfo_nalpha)
822 4 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbeta', qcschema_env%properties%calcinfo_nbeta)
823 :
824 : ! CALL h5dwrite_double_simple(group_id, 'scf_dipole_moment', &
825 : ! qcschema_env%properties%scf_dipole_moment)
826 :
827 : ! energies, scf, mp2, ...
828 4 : CALL h5awrite_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy)
829 4 : CALL h5awrite_double_scalar(group_id, 'scf_total_energy', qcschema_env%properties%scf_total_energy)
830 : CALL h5awrite_double_scalar(group_id, 'nuclear_repulsion_energy', &
831 4 : qcschema_env%properties%nuclear_repulsion_energy)
832 :
833 4 : IF (qcschema_env%properties%scf_iterations /= 0) THEN
834 4 : CALL h5awrite_integer_scalar(group_id, 'scf_iterations', qcschema_env%properties%scf_iterations)
835 : END IF
836 :
837 4 : IF (qcschema_env%properties%scf_one_electron_energy /= 0.0_dp) THEN
838 : CALL h5awrite_double_scalar(group_id, 'scf_one_electron_energy', &
839 4 : qcschema_env%properties%scf_one_electron_energy)
840 : END IF
841 :
842 4 : IF (qcschema_env%properties%scf_two_electron_energy /= 0.0_dp) THEN
843 : CALL h5awrite_double_scalar(group_id, 'scf_two_electron_energy', &
844 4 : qcschema_env%properties%scf_two_electron_energy)
845 : END IF
846 :
847 4 : IF (qcschema_env%properties%scf_xc_energy /= 0.0_dp) THEN
848 : CALL h5awrite_double_scalar(group_id, 'scf_xc_energy', &
849 2 : qcschema_env%properties%scf_xc_energy)
850 : END IF
851 :
852 4 : IF (qcschema_env%properties%scf_dispersion_correction_energy /= 0.0_dp) THEN
853 : CALL h5awrite_double_scalar(group_id, 'scf_dispersion_correction_energy', &
854 0 : qcschema_env%properties%scf_dispersion_correction_energy)
855 : END IF
856 :
857 4 : IF (qcschema_env%properties%mp2) THEN
858 : CALL h5awrite_double_scalar(group_id, 'mp2_correlation_energy', &
859 0 : qcschema_env%properties%mp2_correlation_energy)
860 : END IF
861 :
862 : ! close properties group
863 4 : CALL h5gclose(group_id)
864 :
865 : !========================================================================================!
866 : ! *** QCSchema wavefunction ***
867 : !========================================================================================!
868 : ! create the wavefunction group
869 4 : CALL h5gcreate(file_id, 'wavefunction', group_id)
870 :
871 4 : CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name))
872 :
873 : CALL h5dwrite_double_simple(group_id, 'scf_orbitals_a', &
874 4 : qcschema_env%wavefunction%scf_orbitals_a)
875 :
876 : CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_a', &
877 4 : qcschema_env%wavefunction%scf_eigenvalues_a)
878 :
879 : CALL h5dwrite_double_simple(group_id, 'scf_occupations_a', &
880 4 : qcschema_env%wavefunction%scf_occupations_a)
881 :
882 : CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_a', &
883 4 : qcschema_env%wavefunction%scf_fock_mo_a)
884 :
885 : CALL h5dwrite_double_simple(group_id, 'scf_density_mo_a', &
886 4 : qcschema_env%wavefunction%scf_density_mo_a)
887 :
888 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_aa', &
889 4 : qcschema_env%wavefunction%scf_eri_mo_aa)
890 :
891 4 : IF (.NOT. qcschema_env%wavefunction%restricted) THEN
892 : CALL h5dwrite_double_simple(group_id, 'scf_orbitals_b', &
893 0 : qcschema_env%wavefunction%scf_orbitals_b)
894 :
895 : CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_b', &
896 0 : qcschema_env%wavefunction%scf_eigenvalues_b)
897 :
898 : CALL h5dwrite_double_simple(group_id, 'scf_occupations_b', &
899 0 : qcschema_env%wavefunction%scf_occupations_b)
900 :
901 : CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_b', &
902 0 : qcschema_env%wavefunction%scf_fock_mo_b)
903 :
904 : CALL h5dwrite_double_simple(group_id, 'scf_density_mo_b', &
905 0 : qcschema_env%wavefunction%scf_density_mo_b)
906 :
907 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_bb', &
908 0 : qcschema_env%wavefunction%scf_eri_mo_bb)
909 :
910 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_ba', &
911 0 : qcschema_env%wavefunction%scf_eri_mo_ab)
912 :
913 : END IF
914 :
915 : ! close wavefunction group
916 4 : CALL h5gclose(group_id)
917 :
918 : !========================================================================================!
919 : ! *** QCSchema model ***
920 : !========================================================================================!
921 : ! create the model group
922 4 : CALL h5gcreate(file_id, 'model', group_id)
923 4 : CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name))
924 4 : CALL h5awrite_fixlen_string(group_id, 'method', TRIM(qcschema_env%wavefunction%method))
925 : ! close model group
926 4 : CALL h5gclose(group_id)
927 :
928 : ! create the keywords group
929 4 : CALL h5gcreate(file_id, 'keywords', group_id)
930 : ! close keywords group
931 4 : CALL h5gclose(group_id)
932 :
933 4 : CALL h5fclose(file_id)
934 4 : CALL h5close()
935 : #endif
936 :
937 4 : END SUBROUTINE qcschema_to_hdf5
938 :
939 : #ifdef __HDF5
940 : ! **************************************************************************************************
941 : !> \brief Reads the electron density from a qcschema hdf5 file
942 : !> \param filename the path to the qcschema hdf5 file
943 : !> \param qcschema_env the qcschema environment onto which it writes the density
944 : ! **************************************************************************************************
945 0 : SUBROUTINE read_pmat_from_hdf5(filename, qcschema_env)
946 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
947 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
948 :
949 : INTEGER :: nmo
950 : INTEGER(KIND=hdf5_id) :: file_id, group_id
951 :
952 : ! initialize HDF5 Fortran API
953 0 : CALL h5open()
954 :
955 : ! open qcschema hdf5 file
956 0 : CALL h5fopen(TRIM(filename), file_id)
957 :
958 : ! open the wave function group
959 0 : CALL h5gopen(file_id, 'wavefunction', group_id)
960 :
961 : ! allocate the space for the array containing the density
962 0 : nmo = qcschema_env%properties%calcinfo_nmo
963 0 : IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN
964 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo))
965 : END IF
966 :
967 : ! read the alpha density
968 0 : CALL h5dread_double_simple(group_id, 'scf_density_mo_a', qcschema_env%wavefunction%scf_density_mo_a)
969 :
970 0 : IF (.NOT. qcschema_env%wavefunction%restricted) THEN
971 0 : IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN
972 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo))
973 : END IF
974 : ! read the beta density
975 0 : CALL h5dread_double_simple(group_id, 'scf_density_mo_b', qcschema_env%wavefunction%scf_density_mo_b)
976 : END IF
977 :
978 : ! close everything
979 0 : CALL h5gclose(group_id)
980 0 : CALL h5fclose(file_id)
981 0 : CALL h5close()
982 :
983 0 : END SUBROUTINE read_pmat_from_hdf5
984 :
985 : ! **************************************************************************************************
986 : !> \brief Reads the return energy from a qcschema hdf5 file
987 : !> \param filename the path to the qcschema hdf5 file
988 : !> \param qcschema_env the qcschema environment onto which it writes the energy
989 : ! **************************************************************************************************
990 0 : SUBROUTINE read_return_energy_from_hdf5(filename, qcschema_env)
991 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
992 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
993 :
994 : INTEGER(KIND=hdf5_id) :: file_id, group_id
995 :
996 : ! initialize HDF5 Fortran API
997 0 : CALL h5open()
998 :
999 : ! open qcschema hdf5 file
1000 0 : CALL h5fopen(TRIM(filename), file_id)
1001 :
1002 : ! open the properties group
1003 0 : CALL h5gopen(file_id, 'properties', group_id)
1004 :
1005 : ! read the return energy
1006 0 : CALL h5aread_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy)
1007 :
1008 : ! close everything
1009 0 : CALL h5gclose(group_id)
1010 0 : CALL h5fclose(file_id)
1011 0 : CALL h5close()
1012 :
1013 0 : END SUBROUTINE read_return_energy_from_hdf5
1014 :
1015 : ! **************************************************************************************************
1016 : !> \brief Reads the active space energy from a qcschema file and stores it in active_space_env
1017 : !> \param active_space_env ...
1018 : !> \param qcschema_env ...
1019 : !> \author Stefano Battaglia
1020 : ! **************************************************************************************************
1021 0 : SUBROUTINE read_active_energy_from_hdf5(active_space_env, qcschema_env)
1022 : TYPE(active_space_type), POINTER :: active_space_env
1023 : TYPE(qcschema_type) :: qcschema_env
1024 :
1025 : CHARACTER(LEN=default_path_length) :: qcschema_filename
1026 :
1027 : ! File name
1028 0 : qcschema_filename = active_space_env%qcschema_filename
1029 : ! read active space energy
1030 0 : CALL read_return_energy_from_hdf5(qcschema_filename, qcschema_env)
1031 :
1032 0 : active_space_env%energy_active = qcschema_env%properties%return_energy
1033 0 : active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
1034 :
1035 0 : END SUBROUTINE read_active_energy_from_hdf5
1036 : #endif
1037 :
1038 0 : END MODULE qcschema
|
