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 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 0 : 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 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
290 : TYPE(qs_energy_type), POINTER :: energy
291 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
292 0 : 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 0 : CALL timeset(routineN, handle)
298 :
299 0 : logger => cp_get_default_logger()
300 0 : output_unit = cp_logger_get_default_io_unit(logger)
301 :
302 : ! reset everything
303 0 : CALL qcschema_env_release(qcschema_env)
304 :
305 : ! collect environment info
306 0 : 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 0 : 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 0 : 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 0 : qcschema_env%provenance%creator = 'CP2K'
327 0 : qcschema_env%provenance%version = cp2k_version
328 0 : qcschema_env%provenance%routine = routineN
329 :
330 : !========================================================================================!
331 : ! *** QCSchema topology ***
332 : !========================================================================================!
333 :
334 0 : qcschema_env%topology%schema_name = 'qcschema'
335 0 : qcschema_env%topology%schema_version = 3
336 :
337 0 : natoms = SIZE(particle_set)
338 :
339 0 : ALLOCATE (qcschema_env%topology%geometry(3*natoms))
340 0 : ALLOCATE (qcschema_env%topology%symbols(natoms))
341 0 : ALLOCATE (qcschema_env%topology%atomic_numbers(natoms))
342 0 : ALLOCATE (qcschema_env%topology%masses(natoms))
343 :
344 0 : DO iatom = 1, natoms
345 : ! set the geometry as a flat array
346 0 : qcschema_env%topology%geometry((iatom - 1)*3 + 1:(iatom)*3) = particle_set(iatom)%r(1:3)
347 :
348 : ! set the atomic symbols
349 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, element_symbol=atomic_symbol)
350 0 : qcschema_env%topology%symbols(iatom) = atomic_symbol
351 :
352 : ! set the atomic numbers and masses
353 0 : CALL get_ptable_info(atomic_symbol, number=atomic_number, amass=mass)
354 0 : qcschema_env%topology%atomic_numbers(iatom) = atomic_number
355 0 : qcschema_env%topology%masses(iatom) = mass
356 : END DO
357 :
358 0 : qcschema_env%topology%molecular_charge = dft_control%charge
359 0 : qcschema_env%topology%molecular_multiplicity = dft_control%multiplicity
360 :
361 : !========================================================================================!
362 : ! *** QCSchema properties ***
363 : !========================================================================================!
364 :
365 0 : nspins = active_space_env%nspins
366 :
367 0 : nao = active_space_env%mos_active(1)%nao
368 0 : nmo = active_space_env%nmo_active
369 0 : nel = active_space_env%nelec_active
370 :
371 0 : IF (nspins == 1) THEN
372 0 : nalpha = active_space_env%nelec_active/2
373 0 : 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 0 : qcschema_env%properties%calcinfo_natom = natoms
380 0 : qcschema_env%properties%calcinfo_nbasis = nao
381 0 : qcschema_env%properties%calcinfo_nmo = nmo
382 0 : qcschema_env%properties%calcinfo_nalpha = nalpha
383 0 : qcschema_env%properties%calcinfo_nbeta = nbeta
384 :
385 : ! energy results
386 0 : qcschema_env%properties%return_energy = energy%total
387 0 : qcschema_env%properties%scf_total_energy = energy%total
388 : ! here we abuse the nuclear repulsion energy to store the inactive energy
389 0 : qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive
390 : ! SCF info
391 0 : qcschema_env%properties%scf_iterations = scf_env%iter_count
392 : ! one-electron energy is the sum of all core terms
393 0 : one_el_en = energy%core_overlap + energy%core_self + energy%core
394 0 : 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 0 : two_el_en = energy%hartree + energy%ex + energy%hartree_1c
397 0 : qcschema_env%properties%scf_one_electron_energy = two_el_en
398 : ! xc energy
399 : qcschema_env%properties%scf_xc_energy = &
400 0 : energy%exc + energy%exc_aux_fit + energy%exc1 + energy%exc1_aux_fit
401 : ! dispersion energy
402 0 : dispersion = energy%dispersion + energy%gcp
403 0 : 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 0 : IF (dft_control%smear) CPABORT('WARNING: smearing not supported in QCSchema')
407 0 : IF (dft_control%dft_plus_u) CPABORT('WARNING: DFT+U not supported in QCSchema')
408 0 : IF (dft_control%do_sccs) CPABORT('WARNING: SCCS not supported in QCSchema')
409 0 : IF (qs_env%qmmm) CPABORT('WARNING: QM/MM not supported in QCSchema')
410 0 : IF (dft_control%qs_control%mulliken_restraint) &
411 0 : CPABORT('WARNING: Mulliken restrains not supported in QCSchema')
412 0 : IF (dft_control%qs_control%semi_empirical) &
413 0 : CPABORT('WARNING: semi_empirical methods not supported in QCSchema')
414 0 : IF (dft_control%qs_control%dftb) CPABORT('WARNING: DFTB not supported in QCSchema')
415 0 : IF (dft_control%qs_control%xtb) CPABORT('WARNING: xTB not supported in QCSchema')
416 :
417 : ! MP2 info
418 0 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
419 0 : qcschema_env%properties%mp2 = .TRUE.
420 : ! this info is computed on the fly, but not stored!
421 : ! qcschema_env%properties%mp2_same_spin_correlation_energy
422 : ! qcschema_env%properties%mp2_opposite_spin_correlation_energy
423 :
424 0 : qcschema_env%properties%mp2_correlation_energy = energy%mp2
425 0 : qcschema_env%properties%mp2_total_energy = energy%total
426 :
427 : ! update the scf energy
428 0 : qcschema_env%properties%scf_total_energy = energy%total - energy%mp2
429 : END IF
430 :
431 : !========================================================================================!
432 : ! *** QCSchema wavefunction ***
433 : !========================================================================================!
434 :
435 0 : IF (nspins == 1) THEN
436 0 : qcschema_env%wavefunction%restricted = .TRUE.
437 : ELSE
438 0 : qcschema_env%wavefunction%restricted = .FALSE.
439 : END IF
440 :
441 : ! alpha MO energies
442 0 : ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a(nmo))
443 0 : DO i = 1, nmo
444 0 : i_glb = active_space_env%active_orbitals(i, 1)
445 : qcschema_env%wavefunction%scf_eigenvalues_a(i) = &
446 0 : active_space_env%mos_active(1)%eigenvalues(i_glb)
447 : END DO
448 :
449 : ! alpha MO occupations
450 0 : ALLOCATE (qcschema_env%wavefunction%scf_occupations_a(nmo))
451 0 : DO i = 1, nmo
452 0 : i_glb = active_space_env%active_orbitals(i, 1)
453 : qcschema_env%wavefunction%scf_occupations_a(i) = &
454 0 : active_space_env%mos_active(1)%occupation_numbers(i_glb)
455 : END DO
456 :
457 : ! alpha Fock matrix
458 0 : ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a(nmo*nmo))
459 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), &
460 : qcschema_env%wavefunction%scf_fock_mo_a, &
461 : active_space_env%active_orbitals(:, 1), &
462 0 : active_space_env%active_orbitals(:, 1))
463 :
464 : ! alpha density matrix
465 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo))
466 : CALL subspace_matrix_to_array(active_space_env%p_active(1), &
467 : qcschema_env%wavefunction%scf_density_mo_a, &
468 : active_space_env%active_orbitals(:, 1), &
469 0 : active_space_env%active_orbitals(:, 1))
470 :
471 : ! alpha MOs coefficients
472 0 : ALLOCATE (qcschema_env%wavefunction%scf_orbitals_a(nao*nmo))
473 : CALL subspace_matrix_to_array(active_space_env%mos_active(1)%mo_coeff, &
474 : qcschema_env%wavefunction%scf_orbitals_a, &
475 0 : (/(i, i=1, nao)/), active_space_env%active_orbitals(:, 1))
476 :
477 0 : IF (nspins == 2) THEN
478 : ! beta MO energies
479 0 : ALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b(nmo))
480 0 : DO i = 1, nmo
481 0 : i_glb = active_space_env%active_orbitals(i, 2)
482 : qcschema_env%wavefunction%scf_eigenvalues_b(i) = &
483 0 : active_space_env%mos_active(2)%eigenvalues(i_glb)
484 : END DO
485 :
486 : ! beta MO occupations
487 0 : ALLOCATE (qcschema_env%wavefunction%scf_occupations_b(nmo))
488 0 : DO i = 1, nmo
489 0 : i_glb = active_space_env%active_orbitals(i, 2)
490 : qcschema_env%wavefunction%scf_occupations_b(i) = &
491 0 : active_space_env%mos_active(2)%occupation_numbers(i_glb)
492 : END DO
493 :
494 : ! beta Fock matrix
495 0 : ALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b(nmo*nmo))
496 : CALL subspace_matrix_to_array(active_space_env%fock_sub(2), &
497 : qcschema_env%wavefunction%scf_fock_mo_b, &
498 : active_space_env%active_orbitals(:, 2), &
499 0 : active_space_env%active_orbitals(:, 2))
500 :
501 : ! beta density matrix
502 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo))
503 : CALL subspace_matrix_to_array(active_space_env%p_active(2), &
504 : qcschema_env%wavefunction%scf_density_mo_b, &
505 : active_space_env%active_orbitals(:, 2), &
506 0 : active_space_env%active_orbitals(:, 2))
507 :
508 : ! beta MOs coefficients
509 0 : ALLOCATE (qcschema_env%wavefunction%scf_orbitals_b(nao*nmo))
510 : CALL subspace_matrix_to_array(active_space_env%mos_active(2)%mo_coeff, &
511 : qcschema_env%wavefunction%scf_orbitals_b, &
512 0 : (/(i, i=1, nao)/), active_space_env%active_orbitals(:, 2))
513 : END IF
514 :
515 : ! get the alpha-alpha eri
516 0 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa(nmo**4))
517 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_aa, &
518 0 : active_space_env%active_orbitals, 1, 1)
519 :
520 0 : IF (nspins == 2) THEN
521 : ! get the alpha-beta eri
522 0 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab(nmo**4))
523 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_ab, &
524 0 : active_space_env%active_orbitals, 1, 2)
525 :
526 : ! get the beta-beta eri
527 0 : ALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb(nmo**4))
528 : CALL eri_to_array(active_space_env%eri, qcschema_env%wavefunction%scf_eri_mo_bb, &
529 0 : active_space_env%active_orbitals, 2, 2)
530 : END IF
531 :
532 : !========================================================================================!
533 : ! *** QCSchema model ***
534 : !========================================================================================!
535 :
536 0 : DO iatom = 1, natoms
537 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
538 0 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_set)
539 :
540 0 : basis_set_name = basis_set%name
541 :
542 : ! make sure that we do not run a mixed basis set
543 0 : IF (iatom > 1) THEN
544 0 : CPASSERT(basis_set_name == basis_set%name)
545 : END IF
546 : END DO
547 0 : qcschema_env%wavefunction%basis_set%name = basis_set_name
548 :
549 : ! figure out which method was used for the calculation
550 0 : IF (dft_control%uks) THEN
551 0 : method = 'U'
552 0 : ELSE IF (dft_control%roks) THEN
553 0 : method = 'RO'
554 : ELSE
555 0 : method = 'R'
556 : END IF
557 :
558 0 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
559 0 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
560 :
561 0 : IF (do_hfx) THEN
562 0 : method = TRIM(method)//'HF'
563 0 : ELSE IF (qcschema_env%properties%mp2) THEN
564 0 : method = TRIM(method)//'MP2'
565 : ELSE
566 0 : method = TRIM(method)//'KS'
567 : END IF
568 :
569 0 : qcschema_env%wavefunction%method = TRIM(method)
570 :
571 : !========================================================================================!
572 : ! *** QCSchema root ***
573 : !========================================================================================!
574 :
575 : ! driver
576 0 : IF (ASSOCIATED(force)) THEN
577 0 : qcschema_env%driver = 'gradient'
578 : ELSE
579 0 : qcschema_env%driver = 'energy'
580 : END IF
581 :
582 : ! success
583 : ! TODO: how to check if the calculation was succesful?
584 0 : qcschema_env%success = .TRUE.
585 :
586 : ! return result
587 : IF (qcschema_env%success) THEN
588 0 : IF (qcschema_env%driver == 'energy') THEN
589 0 : ALLOCATE (qcschema_env%return_result(1))
590 0 : qcschema_env%return_result(1) = energy%total
591 : ELSE
592 0 : ALLOCATE (qcschema_env%return_result(3*SIZE(particle_set)))
593 : ! TODO: populate with forces!!
594 0 : qcschema_env%return_result = 0.0_dp
595 : END IF
596 : ELSE
597 : CPABORT("The calculation to build the AS is unsuccessful")
598 : END IF
599 :
600 0 : CALL timestop(handle)
601 :
602 0 : END SUBROUTINE qcschema_env_create
603 :
604 : ! **************************************************************************************************
605 : !> \brief Releases the allocated memory of a qcschema environment
606 : !> \param qcschema_env the qcschema environment to release
607 : ! **************************************************************************************************
608 0 : SUBROUTINE qcschema_env_release(qcschema_env)
609 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
610 :
611 0 : IF (ALLOCATED(qcschema_env%return_result)) THEN
612 0 : DEALLOCATE (qcschema_env%return_result)
613 : END IF
614 :
615 0 : IF (ALLOCATED(qcschema_env%topology%atomic_numbers)) THEN
616 0 : DEALLOCATE (qcschema_env%topology%atomic_numbers)
617 : END IF
618 :
619 0 : IF (ALLOCATED(qcschema_env%topology%masses)) THEN
620 0 : DEALLOCATE (qcschema_env%topology%masses)
621 : END IF
622 :
623 0 : IF (ALLOCATED(qcschema_env%topology%geometry)) THEN
624 0 : DEALLOCATE (qcschema_env%topology%geometry)
625 : END IF
626 :
627 0 : IF (ALLOCATED(qcschema_env%topology%symbols)) THEN
628 0 : DEALLOCATE (qcschema_env%topology%symbols)
629 : END IF
630 :
631 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN
632 0 : DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_a)
633 : END IF
634 :
635 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN
636 0 : DEALLOCATE (qcschema_env%wavefunction%scf_density_mo_b)
637 : END IF
638 :
639 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_a)) THEN
640 0 : DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_a)
641 : END IF
642 :
643 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_fock_mo_b)) THEN
644 0 : DEALLOCATE (qcschema_env%wavefunction%scf_fock_mo_b)
645 : END IF
646 :
647 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_a)) THEN
648 0 : DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_a)
649 : END IF
650 :
651 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_orbitals_b)) THEN
652 0 : DEALLOCATE (qcschema_env%wavefunction%scf_orbitals_b)
653 : END IF
654 :
655 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_a)) THEN
656 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_a)
657 : END IF
658 :
659 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eigenvalues_b)) THEN
660 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eigenvalues_b)
661 : END IF
662 :
663 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_a)) THEN
664 0 : DEALLOCATE (qcschema_env%wavefunction%scf_occupations_a)
665 : END IF
666 :
667 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_occupations_b)) THEN
668 0 : DEALLOCATE (qcschema_env%wavefunction%scf_occupations_b)
669 : END IF
670 :
671 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri)) THEN
672 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri)
673 : END IF
674 :
675 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_aa)) THEN
676 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_aa)
677 : END IF
678 :
679 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_bb)) THEN
680 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_bb)
681 : END IF
682 :
683 0 : IF (ALLOCATED(qcschema_env%wavefunction%scf_eri_mo_ab)) THEN
684 0 : DEALLOCATE (qcschema_env%wavefunction%scf_eri_mo_ab)
685 : END IF
686 :
687 0 : IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_a)) THEN
688 0 : DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_a)
689 : END IF
690 :
691 0 : IF (ALLOCATED(qcschema_env%wavefunction%localized_orbitals_b)) THEN
692 0 : DEALLOCATE (qcschema_env%wavefunction%localized_orbitals_b)
693 : END IF
694 :
695 0 : IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_a)) THEN
696 0 : DEALLOCATE (qcschema_env%wavefunction%localized_fock_a)
697 : END IF
698 :
699 0 : IF (ALLOCATED(qcschema_env%wavefunction%localized_fock_b)) THEN
700 0 : DEALLOCATE (qcschema_env%wavefunction%localized_fock_b)
701 : END IF
702 :
703 0 : END SUBROUTINE qcschema_env_release
704 :
705 : ! **************************************************************************************************
706 : !> \brief Updates the Fock matrix and the inactive energy in a qcschema object
707 : !> \param qcschema_env the qcschema environment
708 : !> \param active_space_env the active space environment with the updated data
709 : ! **************************************************************************************************
710 0 : SUBROUTINE qcschema_update_fock(qcschema_env, active_space_env)
711 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
712 : TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env
713 :
714 : ! alpha Fock matrix
715 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), &
716 : qcschema_env%wavefunction%scf_fock_mo_a, &
717 : active_space_env%active_orbitals(:, 1), &
718 0 : active_space_env%active_orbitals(:, 1))
719 :
720 : ! beta Fock matrix
721 0 : IF (active_space_env%nspins == 2) THEN
722 : CALL subspace_matrix_to_array(active_space_env%fock_sub(2), &
723 : qcschema_env%wavefunction%scf_fock_mo_b, &
724 : active_space_env%active_orbitals(:, 2), &
725 0 : active_space_env%active_orbitals(:, 2))
726 : END IF
727 :
728 : ! update inactive energy
729 0 : qcschema_env%properties%nuclear_repulsion_energy = active_space_env%energy_inactive
730 :
731 0 : END SUBROUTINE qcschema_update_fock
732 :
733 : ! **************************************************************************************************
734 : !> \brief Writes a qcschema object to an hdf5 file
735 : !> \param qcschema_env the qcschema environment to write to file
736 : !> \param filename ...
737 : ! **************************************************************************************************
738 0 : SUBROUTINE qcschema_to_hdf5(qcschema_env, filename)
739 : TYPE(qcschema_type), INTENT(IN) :: qcschema_env
740 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
741 : #ifndef __HDF5
742 : CPABORT("CP2K was compiled without the HDF5 library")
743 : MARK_USED(filename)
744 : MARK_USED(qcschema_env)
745 : #else
746 : INTEGER :: output_unit
747 : INTEGER(KIND=hdf5_id) :: file_id, group_id
748 : INTEGER(KIND=int_8) :: nresult
749 : TYPE(cp_logger_type), POINTER :: logger
750 :
751 0 : logger => cp_get_default_logger()
752 0 : output_unit = cp_logger_get_default_io_unit(logger)
753 :
754 : ! initialize HDF5 Fortran API
755 0 : CALL h5open()
756 :
757 : ! create qcschema hdf5 file
758 : ! filename = TRIM(logger%iter_info%project_name) // 'hdf5'
759 0 : CALL h5fcreate(TRIM(filename), file_id)
760 :
761 : ! !===========================================================================!
762 : ! *** Root group ***
763 : ! !===========================================================================!
764 : ! driver
765 0 : CALL h5awrite_fixlen_string(file_id, 'driver', TRIM(qcschema_env%driver))
766 : ! return result
767 0 : nresult = SIZE(qcschema_env%return_result)
768 0 : IF (SIZE(qcschema_env%return_result) == 1) THEN
769 0 : CALL h5awrite_double_scalar(file_id, 'return_result', qcschema_env%return_result(1))
770 : ELSE
771 0 : CALL h5awrite_double_simple(file_id, 'return_result', qcschema_env%return_result)
772 : END IF
773 : ! schema name
774 0 : CALL h5awrite_fixlen_string(file_id, 'schema_name', TRIM(qcschema_env%topology%schema_name))
775 : ! schema version
776 0 : CALL h5awrite_integer_scalar(file_id, 'schema_version', qcschema_env%topology%schema_version)
777 : ! success
778 0 : CALL h5awrite_boolean(file_id, 'success', qcschema_env%success)
779 :
780 : !========================================================================================!
781 : ! *** QCSchema provenance ***
782 : !========================================================================================!
783 : ! create the provenance group
784 0 : CALL h5gcreate(file_id, 'provenance', group_id)
785 : ! populate provenance
786 0 : CALL h5awrite_fixlen_string(group_id, 'creator', TRIM(qcschema_env%provenance%creator))
787 0 : CALL h5awrite_fixlen_string(group_id, 'routine', TRIM(qcschema_env%provenance%routine))
788 0 : CALL h5awrite_fixlen_string(group_id, 'version', TRIM(qcschema_env%provenance%version))
789 : ! close provenance group
790 0 : CALL h5gclose(group_id)
791 :
792 : !========================================================================================!
793 : ! *** QCSchema molecule ***
794 : !========================================================================================!
795 : ! create the molecule group
796 0 : CALL h5gcreate(file_id, 'molecule', group_id)
797 : ! populate molecule
798 0 : CALL h5awrite_double_simple(group_id, 'geometry', qcschema_env%topology%geometry)
799 0 : CALL h5awrite_integer_simple(group_id, 'atomic_numbers', qcschema_env%topology%atomic_numbers)
800 0 : CALL h5awrite_double_simple(group_id, 'masses', qcschema_env%topology%masses)
801 0 : CALL h5awrite_integer_scalar(group_id, 'molecular_charge', qcschema_env%topology%molecular_charge)
802 0 : CALL h5awrite_integer_scalar(group_id, 'molecular_multiplicity', qcschema_env%topology%molecular_multiplicity)
803 0 : CALL h5awrite_string_simple(group_id, 'symbols', qcschema_env%topology%symbols)
804 :
805 0 : CALL h5awrite_fixlen_string(group_id, 'schema_name', 'qcschema_molecule')
806 0 : CALL h5awrite_integer_scalar(group_id, 'schema_version', 2)
807 : ! close molecule group
808 0 : CALL h5gclose(group_id)
809 :
810 : !========================================================================================!
811 : ! *** QCSchema properties ***
812 : !========================================================================================!
813 : ! create the properties group
814 0 : CALL h5gcreate(file_id, 'properties', group_id)
815 : ! populate properties
816 0 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_natom', qcschema_env%properties%calcinfo_natom)
817 0 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbasis', qcschema_env%properties%calcinfo_nbasis)
818 0 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nmo', qcschema_env%properties%calcinfo_nmo)
819 0 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nalpha', qcschema_env%properties%calcinfo_nalpha)
820 0 : CALL h5awrite_integer_scalar(group_id, 'calcinfo_nbeta', qcschema_env%properties%calcinfo_nbeta)
821 :
822 : ! CALL h5dwrite_double_simple(group_id, 'scf_dipole_moment', &
823 : ! qcschema_env%properties%scf_dipole_moment)
824 :
825 : ! energies, scf, mp2, ...
826 0 : CALL h5awrite_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy)
827 0 : CALL h5awrite_double_scalar(group_id, 'scf_total_energy', qcschema_env%properties%scf_total_energy)
828 : CALL h5awrite_double_scalar(group_id, 'nuclear_repulsion_energy', &
829 0 : qcschema_env%properties%nuclear_repulsion_energy)
830 :
831 0 : IF (qcschema_env%properties%scf_iterations /= 0) THEN
832 0 : CALL h5awrite_integer_scalar(group_id, 'scf_iterations', qcschema_env%properties%scf_iterations)
833 : END IF
834 :
835 0 : IF (qcschema_env%properties%scf_one_electron_energy /= 0.0_dp) THEN
836 : CALL h5awrite_double_scalar(group_id, 'scf_one_electron_energy', &
837 0 : qcschema_env%properties%scf_one_electron_energy)
838 : END IF
839 :
840 0 : IF (qcschema_env%properties%scf_two_electron_energy /= 0.0_dp) THEN
841 : CALL h5awrite_double_scalar(group_id, 'scf_two_electron_energy', &
842 0 : qcschema_env%properties%scf_two_electron_energy)
843 : END IF
844 :
845 0 : IF (qcschema_env%properties%scf_xc_energy /= 0.0_dp) THEN
846 : CALL h5awrite_double_scalar(group_id, 'scf_xc_energy', &
847 0 : qcschema_env%properties%scf_xc_energy)
848 : END IF
849 :
850 0 : IF (qcschema_env%properties%scf_dispersion_correction_energy /= 0.0_dp) THEN
851 : CALL h5awrite_double_scalar(group_id, 'scf_dispersion_correction_energy', &
852 0 : qcschema_env%properties%scf_dispersion_correction_energy)
853 : END IF
854 :
855 0 : IF (qcschema_env%properties%mp2) THEN
856 : CALL h5awrite_double_scalar(group_id, 'mp2_correlation_energy', &
857 0 : qcschema_env%properties%mp2_correlation_energy)
858 : END IF
859 :
860 : ! close properties group
861 0 : CALL h5gclose(group_id)
862 :
863 : !========================================================================================!
864 : ! *** QCSchema wavefunction ***
865 : !========================================================================================!
866 : ! create the wavefunction group
867 0 : CALL h5gcreate(file_id, 'wavefunction', group_id)
868 :
869 0 : CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name))
870 :
871 : CALL h5dwrite_double_simple(group_id, 'scf_orbitals_a', &
872 0 : qcschema_env%wavefunction%scf_orbitals_a)
873 :
874 : CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_a', &
875 0 : qcschema_env%wavefunction%scf_eigenvalues_a)
876 :
877 : CALL h5dwrite_double_simple(group_id, 'scf_occupations_a', &
878 0 : qcschema_env%wavefunction%scf_occupations_a)
879 :
880 : CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_a', &
881 0 : qcschema_env%wavefunction%scf_fock_mo_a)
882 :
883 : CALL h5dwrite_double_simple(group_id, 'scf_density_mo_a', &
884 0 : qcschema_env%wavefunction%scf_density_mo_a)
885 :
886 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_aa', &
887 0 : qcschema_env%wavefunction%scf_eri_mo_aa)
888 :
889 0 : IF (.NOT. qcschema_env%wavefunction%restricted) THEN
890 : CALL h5dwrite_double_simple(group_id, 'scf_orbitals_b', &
891 0 : qcschema_env%wavefunction%scf_orbitals_b)
892 :
893 : CALL h5dwrite_double_simple(group_id, 'scf_eigenvalues_b', &
894 0 : qcschema_env%wavefunction%scf_eigenvalues_b)
895 :
896 : CALL h5dwrite_double_simple(group_id, 'scf_occupations_b', &
897 0 : qcschema_env%wavefunction%scf_occupations_b)
898 :
899 : CALL h5dwrite_double_simple(group_id, 'scf_fock_mo_b', &
900 0 : qcschema_env%wavefunction%scf_fock_mo_b)
901 :
902 : CALL h5dwrite_double_simple(group_id, 'scf_density_mo_b', &
903 0 : qcschema_env%wavefunction%scf_density_mo_b)
904 :
905 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_bb', &
906 0 : qcschema_env%wavefunction%scf_eri_mo_bb)
907 :
908 : CALL h5dwrite_double_simple(group_id, 'scf_eri_mo_ba', &
909 0 : qcschema_env%wavefunction%scf_eri_mo_ab)
910 :
911 : END IF
912 :
913 : ! close wavefunction group
914 0 : CALL h5gclose(group_id)
915 :
916 : !========================================================================================!
917 : ! *** QCSchema model ***
918 : !========================================================================================!
919 : ! create the model group
920 0 : CALL h5gcreate(file_id, 'model', group_id)
921 0 : CALL h5awrite_fixlen_string(group_id, 'basis', TRIM(qcschema_env%wavefunction%basis_set%name))
922 0 : CALL h5awrite_fixlen_string(group_id, 'method', TRIM(qcschema_env%wavefunction%method))
923 : ! close model group
924 0 : CALL h5gclose(group_id)
925 :
926 : ! create the keywords group
927 0 : CALL h5gcreate(file_id, 'keywords', group_id)
928 : ! close keywords group
929 0 : CALL h5gclose(group_id)
930 :
931 0 : CALL h5fclose(file_id)
932 0 : CALL h5close()
933 : #endif
934 :
935 0 : END SUBROUTINE qcschema_to_hdf5
936 :
937 : #ifdef __HDF5
938 : ! **************************************************************************************************
939 : !> \brief Reads the electron density from a qcschema hdf5 file
940 : !> \param filename the path to the qcschema hdf5 file
941 : !> \param qcschema_env the qcschema environment onto which it writes the density
942 : ! **************************************************************************************************
943 0 : SUBROUTINE read_pmat_from_hdf5(filename, qcschema_env)
944 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
945 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
946 :
947 : INTEGER :: nmo
948 : INTEGER(KIND=hdf5_id) :: file_id, group_id
949 :
950 : ! initialize HDF5 Fortran API
951 0 : CALL h5open()
952 :
953 : ! open qcschema hdf5 file
954 0 : CALL h5fopen(TRIM(filename), file_id)
955 :
956 : ! open the wave function group
957 0 : CALL h5gopen(file_id, 'wavefunction', group_id)
958 :
959 : ! allocate the space for the array containing the density
960 0 : nmo = qcschema_env%properties%calcinfo_nmo
961 0 : IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_a)) THEN
962 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_a(nmo*nmo))
963 : END IF
964 :
965 : ! read the alpha density
966 0 : CALL h5dread_double_simple(group_id, 'scf_density_mo_a', qcschema_env%wavefunction%scf_density_mo_a)
967 :
968 0 : IF (.NOT. qcschema_env%wavefunction%restricted) THEN
969 0 : IF (.NOT. ALLOCATED(qcschema_env%wavefunction%scf_density_mo_b)) THEN
970 0 : ALLOCATE (qcschema_env%wavefunction%scf_density_mo_b(nmo*nmo))
971 : END IF
972 : ! read the beta density
973 0 : CALL h5dread_double_simple(group_id, 'scf_density_mo_b', qcschema_env%wavefunction%scf_density_mo_b)
974 : END IF
975 :
976 : ! close everything
977 0 : CALL h5gclose(group_id)
978 0 : CALL h5fclose(file_id)
979 0 : CALL h5close()
980 :
981 0 : END SUBROUTINE read_pmat_from_hdf5
982 :
983 : ! **************************************************************************************************
984 : !> \brief Reads the return energy from a qcschema hdf5 file
985 : !> \param filename the path to the qcschema hdf5 file
986 : !> \param qcschema_env the qcschema environment onto which it writes the energy
987 : ! **************************************************************************************************
988 0 : SUBROUTINE read_return_energy_from_hdf5(filename, qcschema_env)
989 : CHARACTER(LEN=default_path_length), INTENT(IN) :: filename
990 : TYPE(qcschema_type), INTENT(INOUT) :: qcschema_env
991 :
992 : INTEGER(KIND=hdf5_id) :: file_id, group_id
993 :
994 : ! initialize HDF5 Fortran API
995 0 : CALL h5open()
996 :
997 : ! open qcschema hdf5 file
998 0 : CALL h5fopen(TRIM(filename), file_id)
999 :
1000 : ! open the properties group
1001 0 : CALL h5gopen(file_id, 'properties', group_id)
1002 :
1003 : ! read the return energy
1004 0 : CALL h5aread_double_scalar(group_id, 'return_energy', qcschema_env%properties%return_energy)
1005 :
1006 : ! close everything
1007 0 : CALL h5gclose(group_id)
1008 0 : CALL h5fclose(file_id)
1009 0 : CALL h5close()
1010 :
1011 0 : END SUBROUTINE read_return_energy_from_hdf5
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief Reads the active space energy from a qcschema file and stores it in active_space_env
1015 : !> \param active_space_env ...
1016 : !> \param qcschema_env ...
1017 : !> \author Stefano Battaglia
1018 : ! **************************************************************************************************
1019 0 : SUBROUTINE read_active_energy_from_hdf5(active_space_env, qcschema_env)
1020 : TYPE(active_space_type), POINTER :: active_space_env
1021 : TYPE(qcschema_type) :: qcschema_env
1022 :
1023 : CHARACTER(LEN=default_path_length) :: qcschema_filename
1024 :
1025 : ! File name
1026 0 : qcschema_filename = active_space_env%qcschema_filename
1027 : ! read active space energy
1028 0 : CALL read_return_energy_from_hdf5(qcschema_filename, qcschema_env)
1029 :
1030 0 : active_space_env%energy_active = qcschema_env%properties%return_energy
1031 0 : active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
1032 :
1033 0 : END SUBROUTINE read_active_energy_from_hdf5
1034 : #endif
1035 :
1036 0 : END MODULE qcschema
|
