(git:852b9a4)
Loading...
Searching...
No Matches
qs_active_space_methods.F
Go to the documentation of this file.
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 Determine active space Hamiltonian
10!> \par History
11!> 04.2016 created [JGH]
12!> \author JGH
13! **************************************************************************************************
14MODULE qs_active_space_methods
15 USE admm_types, ONLY: admm_type, &
16 get_admm_env, &
17 admm_env_release
18 USE atomic_kind_types, ONLY: atomic_kind_type
19 USE basis_set_types, ONLY: allocate_sto_basis_set, &
20 create_gto_from_sto_basis, &
21 deallocate_sto_basis_set, &
22 gto_basis_set_type, &
23 init_orb_basis_set, &
24 set_sto_basis_set, &
25 srules, &
26 sto_basis_set_type
27 USE cell_types, ONLY: cell_type, use_perd_none, use_perd_xyz
28 USE cell_methods, ONLY: init_cell, set_cell_param, write_cell_low
29 USE cp_blacs_env, ONLY: cp_blacs_env_type, cp_blacs_env_create, cp_blacs_env_release, blacs_grid_square
30 USE cp_control_types, ONLY: dft_control_type, qs_control_type
31 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t, &
32 cp_dbcsr_sm_fm_multiply, &
33 dbcsr_allocate_matrix_set, &
34 cp_dbcsr_m_by_n_from_template, copy_dbcsr_to_fm
35 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
36 USE cp_files, ONLY: close_file, &
37 file_exists, &
38 open_file
39 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
40 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
41 cp_fm_struct_release, &
42 cp_fm_struct_type
43 USE cp_fm_types, ONLY: &
44 cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_init_random, cp_fm_release, &
45 cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm, cp_fm_type, cp_fm_write_formatted
46 USE cp_log_handling, ONLY: cp_get_default_logger, &
47 cp_logger_get_default_io_unit, &
48 cp_logger_type
49 USE cp_output_handling, ONLY: &
50 cp_p_file, cp_print_key_finished_output, cp_print_key_should_output, cp_print_key_unit_nr, &
51 debug_print_level, high_print_level, low_print_level, medium_print_level, &
52 silent_print_level
53 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
54 USE cp_dbcsr_api, ONLY: &
55 dbcsr_copy, dbcsr_csr_create, dbcsr_csr_type, dbcsr_p_type, dbcsr_type, dbcsr_release, &
56 dbcsr_type_no_symmetry, dbcsr_create, dbcsr_set, dbcsr_multiply, dbcsr_iterator_next_block, &
57 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_blocks_left, &
58 dbcsr_iterator_type, dbcsr_type_symmetric, dbcsr_get_occupation, dbcsr_get_info
59 USE erf_complex, ONLY: erfz_fast
60 USE group_dist_types, ONLY: get_group_dist, release_group_dist, group_dist_d1_type
61 USE input_constants, ONLY: &
62 casci_canonical, eri_method_full_gpw, eri_method_gpw_ht, eri_operator_coulomb, &
63 eri_operator_erf, eri_operator_erfc, eri_operator_gaussian, eri_operator_yukawa, &
64 eri_operator_trunc, eri_operator_lr_trunc, &
65 fci_solver, manual_selection, mao_projection, no_solver, qiskit_solver, wannier_projection, &
66 eri_poisson_analytic, eri_poisson_periodic, eri_poisson_mt, high_spin_roks
67 USE kpoint_types, ONLY: get_kpoint_info, kpoint_type
68 USE input_section_types, ONLY: section_vals_get, section_vals_get_subs_vals, &
69 section_vals_set_subs_vals, section_vals_type, &
70 section_vals_val_get, &
71 section_vals_val_set
72 USE iso_c_binding, ONLY: c_null_char
73 USE kinds, ONLY: default_path_length, &
74 default_string_length, &
75 dp, &
76 int_8
77 USE hfx_types, ONLY: hfx_create, hfx_release
78 USE machine, ONLY: m_walltime, m_flush
79 USE mathlib, ONLY: diamat_all
80 USE mathconstants, ONLY: fourpi, twopi, pi, rootpi
81 USE memory_utilities, ONLY: reallocate
82 USE message_passing, ONLY: mp_comm_type, &
83 mp_para_env_type, &
84 mp_para_env_release
85 USE mp2_gpw, ONLY: create_mat_munu, grep_rows_in_subgroups, build_dbcsr_from_rows
86 USE mt_util, ONLY: mt0d
87 USE parallel_gemm_api, ONLY: parallel_gemm
88 USE particle_list_types, ONLY: particle_list_type
89 USE particle_types, ONLY: particle_type
90 USE periodic_table, ONLY: ptable
91 USE physcon, ONLY: angstrom, bohr
92 USE preconditioner_types, ONLY: preconditioner_type
93 USE pw_env_methods, ONLY: pw_env_create, &
94 pw_env_rebuild
95 USE pw_env_types, ONLY: pw_env_get, &
96 pw_env_release, &
97 pw_env_type
98 USE pw_methods, ONLY: pw_integrate_function, &
99 pw_multiply, &
100 pw_multiply_with, &
101 pw_transfer, &
102 pw_zero, pw_integral_ab, pw_scale, &
103 pw_gauss_damp, pw_compl_gauss_damp
104 USE pw_poisson_methods, ONLY: pw_poisson_rebuild, &
105 pw_poisson_solve
106 USE pw_poisson_types, ONLY: analytic0d, &
107 periodic3d, &
108 greens_fn_type, &
109 pw_poisson_analytic, &
110 pw_poisson_periodic, &
111 pw_poisson_type
112 USE pw_pool_types, ONLY: &
113 pw_pool_type
114 USE pw_types, ONLY: &
115 pw_c1d_gs_type, &
116 pw_r3d_rs_type
117 USE qcschema, ONLY: qcschema_env_create, &
118 qcschema_env_release, &
119 qcschema_to_hdf5, &
120 qcschema_type
121 USE qs_active_space_fci, ONLY: solve_active_space_fci
122 USE qs_active_space_types, ONLY: active_space_type, &
123 create_active_space_type, &
124 csr_idx_from_combined, &
125 csr_idx_to_combined, &
126 eri_type, &
127 eri_type_eri_element_func
128 USE qs_active_space_mixing, ONLY: active_space_mixing_label, &
129 initialize_active_space_mixing, &
130 update_active_density
131 USE qs_active_space_utils, ONLY: eri_to_array, &
132 subspace_matrix_to_array
133 USE qs_collocate_density, ONLY: calculate_wavefunction
134 USE qs_density_matrices, ONLY: calculate_density_matrix
135 USE qs_energy_types, ONLY: qs_energy_type
136 USE qs_environment_types, ONLY: get_qs_env, &
137 qs_environment_type, &
138 set_qs_env
139 USE qs_integrate_potential, ONLY: integrate_v_rspace
140 USE qs_kind_types, ONLY: qs_kind_type
141 USE qs_ks_methods, ONLY: qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
142 evaluate_core_matrix_traces
143 USE qs_ks_types, ONLY: qs_ks_did_change, &
144 qs_ks_env_type, set_ks_env
145 USE qs_mo_io, ONLY: write_mo_set_to_output_unit
146 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
147 USE qs_mo_types, ONLY: allocate_mo_set, &
148 get_mo_set, &
149 init_mo_set, &
150 mo_set_type
151 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type, release_neighbor_list_sets
152 USE qs_ot_eigensolver, ONLY: ot_eigensolver
153 USE qs_rho_methods, ONLY: qs_rho_update_rho
154 USE qs_rho_types, ONLY: qs_rho_get, &
155 qs_rho_type
156 USE qs_subsys_types, ONLY: qs_subsys_get, &
157 qs_subsys_type
158 USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
159 USE scf_control_types, ONLY: scf_control_type
160#ifndef __NO_SOCKETS
161 USE sockets_interface, ONLY: accept_socket, &
162 close_socket, &
163 listen_socket, &
164 open_bind_socket, &
165 readbuffer, &
166 remove_socket_file, &
167 writebuffer
168#endif
169 USE task_list_methods, ONLY: generate_qs_task_list
170 USE task_list_types, ONLY: allocate_task_list, &
171 deallocate_task_list, &
172 task_list_type
173 USE util, ONLY: get_limit
174#include "./base/base_uses.f90"
175
176 IMPLICIT NONE
177 PRIVATE
178
179 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_methods'
180
181 PUBLIC :: active_space_main
182
183 TYPE, EXTENDS(eri_type_eri_element_func) :: eri_fcidump_print
184 INTEGER :: unit_nr = -1, bra_start = -1, ket_start = -1
185 CONTAINS
186 PROCEDURE :: func => eri_fcidump_print_func
187 END TYPE eri_fcidump_print
188
189 TYPE, EXTENDS(eri_type_eri_element_func) :: eri_fcidump_checksum
190 INTEGER :: bra_start = 0, ket_start = 0
191 REAL(KIND=dp) :: checksum = 0.0_dp
192 CONTAINS
193 PROCEDURE, PASS :: set => eri_fcidump_set
194 PROCEDURE :: func => eri_fcidump_checksum_func
195 END TYPE eri_fcidump_checksum
196
197CONTAINS
198
199! **************************************************************************************************
200!> \brief Sets the starting indices of the bra and ket.
201!> \param this object reference
202!> \param bra_start starting index of the bra
203!> \param ket_start starting index of the ket
204! **************************************************************************************************
205 SUBROUTINE eri_fcidump_set(this, bra_start, ket_start)
206 CLASS(eri_fcidump_checksum) :: this
207 INTEGER, INTENT(IN) :: bra_start, ket_start
208 this%bra_start = bra_start
209 this%ket_start = ket_start
210 END SUBROUTINE eri_fcidump_set
211
212! **************************************************************************************************
213!> \brief Main method for determining the active space Hamiltonian
214!> \param qs_env ...
215! **************************************************************************************************
216 SUBROUTINE active_space_main(qs_env)
217 TYPE(qs_environment_type), POINTER :: qs_env
218
219 CHARACTER(len=*), PARAMETER :: routinen = 'active_space_main'
220
221 CHARACTER(len=10) :: cshell, lnam(5)
222 CHARACTER(len=default_path_length) :: qcschema_filename
223 CHARACTER(LEN=default_string_length) :: basis_type, kp_scheme
224 INTEGER :: as_solver, eri_method, eri_operator, eri_print, group_size, handle, i, iatom, &
225 ishell, isp, ispin, iw, j, jm, m, max_orb_ind, mselect, n1, n2, nao, natom, nel, &
226 nelec_active, nelec_inactive, nelec_total, nkp, nmo, nmo_active, nmo_available, &
227 nmo_inactive, nmo_inactive_remaining, nmo_occ, nmo_virtual, nn1, nn2, nrow_global, nspins
228 INTEGER, DIMENSION(5) :: nshell
229 INTEGER, DIMENSION(:), POINTER :: invals
230 LOGICAL :: do_ddapc, do_kpoints, ex_omega, &
231 ex_operator, ex_perd, ex_rcut, &
232 explicit, stop_after_print, store_wfn, &
233 use_real_wfn
234 REAL(kind=dp) :: eri_eps_filter, eri_eps_grid, eri_eps_int, eri_gpw_cutoff, eri_op_omega, &
235 eri_rcut, eri_rel_cutoff, fel, focc, maxocc, nze_percentage
236 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
237 REAL(kind=dp), DIMENSION(:), POINTER :: evals_virtual
238 TYPE(active_space_type), POINTER :: active_space_env
239 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
240 TYPE(cell_type), POINTER :: cell
241 TYPE(cp_blacs_env_type), POINTER :: context
242 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
243 TYPE(cp_fm_type) :: fm_dummy, mo_virtual
244 TYPE(cp_fm_type), POINTER :: fm_target_active, fm_target_inactive, &
245 fmat, mo_coeff, mo_ref, mo_target
246 TYPE(cp_logger_type), POINTER :: logger
247 TYPE(dbcsr_csr_type), POINTER :: eri_mat
248 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, rho_ao, s_matrix
249 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix_kp, rho_ao_kp, s_matrix_kp
250 TYPE(dbcsr_type), POINTER :: denmat
251 TYPE(dft_control_type), POINTER :: dft_control
252 TYPE(kpoint_type), POINTER :: kpoints
253 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
254 TYPE(mo_set_type), POINTER :: mo_set, mo_set_active, mo_set_inactive
255 TYPE(mp_para_env_type), POINTER :: para_env
256 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
257 TYPE(preconditioner_type), POINTER :: local_preconditioner
258 TYPE(qcschema_type) :: qcschema_env
259 TYPE(qs_energy_type), POINTER :: energy
260 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
261 TYPE(qs_ks_env_type), POINTER :: ks_env
262 TYPE(qs_rho_type), POINTER :: rho
263 TYPE(scf_control_type), POINTER :: scf_control
264 TYPE(section_vals_type), POINTER :: adiabatic_rescaling, as_input, &
265 hfx_section, input, loc_print, &
266 loc_section, print_orb, xc_section
267
268 !--------------------------------------------------------------------------------------------!
269
270 CALL get_qs_env(qs_env, input=input)
271 as_input => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE")
272 CALL section_vals_get(as_input, explicit=explicit)
273 IF (.NOT. explicit) RETURN
274 CALL timeset(routinen, handle)
275
276 logger => cp_get_default_logger()
277 iw = cp_logger_get_default_io_unit(logger)
278
279 IF (iw > 0) THEN
280 WRITE (iw, '(/,T2,A)') &
281 '!-----------------------------------------------------------------------------!'
282 WRITE (iw, '(T26,A)') "Active Space Embedding Module"
283 WRITE (iw, '(T2,A)') &
284 '!-----------------------------------------------------------------------------!'
285 END IF
286
287 ! k-points?
288 NULLIFY (kpoints)
289 CALL get_qs_env(qs_env, do_kpoints=do_kpoints, dft_control=dft_control, kpoints=kpoints)
290 IF (do_kpoints) THEN
291 IF (.NOT. ASSOCIATED(kpoints)) THEN
292 CALL cp_abort(__location__, "Missing Gamma-point environment for active space module")
293 END IF
294 CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme, nkp=nkp, use_real_wfn=use_real_wfn)
295 IF (trim(kp_scheme) /= "GAMMA" .OR. nkp /= 1 .OR. .NOT. use_real_wfn) THEN
296 CALL cp_abort(__location__, &
297 "Only Gamma-point DFT%KPOINTS are supported in the active space module")
298 END IF
299 IF (.NOT. ASSOCIATED(kpoints%kp_env)) THEN
300 CALL cp_abort(__location__, "Missing Gamma-point environment for active space module")
301 END IF
302 IF (.NOT. ASSOCIATED(kpoints%kp_env(1)%kpoint_env)) THEN
303 CALL cp_abort(__location__, "Missing Gamma-point environment for active space module")
304 END IF
305 IF (.NOT. ASSOCIATED(kpoints%kp_env(1)%kpoint_env%mos)) THEN
306 CALL cp_abort(__location__, "Missing Gamma-point MOs for active space module")
307 END IF
308 END IF
309
310 ! adiabatic rescaling?
311 adiabatic_rescaling => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
312 CALL section_vals_get(adiabatic_rescaling, explicit=explicit)
313 IF (explicit) THEN
314 CALL cp_abort(__location__, "Adiabatic rescaling not supported in active space module")
315 END IF
316
317 ! Setup the possible usage of DDAPC charges
318 do_ddapc = dft_control%qs_control%ddapc_restraint .OR. &
319 qs_env%cp_ddapc_ewald%do_decoupling .OR. &
320 qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
321 qs_env%cp_ddapc_ewald%do_solvation
322 IF (do_ddapc) THEN
323 CALL cp_abort(__location__, "DDAPC charges are not supported in the active space module")
324 END IF
325 IF (dft_control%do_sccs) THEN
326 CALL cp_abort(__location__, "SCCS is not supported in the active space module")
327 END IF
328 IF (dft_control%correct_surf_dip) THEN
329 IF (dft_control%surf_dip_correct_switch) THEN
330 CALL cp_abort(__location__, "Surface dipole correction not supported in the AS module")
331 END IF
332 END IF
333 IF (dft_control%smeagol_control%smeagol_enabled) THEN
334 CALL cp_abort(__location__, "SMEAGOL is not supported in the active space module")
335 END IF
336 IF (dft_control%qs_control%do_kg) THEN
337 CALL cp_abort(__location__, "KG correction not supported in the active space module")
338 END IF
339
340 NULLIFY (active_space_env)
341 CALL create_active_space_type(active_space_env)
342 active_space_env%energy_total = 0.0_dp
343 active_space_env%energy_ref = 0.0_dp
344 active_space_env%energy_inactive = 0.0_dp
345 active_space_env%energy_active = 0.0_dp
346
347 ! input options
348
349 ! figure out what needs to be printed/stored
350 IF (btest(cp_print_key_should_output(logger%iter_info, as_input, "FCIDUMP"), cp_p_file)) THEN
351 active_space_env%fcidump = .true.
352 END IF
353
354 CALL section_vals_val_get(as_input, "QCSCHEMA", c_val=qcschema_filename, explicit=explicit)
355 IF (explicit) THEN
356 active_space_env%qcschema = .true.
357 active_space_env%qcschema_filename = qcschema_filename
358 END IF
359
360 CALL section_vals_val_get(as_input, "ACTIVE_ELECTRONS", i_val=nelec_active)
361 CALL get_qs_env(qs_env, nelectron_total=nelec_total)
362
363 IF (nelec_active <= 0) cpabort("Specify a positive number of active electrons.")
364 IF (nelec_active > nelec_total) cpabort("More active electrons than total electrons.")
365
366 nelec_inactive = nelec_total - nelec_active
367 IF (mod(nelec_inactive, 2) /= 0) THEN
368 cpabort("The remaining number of inactive electrons has to be even.")
369 END IF
370
371 IF (iw > 0) THEN
372 WRITE (iw, '(T3,A,T70,I10)') "Total number of electrons", nelec_total
373 WRITE (iw, '(T3,A,T70,I10)') "Number of inactive electrons", nelec_inactive
374 WRITE (iw, '(T3,A,T70,I10)') "Number of active electrons", nelec_active
375 END IF
376
377 CALL get_qs_env(qs_env, dft_control=dft_control)
378 nspins = dft_control%nspins
379
380 active_space_env%nelec_active = nelec_active
381 active_space_env%nelec_inactive = nelec_inactive
382 active_space_env%nelec_total = nelec_total
383 active_space_env%nspins = nspins
384 active_space_env%multiplicity = dft_control%multiplicity
385 active_space_env%restricted_orbitals = dft_control%roks
386
387 ! define the active/inactive space orbitals
388 CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit, i_val=nmo_active)
389 IF (.NOT. explicit) THEN
390 CALL cp_abort(__location__, "Number of Active Orbitals has to be specified.")
391 END IF
392 active_space_env%nmo_active = nmo_active
393 ! this is safe because nelec_inactive is always even
394 nmo_inactive = nelec_inactive/2
395 active_space_env%nmo_inactive = nmo_inactive
396
397 CALL initialize_active_space_mixing(active_space_env, as_input)
398
399 CALL section_vals_val_get(as_input, "ORBITAL_SELECTION", i_val=mselect)
400 IF (iw > 0) THEN
401 SELECT CASE (mselect)
402 CASE DEFAULT
403 cpabort("Unknown orbital selection method")
404 CASE (casci_canonical)
405 WRITE (iw, '(/,T3,A)') &
406 "Active space orbitals selected using energy ordered canonical orbitals"
407 CASE (wannier_projection)
408 WRITE (iw, '(/,T3,A)') &
409 "Active space orbitals selected using projected Wannier orbitals"
410 CASE (mao_projection)
411 WRITE (iw, '(/,T3,A)') &
412 "Active space orbitals selected using modified atomic orbitals (MAO)"
413 CASE (manual_selection)
414 WRITE (iw, '(/,T3,A)') &
415 "Active space orbitals selected manually"
416 END SELECT
417
418 WRITE (iw, '(T3,A,T70,I10)') "Number of inactive orbitals", nmo_inactive
419 WRITE (iw, '(T3,A,T70,I10)') "Number of active orbitals", nmo_active
420 END IF
421
422 ! get projection spaces
423 CALL section_vals_val_get(as_input, "SUBSPACE_ATOM", i_val=iatom, explicit=explicit)
424 IF (explicit) THEN
425 CALL get_qs_env(qs_env, natom=natom)
426 IF (iatom <= 0 .OR. iatom > natom) THEN
427 IF (iw > 0) THEN
428 WRITE (iw, '(/,T3,A,I3)') "ERROR: SUBSPACE_ATOM number is not valid", iatom
429 END IF
430 cpabort("Select a valid SUBSPACE_ATOM")
431 END IF
432 END IF
433 CALL section_vals_val_get(as_input, "SUBSPACE_SHELL", c_val=cshell, explicit=explicit)
434 nshell = 0
435 lnam = ""
436 IF (explicit) THEN
437 cshell = adjustl(cshell)
438 n1 = 1
439 DO i = 1, 5
440 ishell = i
441 IF (cshell(n1:n1) == " ") THEN
442 ishell = ishell - 1
443 EXIT
444 END IF
445 READ (cshell(n1:), "(I1,A1)") nshell(i), lnam(i)
446 n1 = n1 + 2
447 END DO
448 END IF
449
450 ! generate orbitals
451 SELECT CASE (mselect)
452 CASE DEFAULT
453 cpabort("Unknown orbital selection method")
454 CASE (casci_canonical)
455 IF (do_kpoints) THEN
456 mos => kpoints%kp_env(1)%kpoint_env%mos(1, :)
457 ELSE
458 CALL get_qs_env(qs_env, mos=mos)
459 END IF
460
461 ! total number of occupied orbitals, i.e. inactive plus active MOs
462 nmo_occ = nmo_inactive + nmo_active
463
464 ! set inactive orbital indices, these are trivially 1...nmo_inactive
465 ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
466 DO ispin = 1, nspins
467 DO i = 1, nmo_inactive
468 active_space_env%inactive_orbitals(i, ispin) = i
469 END DO
470 END DO
471
472 ! set active orbital indices, these are shifted by nmo_inactive
473 ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
474 DO ispin = 1, nspins
475 DO i = 1, nmo_active
476 active_space_env%active_orbitals(i, ispin) = nmo_inactive + i
477 END DO
478 END DO
479
480 ! allocate and initialize inactive and active mo coefficients.
481 ! These are stored in a data structure for the full occupied space:
482 ! for inactive mos, the active subset is set to zero, vice versa for the active mos
483 ! TODO: allocate data structures only for the eaxct number MOs
484 maxocc = 2.0_dp
485 IF (nspins > 1) maxocc = 1.0_dp
486 ALLOCATE (active_space_env%mos_active(nspins))
487 ALLOCATE (active_space_env%mos_inactive(nspins))
488 DO ispin = 1, nspins
489 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
490 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
491 ! the right number of active electrons per spin channel is initialized further down
492 CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo_occ, 0, 0.0_dp, maxocc, 0.0_dp)
493 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
494 nrow_global=nrow_global, ncol_global=nmo_occ)
495 CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
496 CALL cp_fm_struct_release(fm_struct_tmp)
497 IF (nspins == 2) THEN
498 nel = nelec_inactive/2
499 ELSE
500 nel = nelec_inactive
501 END IF
502 CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo_occ, nel, &
503 REAL(nel, kind=dp), maxocc, 0.0_dp)
504 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
505 nrow_global=nrow_global, ncol_global=nmo_occ)
506 CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO")
507 CALL cp_fm_struct_release(fm_struct_tmp)
508 END DO
509
510 ! create canonical orbitals
511 CALL get_qs_env(qs_env, scf_control=scf_control)
512 IF (dft_control%roks .AND. scf_control%roks_scheme /= high_spin_roks) THEN
513 CALL cp_abort(__location__, &
514 "Only high-spin ROKS is supported for ACTIVE_SPACE FCI; "// &
515 "general ROKS MO definitions are not implemented.")
516 ELSE
517 IF (dft_control%do_admm) THEN
518 IF (dft_control%do_admm_mo) THEN
519 cpabort("ADMM currently possible only with purification none_dm")
520 END IF
521 END IF
522
523 ALLOCATE (eigenvalues(nmo_occ, nspins))
524 eigenvalues = 0.0_dp
525 IF (do_kpoints) THEN
526 CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix_kp, matrix_s_kp=s_matrix_kp, &
527 scf_control=scf_control)
528 ks_matrix => ks_matrix_kp(:, 1)
529 s_matrix => s_matrix_kp(:, 1)
530 ELSE
531 CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
532 END IF
533
534 ! calculate virtual MOs and copy inactive and active orbitals
535 IF (iw > 0) THEN
536 WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
537 END IF
538 DO ispin = 1, nspins
539 ! nmo_available is the number of MOs available from the SCF calculation:
540 ! this is at least the number of occupied orbitals in the SCF, plus
541 ! any number of added MOs (virtuals) requested in the SCF section
542 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
543
544 ! calculate how many extra MOs we still have to compute
545 nmo_virtual = nmo_occ - nmo_available
546 nmo_virtual = max(nmo_virtual, 0)
547
548 NULLIFY (evals_virtual)
549 ALLOCATE (evals_virtual(nmo_virtual))
550
551 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
552 nrow_global=nrow_global)
553
554 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
555 nrow_global=nrow_global, ncol_global=nmo_virtual)
556 CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
557 CALL cp_fm_struct_release(fm_struct_tmp)
558 CALL cp_fm_init_random(mo_virtual, nmo_virtual)
559
560 NULLIFY (local_preconditioner)
561
562 ! compute missing virtual MOs
563 CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
564 matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
565 eps_gradient=scf_control%eps_lumos, &
566 preconditioner=local_preconditioner, &
567 iter_max=scf_control%max_iter_lumos, &
568 size_ortho_space=nmo_available)
569
570 ! get the eigenvalues
571 CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, evals_virtual)
572
573 ! we need to send the copy of MOs to preserve the sign
574 CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct)
575 CALL cp_fm_to_fm(mo_ref, fm_dummy)
576 CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, &
577 evals_arg=eigenvalues(:, ispin), do_rotation=.true.)
578
579 ! copy inactive orbitals
580 mo_set => active_space_env%mos_inactive(ispin)
581 CALL get_mo_set(mo_set, mo_coeff=mo_target)
582 DO i = 1, SIZE(active_space_env%inactive_orbitals, 1)
583 m = active_space_env%inactive_orbitals(i, ispin)
584 CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
585 mo_set%eigenvalues(m) = eigenvalues(m, ispin)
586 IF (nspins > 1) THEN
587 mo_set%occupation_numbers(m) = 1.0
588 ELSE
589 mo_set%occupation_numbers(m) = 2.0
590 END IF
591 END DO
592
593 ! copy active orbitals
594 mo_set => active_space_env%mos_active(ispin)
595 CALL get_mo_set(mo_set, mo_coeff=mo_target)
596 ! for mult > 1, put the polarized electrons in the alpha channel
597 IF (nspins == 2) THEN
598 IF (ispin == 1) THEN
599 nel = (nelec_active + active_space_env%multiplicity - 1)/2
600 ELSE
601 nel = (nelec_active - active_space_env%multiplicity + 1)/2
602 END IF
603 ELSE
604 nel = nelec_active
605 END IF
606 mo_set%nelectron = nel
607 mo_set%n_el_f = real(nel, kind=dp)
608 DO i = 1, nmo_active
609 m = active_space_env%active_orbitals(i, ispin)
610 IF (m > nmo_available) THEN
611 CALL cp_fm_to_fm(mo_virtual, mo_target, 1, m - nmo_available, m)
612 eigenvalues(m, ispin) = evals_virtual(m - nmo_available)
613 mo_set%occupation_numbers(m) = 0.0
614 ELSE
615 CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
616 mo_set%occupation_numbers(m) = mos(ispin)%occupation_numbers(m)
617 END IF
618 mo_set%eigenvalues(m) = eigenvalues(m, ispin)
619 END DO
620 ! Release
621 DEALLOCATE (evals_virtual)
622 CALL cp_fm_release(fm_dummy)
623 CALL cp_fm_release(mo_virtual)
624 END DO
625
626 IF (iw > 0) THEN
627 DO ispin = 1, nspins
628 WRITE (iw, '(/,T3,A,I3,T66,A)') "Canonical Orbital Selection for spin", ispin, &
629 "[atomic units]"
630 DO i = 1, nmo_inactive, 4
631 jm = min(3, nmo_inactive - i)
632 WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [I]", j=0, jm)
633 END DO
634 DO i = nmo_inactive + 1, nmo_inactive + nmo_active, 4
635 jm = min(3, nmo_inactive + nmo_active - i)
636 WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [A]", j=0, jm)
637 END DO
638 WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
639 DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
640 jm = min(3, SIZE(active_space_env%active_orbitals, 1) - i)
641 WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
642 END DO
643 END DO
644 END IF
645 DEALLOCATE (eigenvalues)
646 END IF
647
648 CASE (manual_selection)
649 ! create canonical orbitals
650 IF (dft_control%roks) THEN
651 CALL cp_abort(__location__, &
652 "Manual ACTIVE_SPACE orbital selection is not supported for ROKS; "// &
653 "use canonical high-spin ROKS.")
654 ELSE
655 IF (dft_control%do_admm) THEN
656 ! For admm_mo, the auxiliary density is computed from the MOs, which never change
657 ! in the rs-dft embedding, therefore the energy is wrong as the LR HFX never changes.
658 ! For admm_dm, the auxiliary density is computed from the density matrix, which is
659 ! updated at each iteration and therefore works.
660 IF (dft_control%do_admm_mo) THEN
661 cpabort("ADMM currently possible only with purification none_dm")
662 END IF
663 END IF
664
665 CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", explicit=explicit, i_vals=invals)
666 IF (.NOT. explicit) THEN
667 CALL cp_abort(__location__, "Manual orbital selection requires to explicitly "// &
668 "set the active orbital indices via ACTIVE_ORBITAL_INDICES")
669 END IF
670
671 IF (nspins == 1) THEN
672 cpassert(SIZE(invals) == nmo_active)
673 ELSE
674 cpassert(SIZE(invals) == 2*nmo_active)
675 END IF
676 ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
677 ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
678
679 DO ispin = 1, nspins
680 DO i = 1, nmo_active
681 active_space_env%active_orbitals(i, ispin) = invals(i + (ispin - 1)*nmo_active)
682 END DO
683 END DO
684
685 IF (do_kpoints) THEN
686 mos => kpoints%kp_env(1)%kpoint_env%mos(1, :)
687 ELSE
688 CALL get_qs_env(qs_env, mos=mos)
689 END IF
690
691 ! include MOs up to the largest index in the list
692 max_orb_ind = maxval(invals)
693 maxocc = 2.0_dp
694 IF (nspins > 1) maxocc = 1.0_dp
695 ALLOCATE (active_space_env%mos_active(nspins))
696 ALLOCATE (active_space_env%mos_inactive(nspins))
697 DO ispin = 1, nspins
698 ! init active orbitals
699 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
700 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
701 CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, max_orb_ind, 0, 0.0_dp, maxocc, 0.0_dp)
702 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
703 nrow_global=nrow_global, ncol_global=max_orb_ind)
704 CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
705 CALL cp_fm_struct_release(fm_struct_tmp)
706
707 ! init inactive orbitals
708 IF (nspins == 2) THEN
709 nel = nelec_inactive/2
710 ELSE
711 nel = nelec_inactive
712 END IF
713 CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, max_orb_ind, nel, real(nel, kind=dp), maxocc, 0.0_dp)
714 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
715 nrow_global=nrow_global, ncol_global=max_orb_ind)
716 CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO")
717 ! small hack: set the correct inactive occupations down below
718 active_space_env%mos_inactive(ispin)%occupation_numbers = 0.0_dp
719 CALL cp_fm_struct_release(fm_struct_tmp)
720 END DO
721
722 ALLOCATE (eigenvalues(max_orb_ind, nspins))
723 eigenvalues = 0.0_dp
724 IF (do_kpoints) THEN
725 CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix_kp, matrix_s_kp=s_matrix_kp, &
726 scf_control=scf_control)
727 ks_matrix => ks_matrix_kp(:, 1)
728 s_matrix => s_matrix_kp(:, 1)
729 ELSE
730 CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
731 END IF
732
733 ! calculate virtual MOs and copy inactive and active orbitals
734 IF (iw > 0) THEN
735 WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
736 END IF
737 DO ispin = 1, nspins
738 CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
739 nmo_virtual = max_orb_ind - nmo_available
740 nmo_virtual = max(nmo_virtual, 0)
741
742 NULLIFY (evals_virtual)
743 ALLOCATE (evals_virtual(nmo_virtual))
744
745 CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
746 nrow_global=nrow_global)
747
748 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
749 nrow_global=nrow_global, ncol_global=nmo_virtual)
750 CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
751 CALL cp_fm_struct_release(fm_struct_tmp)
752 CALL cp_fm_init_random(mo_virtual, nmo_virtual)
753
754 NULLIFY (local_preconditioner)
755
756 CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
757 matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
758 eps_gradient=scf_control%eps_lumos, &
759 preconditioner=local_preconditioner, &
760 iter_max=scf_control%max_iter_lumos, &
761 size_ortho_space=nmo_available)
762
763 CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, &
764 evals_virtual)
765
766 ! We need to send the copy of MOs to preserve the sign
767 CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct)
768 CALL cp_fm_to_fm(mo_ref, fm_dummy)
769
770 CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, &
771 evals_arg=eigenvalues(:, ispin), do_rotation=.true.)
772
773 mo_set_active => active_space_env%mos_active(ispin)
774 CALL get_mo_set(mo_set_active, mo_coeff=fm_target_active)
775 mo_set_inactive => active_space_env%mos_inactive(ispin)
776 CALL get_mo_set(mo_set_inactive, mo_coeff=fm_target_inactive)
777
778 ! copy orbitals
779 nmo_inactive_remaining = nmo_inactive
780 DO i = 1, max_orb_ind
781 ! case for i being an active orbital
782 IF (any(active_space_env%active_orbitals(:, ispin) == i)) THEN
783 IF (i > nmo_available) THEN
784 CALL cp_fm_to_fm(mo_virtual, fm_target_active, 1, i - nmo_available, i)
785 eigenvalues(i, ispin) = evals_virtual(i - nmo_available)
786 mo_set_active%occupation_numbers(i) = 0.0
787 ELSE
788 CALL cp_fm_to_fm(fm_dummy, fm_target_active, 1, i, i)
789 mo_set_active%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
790 END IF
791 mo_set_active%eigenvalues(i) = eigenvalues(i, ispin)
792 ! if it was not an active orbital, check whether it is an inactive orbital
793 ELSE IF (nmo_inactive_remaining > 0) THEN
794 CALL cp_fm_to_fm(fm_dummy, fm_target_inactive, 1, i, i)
795 ! store on the fly the mapping of inactive orbitals
796 active_space_env%inactive_orbitals(nmo_inactive - nmo_inactive_remaining + 1, ispin) = i
797 mo_set_inactive%eigenvalues(i) = eigenvalues(i, ispin)
798 mo_set_inactive%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
799 ! hack: set homo and lumo manually
800 IF (nmo_inactive_remaining == 1) THEN
801 mo_set_inactive%homo = i
802 mo_set_inactive%lfomo = i + 1
803 END IF
804 nmo_inactive_remaining = nmo_inactive_remaining - 1
805 ELSE
806 cycle
807 END IF
808 END DO
809
810 ! Release
811 DEALLOCATE (evals_virtual)
812 CALL cp_fm_release(fm_dummy)
813 CALL cp_fm_release(mo_virtual)
814 END DO
815
816 IF (iw > 0) THEN
817 DO ispin = 1, nspins
818 WRITE (iw, '(/,T3,A,I3,T66,A)') "Orbital Energies and Selection for spin", ispin, "[atomic units]"
819
820 DO i = 1, max_orb_ind, 4
821 jm = min(3, max_orb_ind - i)
822 WRITE (iw, '(T4)', advance="no")
823 DO j = 0, jm
824 IF (any(active_space_env%active_orbitals(:, ispin) == i + j)) THEN
825 WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [A]"
826 ELSE IF (any(active_space_env%inactive_orbitals(:, ispin) == i + j)) THEN
827 WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [I]"
828 ELSE
829 WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [V]"
830 END IF
831 END DO
832 WRITE (iw, *)
833 END DO
834 WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
835 DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
836 jm = min(3, SIZE(active_space_env%active_orbitals, 1) - i)
837 WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
838 END DO
839 END DO
840 END IF
841 DEALLOCATE (eigenvalues)
842 END IF
843
844 CASE (wannier_projection)
845 NULLIFY (loc_section, loc_print)
846 loc_section => section_vals_get_subs_vals(as_input, "LOCALIZE")
847 cpassert(ASSOCIATED(loc_section))
848 loc_print => section_vals_get_subs_vals(as_input, "LOCALIZE%PRINT")
849 !
850 cpabort("not yet available")
851 !
852 CASE (mao_projection)
853 !
854 cpabort("not yet available")
855 !
856 END SELECT
857
858 ! Print orbitals on Cube files
859 print_orb => section_vals_get_subs_vals(as_input, "PRINT_ORBITAL_CUBES")
860 CALL section_vals_get(print_orb, explicit=explicit)
861 CALL section_vals_val_get(print_orb, "STOP_AFTER_CUBES", l_val=stop_after_print)
862 IF (explicit) THEN
863 !
864 CALL print_orbital_cubes(print_orb, qs_env, active_space_env%mos_active)
865 !
866 IF (stop_after_print) THEN
867
868 IF (iw > 0) THEN
869 WRITE (iw, '(/,T2,A)') &
870 '!----------------- Early End of Active Space Interface -----------------------!'
871 END IF
872
873 CALL timestop(handle)
874
875 RETURN
876 END IF
877 END IF
878
879 ! calculate inactive density matrix
880 CALL get_qs_env(qs_env, rho=rho)
881 IF (do_kpoints) THEN
882 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
883 rho_ao => rho_ao_kp(:, 1)
884 ELSE
885 CALL qs_rho_get(rho, rho_ao=rho_ao)
886 END IF
887 cpassert(ASSOCIATED(rho_ao))
888 CALL dbcsr_allocate_matrix_set(active_space_env%pmat_inactive, nspins)
889 DO ispin = 1, nspins
890 ALLOCATE (denmat)
891 CALL dbcsr_copy(denmat, rho_ao(ispin)%matrix)
892 mo_set => active_space_env%mos_inactive(ispin)
893 CALL calculate_density_matrix(mo_set, denmat)
894 active_space_env%pmat_inactive(ispin)%matrix => denmat
895 END DO
896
897 ! read in ERI parameters
898 CALL section_vals_val_get(as_input, "ERI%METHOD", i_val=eri_method)
899 active_space_env%eri%method = eri_method
900 CALL section_vals_val_get(as_input, "ERI%OPERATOR", i_val=eri_operator, explicit=ex_operator)
901 active_space_env%eri%operator = eri_operator
902 CALL section_vals_val_get(as_input, "ERI%OMEGA", r_val=eri_op_omega, explicit=ex_omega)
903 active_space_env%eri%omega = eri_op_omega
904 CALL section_vals_val_get(as_input, "ERI%CUTOFF_RADIUS", r_val=eri_rcut, explicit=ex_rcut)
905 active_space_env%eri%cutoff_radius = eri_rcut ! this is already converted to bohr!
906 CALL section_vals_val_get(as_input, "ERI%PERIODICITY", i_vals=invals, explicit=ex_perd)
907 CALL section_vals_val_get(as_input, "ERI%EPS_INTEGRAL", r_val=eri_eps_int)
908 active_space_env%eri%eps_integral = eri_eps_int
909 ! if eri periodicity is explicitly set, we use it, otherwise we use the cell periodicity
910 IF (ex_perd) THEN
911 IF (SIZE(invals) == 1) THEN
912 active_space_env%eri%periodicity(1:3) = invals(1)
913 ELSE
914 active_space_env%eri%periodicity(1:3) = invals(1:3)
915 END IF
916 ELSE
917 CALL get_qs_env(qs_env, cell=cell)
918 active_space_env%eri%periodicity(1:3) = cell%perd(1:3)
919 END IF
920 IF (iw > 0) THEN
921 WRITE (iw, '(/,T3,A)') "Calculation of Electron Repulsion Integrals"
922
923 SELECT CASE (eri_method)
924 CASE (eri_method_full_gpw)
925 WRITE (iw, '(T3,A,T50,A)') "Integration method", "GPW Fourier transform over MOs"
926 CASE (eri_method_gpw_ht)
927 WRITE (iw, '(T3,A,T44,A)') "Integration method", "Half transformed integrals from GPW"
928 CASE DEFAULT
929 cpabort("Unknown ERI method")
930 END SELECT
931
932 SELECT CASE (eri_operator)
933 CASE (eri_operator_coulomb)
934 WRITE (iw, '(T3,A,T73,A)') "ERI operator", "Coulomb"
935
936 CASE (eri_operator_yukawa)
937 WRITE (iw, '(T3,A,T74,A)') "ERI operator", "Yukawa"
938 IF (.NOT. ex_omega) CALL cp_abort(__location__, &
939 "Yukawa operator requires OMEGA to be explicitly set")
940 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
941
942 CASE (eri_operator_erf)
943 WRITE (iw, '(T3,A,T63,A)') "ERI operator", "Longrange Coulomb"
944 IF (.NOT. ex_omega) CALL cp_abort(__location__, &
945 "Longrange operator requires OMEGA to be explicitly set")
946 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
947
948 CASE (eri_operator_erfc)
949 WRITE (iw, '(T3,A,T62,A)') "ERI operator", "Shortrange Coulomb"
950 IF (.NOT. ex_omega) CALL cp_abort(__location__, &
951 "Shortrange operator requires OMEGA to be explicitly set")
952 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
953
954 CASE (eri_operator_trunc)
955 WRITE (iw, '(T3,A,T63,A)') "ERI operator", "Truncated Coulomb"
956 IF (.NOT. ex_rcut) CALL cp_abort(__location__, &
957 "Cutoff radius not specified for trunc. Coulomb operator")
958 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator cutoff radius (au)", eri_rcut
959
960 CASE (eri_operator_lr_trunc)
961 WRITE (iw, '(T3,A,T53,A)') "ERI operator", "Longrange truncated Coulomb"
962 IF (.NOT. ex_rcut) CALL cp_abort(__location__, &
963 "Cutoff radius not specified for trunc. longrange operator")
964 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator cutoff radius (au)", eri_rcut
965 IF (.NOT. ex_omega) CALL cp_abort(__location__, &
966 "LR truncated operator requires OMEGA to be explicitly set")
967 WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
968 IF (eri_op_omega < 0.01_dp) THEN
969 cpabort("LR truncated operator requires OMEGA >= 0.01 to be stable")
970 END IF
971
972 CASE DEFAULT
973 cpabort("Unknown ERI operator")
974
975 END SELECT
976
977 WRITE (iw, '(T3,A,T68,E12.4)') "Accuracy of ERIs", eri_eps_int
978 WRITE (iw, '(T3,A,T71,3I3)') "Periodicity", active_space_env%eri%periodicity(1:3)
979
980 ! TODO: should be moved after ERI calculation, as it depends on screening
981 IF (nspins < 2) THEN
982 WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI", (nmo_active**4)/8
983 ELSE
984 WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|aa)", (nmo_active**4)/8
985 WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (bb|bb)", (nmo_active**4)/8
986 WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|bb)", (nmo_active**4)/4
987 END IF
988 END IF
989
990 ! allocate container for integrals (CSR matrix)
991 CALL get_qs_env(qs_env, para_env=para_env)
992 m = (nspins*(nspins + 1))/2
993 ! With ROHF/ROKS, we need ERIs from only a single set of orbitals
994 IF (dft_control%roks) m = 1
995 ALLOCATE (active_space_env%eri%eri(m))
996 DO i = 1, m
997 CALL get_mo_set(active_space_env%mos_active(1), nmo=nmo)
998 ALLOCATE (active_space_env%eri%eri(i)%csr_mat)
999 eri_mat => active_space_env%eri%eri(i)%csr_mat
1000 IF (i == 1) THEN
1001 n1 = nmo
1002 n2 = nmo
1003 ELSE IF (i == 2) THEN
1004 n1 = nmo
1005 n2 = nmo
1006 ELSE
1007 n1 = nmo
1008 n2 = nmo
1009 END IF
1010 nn1 = (n1*(n1 + 1))/2
1011 nn2 = (n2*(n2 + 1))/2
1012 CALL dbcsr_csr_create(eri_mat, nn1, nn2, 0_int_8, 0, 0, para_env%get_handle())
1013 active_space_env%eri%norb = nmo
1014 END DO
1015
1016 SELECT CASE (eri_method)
1017 CASE (eri_method_full_gpw, eri_method_gpw_ht)
1018 CALL section_vals_val_get(as_input, "ERI_GPW%EPS_GRID", r_val=eri_eps_grid)
1019 active_space_env%eri%eri_gpw%eps_grid = eri_eps_grid
1020 CALL section_vals_val_get(as_input, "ERI_GPW%EPS_FILTER", r_val=eri_eps_filter)
1021 active_space_env%eri%eri_gpw%eps_filter = eri_eps_filter
1022 CALL section_vals_val_get(as_input, "ERI_GPW%CUTOFF", r_val=eri_gpw_cutoff)
1023 active_space_env%eri%eri_gpw%cutoff = eri_gpw_cutoff
1024 CALL section_vals_val_get(as_input, "ERI_GPW%REL_CUTOFF", r_val=eri_rel_cutoff)
1025 active_space_env%eri%eri_gpw%rel_cutoff = eri_rel_cutoff
1026 CALL section_vals_val_get(as_input, "ERI_GPW%PRINT_LEVEL", i_val=eri_print)
1027 active_space_env%eri%eri_gpw%print_level = eri_print
1028 CALL section_vals_val_get(as_input, "ERI_GPW%STORE_WFN", l_val=store_wfn)
1029 active_space_env%eri%eri_gpw%store_wfn = store_wfn
1030 CALL section_vals_val_get(as_input, "ERI_GPW%GROUP_SIZE", i_val=group_size)
1031 active_space_env%eri%eri_gpw%group_size = group_size
1032 ! Always redo Poisson solver for now
1033 active_space_env%eri%eri_gpw%redo_poisson = .true.
1034 ! active_space_env%eri%eri_gpw%redo_poisson = (ex_operator .OR. ex_perd)
1035 IF (iw > 0) THEN
1036 WRITE (iw, '(/,T2,A,T71,F10.1)') "ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff
1037 WRITE (iw, '(T2,A,T71,F10.1)') "ERI_GPW| Relative energy cutoff [Ry]", eri_rel_cutoff
1038 END IF
1039 !
1040 CALL calculate_eri_gpw(active_space_env%mos_active, active_space_env%active_orbitals, active_space_env%eri, qs_env, iw, &
1041 dft_control%roks)
1042 !
1043 CASE DEFAULT
1044 cpabort("Unknown ERI method")
1045 END SELECT
1046 IF (iw > 0) THEN
1047 DO isp = 1, SIZE(active_space_env%eri%eri)
1048 eri_mat => active_space_env%eri%eri(isp)%csr_mat
1049 nze_percentage = 100.0_dp*(real(eri_mat%nze_total, kind=dp) &
1050 /real(eri_mat%nrows_total, kind=dp))/real(eri_mat%ncols_total, kind=dp)
1051 WRITE (iw, '(/,T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1052 "Number of CSR non-zero elements:", eri_mat%nze_total
1053 WRITE (iw, '(T2,A,I2,T30,A,T68,F12.4)') "ERI_GPW| Spinmatrix:", isp, &
1054 "Percentage CSR non-zero elements:", nze_percentage
1055 WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1056 "nrows_total", eri_mat%nrows_total
1057 WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1058 "ncols_total", eri_mat%ncols_total
1059 WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1060 "nrows_local", eri_mat%nrows_local
1061 END DO
1062 CALL m_flush(iw)
1063 END IF
1064 CALL para_env%sync()
1065
1066 ! set the reference active space density matrix
1067 nspins = active_space_env%nspins
1068 ALLOCATE (active_space_env%p_active(nspins))
1069 DO isp = 1, nspins
1070 mo_set => active_space_env%mos_active(isp)
1071 CALL get_mo_set(mo_set, mo_coeff=mo_coeff, nmo=nmo)
1072 CALL create_subspace_matrix(mo_coeff, active_space_env%p_active(isp), nmo)
1073 END DO
1074 SELECT CASE (mselect)
1075 CASE DEFAULT
1076 cpabort("Unknown orbital selection method")
1077 CASE (casci_canonical, manual_selection)
1078 focc = 2.0_dp
1079 IF (nspins == 2) focc = 1.0_dp
1080 DO isp = 1, nspins
1081 fmat => active_space_env%p_active(isp)
1082 CALL cp_fm_set_all(fmat, alpha=0.0_dp)
1083 IF (nspins == 2) THEN
1084 IF (isp == 1) THEN
1085 nel = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2
1086 ELSE
1087 nel = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2
1088 END IF
1089 ELSE
1090 nel = active_space_env%nelec_active
1091 END IF
1092 DO i = 1, nmo_active
1093 m = active_space_env%active_orbitals(i, isp)
1094 fel = min(focc, real(nel, kind=dp))
1095 CALL cp_fm_set_element(fmat, m, m, fel)
1096 nel = nel - nint(fel)
1097 nel = max(nel, 0)
1098 END DO
1099 END DO
1100 CASE (wannier_projection)
1101 cpabort("NOT IMPLEMENTED")
1102 CASE (mao_projection)
1103 cpabort("NOT IMPLEMENTED")
1104 END SELECT
1105
1106 ! compute alpha-beta overlap matrix in case of spin-polarized calculation
1107 CALL calculate_spin_pol_overlap(active_space_env%mos_active, qs_env, active_space_env)
1108
1109 ! figure out if we have a new xc section for the AS
1110 xc_section => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE%XC")
1111 explicit = .false.
1112 IF (ASSOCIATED(xc_section)) CALL section_vals_get(xc_section, explicit=explicit)
1113
1114 ! rebuild KS matrix if needed
1115 IF (explicit) THEN
1116 ! release the hfx data if it was part of the SCF functional
1117 IF (ASSOCIATED(qs_env%x_data)) CALL hfx_release(qs_env%x_data)
1118 ! also release the admm environment in case we are using admm
1119 IF (ASSOCIATED(qs_env%admm_env)) CALL admm_env_release(qs_env%admm_env)
1120
1121 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1122 particle_set=particle_set, cell=cell, ks_env=ks_env)
1123 IF (dft_control%do_admm) THEN
1124 basis_type = 'AUX_FIT'
1125 ELSE
1126 basis_type = 'ORB'
1127 END IF
1128 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
1129 CALL hfx_create(qs_env%x_data, para_env, hfx_section, atomic_kind_set, &
1130 qs_kind_set, particle_set, dft_control, cell, orb_basis=basis_type, &
1131 nelectron_total=nelec_total)
1132
1133 qs_env%requires_matrix_vxc = .true. ! needs to be set only once
1134
1135 ! a bit of a hack: this forces a new re-init of HFX
1136 CALL set_ks_env(ks_env, s_mstruct_changed=.true.)
1137 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., &
1138 just_energy=.false., &
1139 ext_xc_section=xc_section)
1140 ! we need to reset it to false
1141 CALL set_ks_env(ks_env, s_mstruct_changed=.false.)
1142 ELSE
1143 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1144 END IF
1145 ! set the xc_section
1146 active_space_env%xc_section => xc_section
1147
1148 CALL get_qs_env(qs_env, energy=energy)
1149 ! transform KS/Fock, Vxc and Hcore to AS MO basis
1150 CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env)
1151 ! set the reference energy in the active space
1152 active_space_env%energy_ref = energy%total
1153 ! calculate inactive energy and embedding potential
1154 CALL subspace_fock_matrix(active_space_env, dft_control%roks)
1155
1156 ! associate the active space environment with the qs environment
1157 CALL set_qs_env(qs_env, active_space=active_space_env)
1158
1159 ! Perform the embedding calculation when an active-space solver is specified
1160 CALL section_vals_val_get(as_input, "AS_SOLVER", i_val=as_solver)
1161 SELECT CASE (as_solver)
1162 CASE (no_solver)
1163 IF (iw > 0) THEN
1164 WRITE (iw, '(/,T3,A)') "No active space solver specified, skipping embedding calculation"
1165 CALL m_flush(iw)
1166 END IF
1167 CALL para_env%sync()
1168 CASE (qiskit_solver)
1169 CALL rsdft_embedding(qs_env, active_space_env, as_input)
1170 CALL qs_scf_compute_properties(qs_env, wf_type="MC-DFT", do_mp2=.false.)
1171 CASE (fci_solver)
1172 CALL local_fci_embedding(qs_env, active_space_env, as_input)
1173 CALL qs_scf_compute_properties(qs_env, wf_type="MC-DFT", do_mp2=.false.)
1174 CASE DEFAULT
1175 cpabort("Unknown active space solver")
1176 END SELECT
1177
1178 ! Output a FCIDUMP file if requested
1179 IF (active_space_env%fcidump) CALL fcidump(active_space_env, as_input, dft_control%roks)
1180
1181 ! Output a QCSchema file if requested
1182 IF (active_space_env%qcschema) THEN
1183 CALL qcschema_env_create(qcschema_env, qs_env)
1184 CALL qcschema_to_hdf5(qcschema_env, active_space_env%qcschema_filename)
1185 CALL qcschema_env_release(qcschema_env)
1186 END IF
1187
1188 IF (iw > 0) THEN
1189 WRITE (iw, '(/,T2,A)') &
1190 '!-------------------- End of Active Space Interface --------------------------!'
1191 CALL m_flush(iw)
1192 END IF
1193 CALL para_env%sync()
1194
1195 CALL timestop(handle)
1196
1197 END SUBROUTINE active_space_main
1198
1199! **************************************************************************************************
1200!> \brief computes the alpha-beta overlap within the active subspace
1201!> \param mos the molecular orbital set within the active subspace
1202!> \param qs_env ...
1203!> \param active_space_env ...
1204!> \par History
1205!> 04.2016 created [JGH]
1206! **************************************************************************************************
1207 SUBROUTINE calculate_spin_pol_overlap(mos, qs_env, active_space_env)
1208
1209 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1210 TYPE(qs_environment_type), POINTER :: qs_env
1211 TYPE(active_space_type), POINTER :: active_space_env
1212
1213 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_spin_pol_overlap'
1214
1215 INTEGER :: handle, nmo, nspins
1216 LOGICAL :: do_kpoints
1217 TYPE(cp_fm_type), POINTER :: mo_coeff_a, mo_coeff_b
1218 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_matrix
1219 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: s_matrix_kp
1220
1221 CALL timeset(routinen, handle)
1222
1223 nspins = active_space_env%nspins
1224
1225 ! overlap in AO
1226 IF (nspins > 1) THEN
1227 CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
1228 IF (do_kpoints) THEN
1229 CALL get_qs_env(qs_env, matrix_s_kp=s_matrix_kp)
1230 s_matrix => s_matrix_kp(:, 1)
1231 ELSE
1232 CALL get_qs_env(qs_env, matrix_s=s_matrix)
1233 END IF
1234 ALLOCATE (active_space_env%sab_sub(1))
1235
1236 CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff_a, nmo=nmo)
1237 CALL get_mo_set(mo_set=mos(2), mo_coeff=mo_coeff_b, nmo=nmo)
1238 CALL subspace_operator(mo_coeff_a, nmo, s_matrix(1)%matrix, active_space_env%sab_sub(1), mo_coeff_b)
1239 END IF
1240
1241 CALL timestop(handle)
1242
1243 END SUBROUTINE calculate_spin_pol_overlap
1244
1245! **************************************************************************************************
1246!> \brief computes the one-electron operators in the subspace of the provided orbital set
1247!> \param mos the molecular orbital set within the active subspace
1248!> \param qs_env ...
1249!> \param active_space_env ...
1250!> \par History
1251!> 04.2016 created [JGH]
1252! **************************************************************************************************
1253 SUBROUTINE calculate_operators(mos, qs_env, active_space_env)
1254
1255 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1256 TYPE(qs_environment_type), POINTER :: qs_env
1257 TYPE(active_space_type), POINTER :: active_space_env
1258
1259 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_operators'
1260
1261 INTEGER :: handle, ispin, nmo, nspins
1262 TYPE(cp_fm_type), POINTER :: mo_coeff
1263 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: h_matrix, ks_matrix
1264
1265 CALL timeset(routinen, handle)
1266
1267 nspins = active_space_env%nspins
1268
1269 ! Kohn-Sham / Fock operator
1270 CALL cp_fm_release(active_space_env%ks_sub)
1271 CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix)
1272 ALLOCATE (active_space_env%ks_sub(nspins))
1273 DO ispin = 1, nspins
1274 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1275 CALL subspace_operator(mo_coeff, nmo, ks_matrix(ispin, 1)%matrix, active_space_env%ks_sub(ispin))
1276 END DO
1277
1278 ! Core Hamiltonian
1279 CALL cp_fm_release(active_space_env%h_sub)
1280
1281 NULLIFY (h_matrix)
1282 CALL get_qs_env(qs_env=qs_env, matrix_h_kp=h_matrix)
1283 ALLOCATE (active_space_env%h_sub(nspins))
1284 DO ispin = 1, nspins
1285 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1286 CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(ispin))
1287 END DO
1288
1289 CALL timestop(handle)
1290
1291 END SUBROUTINE calculate_operators
1292
1293! **************************************************************************************************
1294!> \brief computes a one-electron operator in the subspace of the provided orbital set
1295!> \param mo_coeff the orbital coefficient matrix
1296!> \param nmo the number of subspace orbitals
1297!> \param op_matrix operator matrix in AO basis
1298!> \param op_sub operator in orbital basis
1299!> \param mo_coeff_b the beta orbital coefficients
1300!> \par History
1301!> 04.2016 created [JGH]
1302! **************************************************************************************************
1303 SUBROUTINE subspace_operator(mo_coeff, nmo, op_matrix, op_sub, mo_coeff_b)
1304
1305 TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
1306 INTEGER, INTENT(IN) :: nmo
1307 TYPE(dbcsr_type), POINTER :: op_matrix
1308 TYPE(cp_fm_type), INTENT(INOUT) :: op_sub
1309 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: mo_coeff_b
1310
1311 CHARACTER(len=*), PARAMETER :: routinen = 'subspace_operator'
1312
1313 INTEGER :: handle, ncol, nrow
1314 TYPE(cp_fm_type) :: vectors
1315
1316 CALL timeset(routinen, handle)
1317
1318 CALL cp_fm_get_info(matrix=mo_coeff, ncol_global=ncol, nrow_global=nrow)
1319 cpassert(nmo <= ncol)
1320
1321 IF (nmo > 0) THEN
1322 CALL cp_fm_create(vectors, mo_coeff%matrix_struct, "vectors")
1323 CALL create_subspace_matrix(mo_coeff, op_sub, nmo)
1324
1325 IF (PRESENT(mo_coeff_b)) THEN
1326 ! if beta orbitals are present, compute the cross alpha_beta term
1327 CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff_b, vectors, nmo)
1328 ELSE
1329 ! otherwise the same spin, whatever that is
1330 CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff, vectors, nmo)
1331 END IF
1332
1333 CALL parallel_gemm('T', 'N', nmo, nmo, nrow, 1.0_dp, mo_coeff, vectors, 0.0_dp, op_sub)
1334 CALL cp_fm_release(vectors)
1335 END IF
1336
1337 CALL timestop(handle)
1338
1339 END SUBROUTINE subspace_operator
1340
1341! **************************************************************************************************
1342!> \brief creates a matrix of subspace size
1343!> \param orbitals the orbital coefficient matrix
1344!> \param op_sub operator in orbital basis
1345!> \param n the number of orbitals
1346!> \par History
1347!> 04.2016 created [JGH]
1348! **************************************************************************************************
1349 SUBROUTINE create_subspace_matrix(orbitals, op_sub, n)
1350
1351 TYPE(cp_fm_type), INTENT(IN) :: orbitals
1352 TYPE(cp_fm_type), INTENT(OUT) :: op_sub
1353 INTEGER, INTENT(IN) :: n
1354
1355 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1356
1357 IF (n > 0) THEN
1358
1359 NULLIFY (fm_struct)
1360 CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
1361 para_env=orbitals%matrix_struct%para_env, &
1362 context=orbitals%matrix_struct%context)
1363 CALL cp_fm_create(op_sub, fm_struct, name="Subspace operator")
1364 CALL cp_fm_struct_release(fm_struct)
1365
1366 END IF
1367
1368 END SUBROUTINE create_subspace_matrix
1369
1370! **************************************************************************************************
1371!> \brief computes the electron repulsion integrals using the GPW technology
1372!> \param mos the molecular orbital set within the active subspace
1373!> \param orbitals ...
1374!> \param eri_env ...
1375!> \param qs_env ...
1376!> \param iw ...
1377!> \param restricted ...
1378!> \par History
1379!> 04.2016 created [JGH]
1380! **************************************************************************************************
1381 SUBROUTINE calculate_eri_gpw(mos, orbitals, eri_env, qs_env, iw, restricted)
1382
1383 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1384 INTEGER, DIMENSION(:, :), POINTER :: orbitals
1385 TYPE(eri_type) :: eri_env
1386 TYPE(qs_environment_type), POINTER :: qs_env
1387 INTEGER, INTENT(IN) :: iw
1388 LOGICAL, INTENT(IN) :: restricted
1389
1390 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_eri_gpw'
1391
1392 INTEGER :: col_local, color, handle, i1, i2, i3, i4, i_multigrid, icount2, intcount, isp, &
1393 isp1, isp2, ispin, iwa1, iwa12, iwa2, iwb1, iwb12, iwb2, iwbs, iwbt, iwfn, n_multigrid, &
1394 ncol_global, ncol_local, nmm, nmo, nmo1, nmo2, nrow_global, nrow_local, nspins, &
1395 number_of_subgroups, nx, row_local, stored_integrals
1396 INTEGER, ALLOCATABLE, DIMENSION(:) :: eri_index
1397 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1398 LOGICAL :: print1, print2, &
1399 skip_load_balance_distributed
1400 REAL(kind=dp) :: dvol, erint, pair_int, &
1401 progression_factor, rc, rsize, t1, t2
1402 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eri
1403 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1404 TYPE(cell_type), POINTER :: cell
1405 TYPE(cp_blacs_env_type), POINTER :: blacs_env, blacs_env_sub
1406 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1407 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_matrix_pq_rnu, fm_matrix_pq_rs, &
1408 fm_mo_coeff_as
1409 TYPE(cp_fm_type), POINTER :: mo_coeff
1410 TYPE(dbcsr_p_type) :: mat_munu
1411 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_pq_rnu, mo_coeff_as
1412 TYPE(dft_control_type), POINTER :: dft_control
1413 TYPE(mp_para_env_type), POINTER :: para_env
1414 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1415 POINTER :: sab_orb_sub
1416 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1417 TYPE(pw_c1d_gs_type) :: pot_g, rho_g
1418 TYPE(pw_env_type), POINTER :: pw_env_sub
1419 TYPE(pw_poisson_type), POINTER :: poisson_env
1420 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1421 TYPE(pw_r3d_rs_type) :: rho_r, wfn_r
1422 TYPE(pw_r3d_rs_type), ALLOCATABLE, &
1423 DIMENSION(:, :), TARGET :: wfn_a
1424 TYPE(pw_r3d_rs_type), POINTER :: wfn1, wfn2, wfn3, wfn4
1425 TYPE(qs_control_type), POINTER :: qs_control, qs_control_old
1426 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1427 TYPE(qs_ks_env_type), POINTER :: ks_env
1428 TYPE(task_list_type), POINTER :: task_list_sub
1429
1430 CALL timeset(routinen, handle)
1431
1432 IF (iw > 0) t1 = m_walltime()
1433
1434 ! print levels
1435 SELECT CASE (eri_env%eri_gpw%print_level)
1436 CASE (silent_print_level)
1437 print1 = .false.
1438 print2 = .false.
1439 CASE (low_print_level)
1440 print1 = .false.
1441 print2 = .false.
1442 CASE (medium_print_level)
1443 print1 = .true.
1444 print2 = .false.
1445 CASE (high_print_level)
1446 print1 = .true.
1447 print2 = .true.
1448 CASE (debug_print_level)
1449 print1 = .true.
1450 print2 = .true.
1451 CASE DEFAULT
1452 ! do nothing
1453 END SELECT
1454
1455 ! Check the input group
1456 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1457 IF (eri_env%eri_gpw%group_size < 1) eri_env%eri_gpw%group_size = para_env%num_pe
1458 IF (mod(para_env%num_pe, eri_env%eri_gpw%group_size) /= 0) THEN
1459 cpabort("Group size must be a divisor of the total number of processes!")
1460 END IF
1461 ! Create a new para_env or reuse the old one
1462 IF (eri_env%eri_gpw%group_size == para_env%num_pe) THEN
1463 eri_env%para_env_sub => para_env
1464 CALL eri_env%para_env_sub%retain()
1465 blacs_env_sub => blacs_env
1466 CALL blacs_env_sub%retain()
1467 number_of_subgroups = 1
1468 color = 0
1469 ELSE
1470 number_of_subgroups = para_env%num_pe/eri_env%eri_gpw%group_size
1471 color = para_env%mepos/eri_env%eri_gpw%group_size
1472 ALLOCATE (eri_env%para_env_sub)
1473 CALL eri_env%para_env_sub%from_split(para_env, color)
1474 NULLIFY (blacs_env_sub)
1475 CALL cp_blacs_env_create(blacs_env_sub, eri_env%para_env_sub, blacs_grid_square, .true.)
1476 END IF
1477 CALL eri_env%comm_exchange%from_split(para_env, eri_env%para_env_sub%mepos)
1478
1479 ! This should be done differently! Copied from MP2 code
1480 CALL get_qs_env(qs_env, dft_control=dft_control)
1481 ALLOCATE (qs_control)
1482 qs_control_old => dft_control%qs_control
1483 qs_control = qs_control_old
1484 dft_control%qs_control => qs_control
1485 progression_factor = qs_control%progression_factor
1486 n_multigrid = SIZE(qs_control%e_cutoff)
1487 nspins = SIZE(mos)
1488 ! In case of ROHF/ROKS, we assume the orbital coefficients in both spin channels to be the same
1489 ! and save operations by calculating ERIs from only one spin channel
1490 IF (restricted) nspins = 1
1491 ! Allocate new cutoffs (just in private qs_control, not in qs_control_old)
1492 ALLOCATE (qs_control%e_cutoff(n_multigrid))
1493
1494 qs_control%cutoff = eri_env%eri_gpw%cutoff*0.5_dp
1495 qs_control%e_cutoff(1) = qs_control%cutoff
1496 DO i_multigrid = 2, n_multigrid
1497 qs_control%e_cutoff(i_multigrid) = qs_control%e_cutoff(i_multigrid - 1) &
1498 /progression_factor
1499 END DO
1500 qs_control%relative_cutoff = eri_env%eri_gpw%rel_cutoff*0.5_dp
1501
1502 ! For now, we will distribute neighbor lists etc. within the global communicator
1503 CALL get_qs_env(qs_env, ks_env=ks_env)
1504 CALL create_mat_munu(mat_munu, qs_env, eri_env%eri_gpw%eps_grid, blacs_env_sub, sab_orb_sub=sab_orb_sub, &
1505 do_alloc_blocks_from_nbl=.true., dbcsr_sym_type=dbcsr_type_symmetric)
1506 CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
1507
1508 ! Generate the appropriate pw_env
1509 NULLIFY (pw_env_sub)
1510 CALL pw_env_create(pw_env_sub)
1511 CALL pw_env_rebuild(pw_env_sub, qs_env, external_para_env=eri_env%para_env_sub)
1512 CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1513
1514 ! TODO: maybe we can let `pw_env_rebuild` do what we manually overwrite here?
1515 IF (eri_env%eri_gpw%redo_poisson) THEN
1516 ! We need to rebuild the Poisson solver on the fly
1517 IF (sum(eri_env%periodicity) /= 0) THEN
1518 poisson_env%parameters%solver = pw_poisson_periodic
1519 ELSE
1520 poisson_env%parameters%solver = pw_poisson_analytic
1521 END IF
1522 poisson_env%parameters%periodic = eri_env%periodicity
1523
1524 ! Rebuilds the poisson green (influence) function according
1525 ! to the poisson solver and parameters set so far.
1526 ! Also sets the variable poisson_env%rebuild to .FALSE.
1527 CALL pw_poisson_rebuild(poisson_env)
1528
1529 ! set the cutoff radius for the Greens function in case we use ANALYTIC Poisson solver
1530 CALL get_qs_env(qs_env, cell=cell)
1531 rc = cell%hmat(1, 1)
1532 DO iwa1 = 1, 3
1533 ! TODO: I think this is not the largest possible radius inscribed in the cell
1534 rc = min(rc, 0.5_dp*cell%hmat(iwa1, iwa1))
1535 END DO
1536 poisson_env%green_fft%radius = rc
1537
1538 ! Overwrite the Greens function with the one we want
1539 CALL pw_eri_green_create(poisson_env%green_fft, eri_env)
1540
1541 IF (iw > 0) THEN
1542 CALL get_qs_env(qs_env, cell=cell)
1543 IF (sum(cell%perd) /= sum(eri_env%periodicity)) THEN
1544 IF (sum(eri_env%periodicity) /= 0) THEN
1545 WRITE (unit=iw, fmt="(/,T2,A,T51,A30)") &
1546 "ERI_GPW| Switching Poisson solver to", "PERIODIC"
1547 ELSE
1548 WRITE (unit=iw, fmt="(/,T2,A,T51,A30)") &
1549 "ERI_GPW| Switching Poisson solver to", "ANALYTIC"
1550 END IF
1551 END IF
1552 ! print out the Greens function to check it matches the Poisson solver
1553 SELECT CASE (poisson_env%green_fft%method)
1554 CASE (periodic3d)
1555 WRITE (unit=iw, fmt="(T2,A,T51,A30)") &
1556 "ERI_GPW| Poisson Greens function", "PERIODIC"
1557 CASE (analytic0d)
1558 WRITE (unit=iw, fmt="(T2,A,T51,A30)") &
1559 "ERI_GPW| Poisson Greens function", "ANALYTIC"
1560 WRITE (unit=iw, fmt="(T2,A,T71,F10.4)") "ERI_GPW| Poisson cutoff radius", &
1561 poisson_env%green_fft%radius*angstrom
1562 CASE DEFAULT
1563 cpabort("Wrong Greens function setup")
1564 END SELECT
1565 END IF
1566 END IF
1567
1568 ALLOCATE (mo_coeff_as(nspins), fm_mo_coeff_as(nspins))
1569 DO ispin = 1, nspins
1570 block
1571 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: c, c_active
1572 INTEGER :: nmo
1573 TYPE(group_dist_d1_type) :: gd_array
1574 TYPE(cp_fm_type), POINTER :: mo_coeff
1575 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1576 CALL grep_rows_in_subgroups(para_env, eri_env%para_env_sub, mo_coeff, gd_array, c)
1577
1578 ALLOCATE (c_active(SIZE(c, 1), SIZE(orbitals, 1)))
1579 DO i1 = 1, SIZE(orbitals, 1)
1580 c_active(:, i1) = c(:, orbitals(i1, ispin))
1581 END DO
1582 CALL build_dbcsr_from_rows(eri_env%para_env_sub, mo_coeff_as(ispin), &
1583 c_active, mat_munu%matrix, gd_array, eri_env%eri_gpw%eps_filter)
1584 CALL release_group_dist(gd_array)
1585 DEALLOCATE (c, c_active)
1586 END block
1587
1588 CALL dbcsr_get_info(mo_coeff_as(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
1589
1590 NULLIFY (fm_struct)
1591 CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1592 nrow_global=nrow_global, ncol_global=ncol_global)
1593 CALL cp_fm_create(fm_mo_coeff_as(ispin), fm_struct)
1594 CALL cp_fm_struct_release(fm_struct)
1595
1596 CALL copy_dbcsr_to_fm(mo_coeff_as(ispin), fm_mo_coeff_as(ispin))
1597 END DO
1598
1599 IF (eri_env%method == eri_method_gpw_ht) THEN
1600 ! We need a task list
1601 NULLIFY (task_list_sub)
1602 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1603 CALL allocate_task_list(task_list_sub)
1604 CALL generate_qs_task_list(ks_env, task_list_sub, basis_type="ORB", &
1605 reorder_rs_grid_ranks=.true., &
1606 skip_load_balance_distributed=skip_load_balance_distributed, &
1607 pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
1608
1609 ! Create sparse matrices carrying the matrix products, Code borrowed from the MP2 GPW method
1610 ! Create equal distributions for them (no sparsity present)
1611 ! We use the routines from mp2 suggesting that one may replicate the grids later for better performance
1612 ALLOCATE (matrix_pq_rnu(nspins), fm_matrix_pq_rnu(nspins), fm_matrix_pq_rs(nspins))
1613 DO ispin = 1, nspins
1614 CALL dbcsr_create(matrix_pq_rnu(ispin), template=mo_coeff_as(ispin))
1615 CALL dbcsr_set(matrix_pq_rnu(ispin), 0.0_dp)
1616
1617 CALL dbcsr_get_info(matrix_pq_rnu(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
1618
1619 NULLIFY (fm_struct)
1620 CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1621 nrow_global=nrow_global, ncol_global=ncol_global)
1622 CALL cp_fm_create(fm_matrix_pq_rnu(ispin), fm_struct)
1623 CALL cp_fm_struct_release(fm_struct)
1624
1625 NULLIFY (fm_struct)
1626 CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1627 nrow_global=ncol_global, ncol_global=ncol_global)
1628 CALL cp_fm_create(fm_matrix_pq_rs(ispin), fm_struct)
1629 CALL cp_fm_struct_release(fm_struct)
1630 END DO
1631
1632 ! Copy the active space of the MOs into DBCSR matrices
1633 END IF
1634
1635 CALL auxbas_pw_pool%create_pw(wfn_r)
1636 CALL auxbas_pw_pool%create_pw(rho_g)
1637 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, &
1638 particle_set=particle_set, atomic_kind_set=atomic_kind_set)
1639
1640 ! pre-calculate wavefunctions on reals space grid
1641 nspins = SIZE(mos)
1642 ! In case of ROHF/ROKS, we assume the orbital coefficients in both spin channels to be the same
1643 ! and save operations by calculating ERIs from only one spin channel
1644 IF (restricted) nspins = 1
1645 IF (eri_env%eri_gpw%store_wfn) THEN
1646 ! pre-calculate wavefunctions on reals space grid
1647 rsize = 0.0_dp
1648 nmo = 0
1649 DO ispin = 1, nspins
1650 CALL get_mo_set(mo_set=mos(ispin), nmo=nx)
1651 nmo = max(nmo, nx)
1652 rsize = real(SIZE(wfn_r%array), kind=dp)*nx
1653 END DO
1654 IF (print1 .AND. iw > 0) THEN
1655 rsize = rsize*8._dp/1000000._dp
1656 WRITE (iw, "(T2,'ERI_GPW|',' Store active orbitals on real space grid ',T66,F12.3,' MB')") rsize
1657 END IF
1658 ALLOCATE (wfn_a(nmo, nspins))
1659 DO ispin = 1, nspins
1660 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1661 DO i1 = 1, SIZE(orbitals, 1)
1662 iwfn = orbitals(i1, ispin)
1663 CALL auxbas_pw_pool%create_pw(wfn_a(iwfn, ispin))
1664 CALL calculate_wavefunction(mo_coeff, iwfn, wfn_a(iwfn, ispin), rho_g, atomic_kind_set, &
1665 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1666 IF (print2 .AND. iw > 0) THEN
1667 WRITE (iw, "(T2,'ERI_GPW|',' Orbital stored ',I4,' Spin ',i1)") iwfn, ispin
1668 END IF
1669 END DO
1670 END DO
1671 ELSE
1672 ! Even if we do not store all WFNs, we still need containers for the functions to store
1673 ALLOCATE (wfn1, wfn2)
1674 CALL auxbas_pw_pool%create_pw(wfn1)
1675 CALL auxbas_pw_pool%create_pw(wfn2)
1676 IF (eri_env%method /= eri_method_gpw_ht) THEN
1677 ALLOCATE (wfn3, wfn4)
1678 CALL auxbas_pw_pool%create_pw(wfn3)
1679 CALL auxbas_pw_pool%create_pw(wfn4)
1680 END IF
1681 END IF
1682
1683 ! get some of the grids ready
1684 CALL auxbas_pw_pool%create_pw(rho_r)
1685 CALL auxbas_pw_pool%create_pw(pot_g)
1686
1687 ! run the FFT once, to set up buffers and to take into account the memory
1688 CALL pw_zero(rho_r)
1689 CALL pw_transfer(rho_r, rho_g)
1690 dvol = rho_r%pw_grid%dvol
1691
1692 IF (iw > 0) THEN
1693 CALL m_flush(iw)
1694 END IF
1695 ! calculate the integrals
1696 stored_integrals = 0
1697 DO isp1 = 1, nspins
1698 CALL get_mo_set(mo_set=mos(isp1), nmo=nmo1)
1699 nmm = (nmo1*(nmo1 + 1))/2
1700 DO i1 = 1, SIZE(orbitals, 1)
1701 iwa1 = orbitals(i1, isp1)
1702 IF (eri_env%eri_gpw%store_wfn) THEN
1703 wfn1 => wfn_a(iwa1, isp1)
1704 ELSE
1705 CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwa1, wfn1, rho_g, atomic_kind_set, &
1706 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1707 END IF
1708 DO i2 = i1, SIZE(orbitals, 1)
1709 iwa2 = orbitals(i2, isp1)
1710 iwa12 = csr_idx_to_combined(iwa1, iwa2, nmo1)
1711 ! Skip calculation directly if the pair is not part of our subgroup
1712 IF (mod(iwa12 - 1, eri_env%comm_exchange%num_pe) /= eri_env%comm_exchange%mepos) cycle
1713 iwa12 = (iwa12 - 1)/eri_env%comm_exchange%num_pe + 1
1714 IF (eri_env%eri_gpw%store_wfn) THEN
1715 wfn2 => wfn_a(iwa2, isp1)
1716 ELSE
1717 CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwa2, wfn2, rho_g, atomic_kind_set, &
1718 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1719 END IF
1720 ! calculate charge distribution and potential
1721 CALL pw_zero(rho_r)
1722 CALL pw_multiply(rho_r, wfn1, wfn2)
1723 CALL pw_transfer(rho_r, rho_g)
1724 CALL pw_poisson_solve(poisson_env, rho_g, pair_int, pot_g)
1725
1726 ! screening using pair_int
1727 IF (pair_int < eri_env%eps_integral) cycle
1728 CALL pw_transfer(pot_g, rho_r)
1729 !
1730 IF (eri_env%method == eri_method_gpw_ht) THEN
1731 CALL pw_scale(rho_r, dvol)
1732 DO isp2 = isp1, nspins
1733 CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2)
1734 nx = (nmo2*(nmo2 + 1))/2
1735 ALLOCATE (eri(nx), eri_index(nx))
1736 CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
1737 CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, &
1738 calculate_forces=.false., compute_tau=.false., gapw=.false., &
1739 pw_env_external=pw_env_sub, task_list_external=task_list_sub)
1740
1741 CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_as(isp2), &
1742 0.0_dp, matrix_pq_rnu(isp2), filter_eps=eri_env%eri_gpw%eps_filter)
1743 CALL copy_dbcsr_to_fm(matrix_pq_rnu(isp2), fm_matrix_pq_rnu(isp2))
1744
1745 CALL cp_fm_get_info(fm_matrix_pq_rnu(isp2), ncol_global=ncol_global, nrow_global=nrow_global)
1746
1747 CALL parallel_gemm("T", "N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
1748 fm_matrix_pq_rnu(isp2), fm_mo_coeff_as(isp2), &
1749 0.0_dp, fm_matrix_pq_rs(isp2))
1750 CALL parallel_gemm("T", "N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
1751 fm_mo_coeff_as(isp2), fm_matrix_pq_rnu(isp2), &
1752 1.0_dp, fm_matrix_pq_rs(isp2))
1753
1754 CALL cp_fm_get_info(fm_matrix_pq_rs(isp2), ncol_local=ncol_local, nrow_local=nrow_local, &
1755 col_indices=col_indices, row_indices=row_indices)
1756
1757 icount2 = 0
1758 DO col_local = 1, ncol_local
1759 iwb1 = orbitals(col_indices(col_local), isp2)
1760 IF (isp1 == isp2 .AND. iwb1 < iwa1) cycle
1761 DO row_local = 1, nrow_local
1762 iwb2 = orbitals(row_indices(row_local), isp2)
1763 IF (iwb2 < iwb1) cycle
1764 IF (isp1 == isp2 .AND. iwa1 == iwb1 .AND. iwb2 < iwa2) cycle
1765
1766 iwb12 = csr_idx_to_combined(iwb1, iwb2, nmo2)
1767 erint = fm_matrix_pq_rs(isp2)%local_data(row_local, col_local)
1768 IF (abs(erint) > eri_env%eps_integral) THEN
1769 icount2 = icount2 + 1
1770 eri(icount2) = erint
1771 eri_index(icount2) = iwb12
1772 END IF
1773 END DO
1774 END DO
1775 stored_integrals = stored_integrals + icount2
1776 !
1777 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1778 CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
1779 !
1780 DEALLOCATE (eri, eri_index)
1781 END DO
1782 ELSE IF (eri_env%method == eri_method_full_gpw) THEN
1783 DO isp2 = isp1, nspins
1784 CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2)
1785 nx = (nmo2*(nmo2 + 1))/2
1786 ALLOCATE (eri(nx), eri_index(nx))
1787 icount2 = 0
1788 iwbs = 1
1789 IF (isp1 == isp2) iwbs = i1
1790 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1791 DO i3 = iwbs, SIZE(orbitals, 1)
1792 iwb1 = orbitals(i3, isp2)
1793 IF (eri_env%eri_gpw%store_wfn) THEN
1794 wfn3 => wfn_a(iwb1, isp2)
1795 ELSE
1796 CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwb1, wfn3, rho_g, atomic_kind_set, &
1797 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1798 END IF
1799 CALL pw_zero(wfn_r)
1800 CALL pw_multiply(wfn_r, rho_r, wfn3)
1801 iwbt = i3
1802 IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
1803 DO i4 = iwbt, SIZE(orbitals, 1)
1804 iwb2 = orbitals(i4, isp2)
1805 IF (eri_env%eri_gpw%store_wfn) THEN
1806 wfn4 => wfn_a(iwb2, isp2)
1807 ELSE
1808 CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwb2, wfn4, rho_g, atomic_kind_set, &
1809 qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1810 END IF
1811 ! We reduce the amount of communication by collecting the local sums first and sum globally later
1812 erint = pw_integral_ab(wfn_r, wfn4, local_only=.true.)
1813 icount2 = icount2 + 1
1814 eri(icount2) = erint
1815 eri_index(icount2) = csr_idx_to_combined(iwb1, iwb2, nmo2)
1816 END DO
1817 END DO
1818 ! Now, we sum the integrals globally
1819 CALL eri_env%para_env_sub%sum(eri)
1820 ! and we reorder the integrals to prevent storing too small integrals
1821 intcount = 0
1822 icount2 = 0
1823 iwbs = 1
1824 IF (isp1 == isp2) iwbs = i1
1825 isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
1826 DO i3 = iwbs, SIZE(orbitals, 1)
1827 iwb1 = orbitals(i3, isp2)
1828 iwbt = i3
1829 IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
1830 DO i4 = iwbt, SIZE(orbitals, 1)
1831 iwb2 = orbitals(i4, isp2)
1832 intcount = intcount + 1
1833 erint = eri(intcount)
1834 IF (abs(erint) > eri_env%eps_integral) THEN
1835 IF (mod(intcount, eri_env%para_env_sub%num_pe) == eri_env%para_env_sub%mepos) THEN
1836 icount2 = icount2 + 1
1837 eri(icount2) = erint
1838 eri_index(icount2) = eri_index(intcount)
1839 END IF
1840 END IF
1841 END DO
1842 END DO
1843 stored_integrals = stored_integrals + icount2
1844 !
1845 CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
1846 !
1847 DEALLOCATE (eri, eri_index)
1848 END DO
1849 ELSE
1850 cpabort("Unknown option")
1851 END IF
1852 END DO
1853 END DO
1854 END DO
1855
1856 IF (print1 .AND. iw > 0) THEN
1857 WRITE (iw, "(T2,'ERI_GPW|',' Number of Integrals stored locally',T71,I10)") stored_integrals
1858 END IF
1859
1860 IF (eri_env%eri_gpw%store_wfn) THEN
1861 DO ispin = 1, nspins
1862 DO i1 = 1, SIZE(orbitals, 1)
1863 iwfn = orbitals(i1, ispin)
1864 CALL wfn_a(iwfn, ispin)%release()
1865 END DO
1866 END DO
1867 DEALLOCATE (wfn_a)
1868 ELSE
1869 CALL wfn1%release()
1870 CALL wfn2%release()
1871 DEALLOCATE (wfn1, wfn2)
1872 IF (eri_env%method /= eri_method_gpw_ht) THEN
1873 CALL wfn3%release()
1874 CALL wfn4%release()
1875 DEALLOCATE (wfn3, wfn4)
1876 END IF
1877 END IF
1878 CALL auxbas_pw_pool%give_back_pw(wfn_r)
1879 CALL auxbas_pw_pool%give_back_pw(rho_g)
1880 CALL auxbas_pw_pool%give_back_pw(rho_r)
1881 CALL auxbas_pw_pool%give_back_pw(pot_g)
1882
1883 IF (eri_env%method == eri_method_gpw_ht) THEN
1884 DO ispin = 1, nspins
1885 CALL dbcsr_release(mo_coeff_as(ispin))
1886 CALL dbcsr_release(matrix_pq_rnu(ispin))
1887 CALL cp_fm_release(fm_matrix_pq_rnu(ispin))
1888 CALL cp_fm_release(fm_matrix_pq_rs(ispin))
1889 END DO
1890 DEALLOCATE (matrix_pq_rnu, fm_matrix_pq_rnu, fm_matrix_pq_rs)
1891 CALL deallocate_task_list(task_list_sub)
1892 END IF
1893 DO ispin = 1, nspins
1894 CALL dbcsr_release(mo_coeff_as(ispin))
1895 CALL cp_fm_release(fm_mo_coeff_as(ispin))
1896 END DO
1897 DEALLOCATE (mo_coeff_as, fm_mo_coeff_as)
1898 CALL release_neighbor_list_sets(sab_orb_sub)
1899 CALL cp_blacs_env_release(blacs_env_sub)
1900 CALL dbcsr_release(mat_munu%matrix)
1901 DEALLOCATE (mat_munu%matrix)
1902 CALL pw_env_release(pw_env_sub)
1903 ! Return to the old qs_control
1904 dft_control%qs_control => qs_control_old
1905 DEALLOCATE (qs_control%e_cutoff)
1906 DEALLOCATE (qs_control)
1907
1908 ! print out progress
1909 IF (iw > 0) THEN
1910 t2 = m_walltime()
1911 WRITE (iw, '(/,T2,A,T66,F14.2)') "ERI_GPW| ERI calculation took (sec)", t2 - t1
1912 CALL m_flush(iw)
1913 END IF
1914
1915 CALL timestop(handle)
1916
1917 END SUBROUTINE calculate_eri_gpw
1918
1919! **************************************************************************************************
1920!> \brief Sets the Green's function for the ERI calculation. Here we deal with the G=0 case!
1921!> \param green ...
1922!> \param eri_env ...
1923!> \par History
1924!> 04.2016 created [JGH]
1925!> 08.2025 added support for the LR truncation [SB]
1926! **************************************************************************************************
1927 SUBROUTINE pw_eri_green_create(green, eri_env)
1928
1929 TYPE(greens_fn_type), INTENT(INOUT) :: green
1930 TYPE(eri_type) :: eri_env
1931
1932 COMPLEX(KIND=dp) :: erf_fac_p, z_p
1933 INTEGER :: ig
1934 REAL(kind=dp) :: cossin_fac, ea, erfcos_fac, exp_prefac, &
1935 g, g0, g2, g3d, ga, ginf, omega, &
1936 omega2, rc, rc2
1937
1938 ! initialize influence function
1939 associate(gf => green%influence_fn, grid => green%influence_fn%pw_grid)
1940 SELECT CASE (green%method)
1941 CASE (periodic3d)
1942
1943 SELECT CASE (eri_env%operator)
1944 CASE (eri_operator_coulomb)
1945 DO ig = grid%first_gne0, grid%ngpts_cut_local
1946 g2 = grid%gsq(ig)
1947 gf%array(ig) = fourpi/g2
1948 END DO
1949 IF (grid%have_g0) gf%array(1) = 0.0_dp
1950
1951 CASE (eri_operator_yukawa)
1952 CALL cp_warn(__location__, "Yukawa operator has not been tested")
1953 omega2 = eri_env%omega**2
1954 DO ig = grid%first_gne0, grid%ngpts_cut_local
1955 g2 = grid%gsq(ig)
1956 gf%array(ig) = fourpi/(omega2 + g2)
1957 END DO
1958 IF (grid%have_g0) gf%array(1) = fourpi/omega2
1959
1960 CASE (eri_operator_erf)
1961 omega2 = eri_env%omega**2
1962 DO ig = grid%first_gne0, grid%ngpts_cut_local
1963 g2 = grid%gsq(ig)
1964 gf%array(ig) = fourpi/g2*exp(-0.25_dp*g2/omega2)
1965 END DO
1966 IF (grid%have_g0) gf%array(1) = 0.0_dp
1967
1968 CASE (eri_operator_erfc)
1969 omega2 = eri_env%omega**2
1970 DO ig = grid%first_gne0, grid%ngpts_cut_local
1971 g2 = grid%gsq(ig)
1972 gf%array(ig) = fourpi/g2*(1.0_dp - exp(-0.25_dp*g2/omega2))
1973 END DO
1974 IF (grid%have_g0) gf%array(1) = pi/omega2
1975
1976 CASE (eri_operator_trunc)
1977 rc = eri_env%cutoff_radius
1978 DO ig = grid%first_gne0, grid%ngpts_cut_local
1979 g2 = grid%gsq(ig)
1980 g = sqrt(g2)
1981 ! Taylor expansion around zero
1982 IF (g*rc >= 0.005_dp) THEN
1983 gf%array(ig) = fourpi/g2*(1.0_dp - cos(g*rc))
1984 ELSE
1985 gf%array(ig) = fourpi/g2*(g*rc)**2/2.0_dp*(1.0_dp - (g*rc)**2/12.0_dp)
1986 END IF
1987 END DO
1988 IF (grid%have_g0) gf%array(1) = twopi*rc**2
1989
1990 CASE (eri_operator_lr_trunc)
1991 omega = eri_env%omega
1992 omega2 = omega**2
1993 rc = eri_env%cutoff_radius
1994 rc2 = rc**2
1995 g0 = 0.001_dp ! threshold for the G=0 case
1996 ginf = 20.0_dp ! threshold for the Taylor exapnsion arounf G=∞
1997 DO ig = grid%first_gne0, grid%ngpts_cut_local
1998 g2 = grid%gsq(ig)
1999 g = sqrt(g2)
2000 IF (g <= 2.0_dp*g0) THEN
2001 gf%array(ig) = -pi/omega2*erf(omega*rc) &
2002 + twopi*rc2*erf(omega*rc) &
2003 + 2*rootpi*rc*exp(-omega2*rc2)/omega
2004 ELSE IF (g >= 2.0_dp*ginf*omega) THEN
2005 ! exponential prefactor
2006 exp_prefac = exp(-omega2*rc2)/(rootpi*(omega2*rc2 + 0.25_dp*g2/omega2))
2007 ! cos sin factor
2008 cossin_fac = omega*rc*cos(g*rc) - 0.5_dp*g/omega*sin(g*rc)
2009 ! real erf term with cosine
2010 erfcos_fac = erf(omega*rc)*cos(g*rc)
2011 ! Combine terms
2012 gf%array(ig) = fourpi/g2*(-exp_prefac*cossin_fac - erfcos_fac)
2013 ELSE
2014 ! exponential prefactor
2015 exp_prefac = twopi/g2*exp(-0.25_dp*g2/omega2)
2016 ! Compute complex arguments for erf
2017 z_p = cmplx(omega*rc, 0.5_dp*g/omega, kind=dp)
2018 ! Evaluate complex error functions
2019 erf_fac_p = 2.0_dp*real(erfz_fast(z_p))
2020 ! Real erf term with cosine
2021 erfcos_fac = fourpi/g2*erf(omega*rc)*cos(g*rc)
2022 ! Combine terms
2023 gf%array(ig) = exp_prefac*erf_fac_p - erfcos_fac
2024 END IF
2025 END DO
2026 IF (grid%have_g0) THEN
2027 gf%array(1) = -pi/omega2*erf(omega*rc) &
2028 + twopi*rc2*erf(omega*rc) &
2029 + 2*rootpi*rc*exp(-omega2*rc2)/omega
2030 END IF
2031
2032 CASE DEFAULT
2033 cpabort("Please specify a valid operator for the periodic Poisson solver")
2034 END SELECT
2035
2036 ! The analytic Poisson solver simply limits the domain of integration
2037 ! of the Fourier transform to a sphere of radius Rc, rather than integrating
2038 ! over all space (-∞,∞)
2039 CASE (analytic0d)
2040
2041 SELECT CASE (eri_env%operator)
2042 ! This is identical to the truncated Coulomb operator integrated
2043 ! over all space, when the truncation radius is equal to the radius of
2044 ! the Poisson solver
2045 CASE (eri_operator_coulomb, eri_operator_trunc)
2046 IF (eri_env%operator == eri_operator_coulomb) THEN
2047 rc = green%radius
2048 ELSE
2049 rc = eri_env%cutoff_radius
2050 END IF
2051 DO ig = grid%first_gne0, grid%ngpts_cut_local
2052 g2 = grid%gsq(ig)
2053 g = sqrt(g2)
2054 ! Taylor expansion around zero
2055 IF (g*rc >= 0.005_dp) THEN
2056 gf%array(ig) = fourpi/g2*(1.0_dp - cos(g*rc))
2057 ELSE
2058 gf%array(ig) = fourpi/g2*(g*rc)**2/2.0_dp*(1.0_dp - (g*rc)**2/12.0_dp)
2059 END IF
2060 END DO
2061 IF (grid%have_g0) gf%array(1) = twopi*rc**2
2062
2063 ! Not tested
2064 CASE (eri_operator_yukawa)
2065 CALL cp_warn(__location__, "Yukawa operator has not been tested")
2066 rc = green%radius
2067 omega = eri_env%omega
2068 ea = exp(-omega*rc)
2069 DO ig = grid%first_gne0, grid%ngpts_cut_local
2070 g2 = grid%gsq(ig)
2071 g = sqrt(g2)
2072 g3d = fourpi/(omega**2 + g2)
2073 gf%array(ig) = g3d*(1.0_dp - ea*(cos(g*rc) + omega/g*sin(g*rc)))
2074 END DO
2075 IF (grid%have_g0) gf%array(1) = fourpi/(omega**2)*(1.0_dp - ea*(1.0_dp + omega*rc))
2076
2077 ! Long-range Coulomb
2078 ! TODO: this should be equivalent to LR truncated Coulomb from above!
2079 CASE (eri_operator_erf, eri_operator_lr_trunc)
2080 IF (eri_env%operator == eri_operator_erf) THEN
2081 rc = green%radius
2082 ELSE
2083 rc = eri_env%cutoff_radius
2084 END IF
2085 omega2 = eri_env%omega**2
2086 DO ig = grid%first_gne0, grid%ngpts_cut_local
2087 g2 = grid%gsq(ig)
2088 g = sqrt(g2)
2089 ga = -0.25_dp*g2/omega2
2090 gf%array(ig) = fourpi/g2*exp(ga)*(1.0_dp - cos(g*rc))
2091 END DO
2092 IF (grid%have_g0) gf%array(1) = twopi*rc**2
2093
2094 ! Short-range Coulomb
2095 ! TODO: this should actually be properly derived and see whether it is correct
2096 CASE (eri_operator_erfc)
2097 CALL cp_warn(__location__, &
2098 "Short-range Coulomb operator may be incorrect with ANALYTIC0D Poisson solver")
2099 rc = green%radius
2100 omega2 = eri_env%omega**2
2101 DO ig = grid%first_gne0, grid%ngpts_cut_local
2102 g2 = grid%gsq(ig)
2103 g = sqrt(g2)
2104 ga = -0.25_dp*g2/omega2
2105 gf%array(ig) = fourpi/g2*(1.0_dp - exp(ga))*(1.0_dp - cos(g*rc))
2106 END DO
2107 IF (grid%have_g0) gf%array(1) = pi/omega2
2108
2109 CASE DEFAULT
2110 cpabort("Unsupported operator")
2111 END SELECT
2112
2113 CASE DEFAULT
2114 cpabort("Unsupported Poisson solver")
2115 END SELECT
2116 END associate
2117
2118 END SUBROUTINE pw_eri_green_create
2119
2120! **************************************************************************************************
2121!> \brief Adds data for a new row to the csr matrix
2122!> \param csr_mat ...
2123!> \param nnz ...
2124!> \param rdat ...
2125!> \param rind ...
2126!> \param irow ...
2127!> \par History
2128!> 04.2016 created [JGH]
2129! **************************************************************************************************
2130 SUBROUTINE update_csr_matrix(csr_mat, nnz, rdat, rind, irow)
2131
2132 TYPE(dbcsr_csr_type), INTENT(INOUT) :: csr_mat
2133 INTEGER, INTENT(IN) :: nnz
2134 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rdat
2135 INTEGER, DIMENSION(:), INTENT(IN) :: rind
2136 INTEGER, INTENT(IN) :: irow
2137
2138 INTEGER :: k, nrow, nze, nze_new
2139
2140 IF (irow /= 0) THEN
2141 nze = csr_mat%nze_local
2142 nze_new = nze + nnz
2143 ! values
2144 CALL reallocate(csr_mat%nzval_local%r_dp, 1, nze_new)
2145 csr_mat%nzval_local%r_dp(nze + 1:nze_new) = rdat(1:nnz)
2146 ! col indices
2147 CALL reallocate(csr_mat%colind_local, 1, nze_new)
2148 csr_mat%colind_local(nze + 1:nze_new) = rind(1:nnz)
2149 ! rows
2150 nrow = csr_mat%nrows_local
2151 CALL reallocate(csr_mat%rowptr_local, 1, irow + 1)
2152 csr_mat%rowptr_local(nrow + 1:irow) = nze + 1
2153 csr_mat%rowptr_local(irow + 1) = nze_new + 1
2154 ! nzerow
2155 CALL reallocate(csr_mat%nzerow_local, 1, irow)
2156 DO k = nrow + 1, irow
2157 csr_mat%nzerow_local(k) = csr_mat%rowptr_local(k + 1) - csr_mat%rowptr_local(k)
2158 END DO
2159 csr_mat%nrows_local = irow
2160 csr_mat%nze_local = csr_mat%nze_local + nnz
2161 END IF
2162 csr_mat%nze_total = csr_mat%nze_total + nnz
2163 csr_mat%has_indices = .true.
2164
2165 END SUBROUTINE update_csr_matrix
2166
2167! **************************************************************************************************
2168!> \brief Computes and prints the active orbitals on Cube Files
2169!> \param input ...
2170!> \param qs_env the qs_env in which the qs_env lives
2171!> \param mos ...
2172! **************************************************************************************************
2173 SUBROUTINE print_orbital_cubes(input, qs_env, mos)
2174 TYPE(section_vals_type), POINTER :: input
2175 TYPE(qs_environment_type), POINTER :: qs_env
2176 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2177
2178 CHARACTER(LEN=default_path_length) :: filebody, filename, title
2179 INTEGER :: i, imo, isp, nmo, str(3), unit_nr
2180 INTEGER, DIMENSION(:), POINTER :: alist, blist, istride
2181 LOGICAL :: do_mo, explicit_a, explicit_b
2182 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2183 TYPE(cell_type), POINTER :: cell
2184 TYPE(cp_fm_type), POINTER :: mo_coeff
2185 TYPE(dft_control_type), POINTER :: dft_control
2186 TYPE(mp_para_env_type), POINTER :: para_env
2187 TYPE(particle_list_type), POINTER :: particles
2188 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2189 TYPE(pw_c1d_gs_type) :: wf_g
2190 TYPE(pw_env_type), POINTER :: pw_env
2191 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2192 TYPE(pw_r3d_rs_type) :: wf_r
2193 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2194 TYPE(qs_subsys_type), POINTER :: subsys
2195 TYPE(section_vals_type), POINTER :: dft_section, scf_input
2196
2197 CALL section_vals_val_get(input, "FILENAME", c_val=filebody)
2198 CALL section_vals_val_get(input, "STRIDE", i_vals=istride)
2199 IF (SIZE(istride) == 1) THEN
2200 str(1:3) = istride(1)
2201 ELSE IF (SIZE(istride) == 3) THEN
2202 str(1:3) = istride(1:3)
2203 ELSE
2204 cpabort("STRIDE arguments inconsistent")
2205 END IF
2206 CALL section_vals_val_get(input, "ALIST", i_vals=alist, explicit=explicit_a)
2207 CALL section_vals_val_get(input, "BLIST", i_vals=blist, explicit=explicit_b)
2208
2209 CALL get_qs_env(qs_env=qs_env, &
2210 dft_control=dft_control, &
2211 para_env=para_env, &
2212 subsys=subsys, &
2213 atomic_kind_set=atomic_kind_set, &
2214 qs_kind_set=qs_kind_set, &
2215 cell=cell, &
2216 particle_set=particle_set, &
2217 pw_env=pw_env, &
2218 input=scf_input)
2219
2220 CALL qs_subsys_get(subsys, particles=particles)
2221 !
2222 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2223 CALL auxbas_pw_pool%create_pw(wf_r)
2224 CALL auxbas_pw_pool%create_pw(wf_g)
2225 !
2226 dft_section => section_vals_get_subs_vals(scf_input, "DFT")
2227 !
2228 DO isp = 1, SIZE(mos)
2229 CALL get_mo_set(mo_set=mos(isp), mo_coeff=mo_coeff, nmo=nmo)
2230
2231 IF (SIZE(mos) > 1) THEN
2232 SELECT CASE (isp)
2233 CASE (1)
2234 CALL write_mo_set_to_output_unit(mos(isp), qs_kind_set, particle_set, &
2235 dft_section, 4, 0, final_mos=.true., spin="ALPHA")
2236 CASE (2)
2237 CALL write_mo_set_to_output_unit(mos(isp), qs_kind_set, particle_set, &
2238 dft_section, 4, 0, final_mos=.true., spin="BETA")
2239 CASE DEFAULT
2240 cpabort("Invalid spin")
2241 END SELECT
2242 ELSE
2243 CALL write_mo_set_to_output_unit(mos(isp), qs_kind_set, particle_set, &
2244 dft_section, 4, 0, final_mos=.true.)
2245 END IF
2246
2247 DO imo = 1, nmo
2248 IF (isp == 1 .AND. explicit_a) THEN
2249 IF (alist(1) == -1) THEN
2250 do_mo = .true.
2251 ELSE
2252 do_mo = .false.
2253 DO i = 1, SIZE(alist)
2254 IF (imo == alist(i)) do_mo = .true.
2255 END DO
2256 END IF
2257 ELSE IF (isp == 2 .AND. explicit_b) THEN
2258 IF (blist(1) == -1) THEN
2259 do_mo = .true.
2260 ELSE
2261 do_mo = .false.
2262 DO i = 1, SIZE(blist)
2263 IF (imo == blist(i)) do_mo = .true.
2264 END DO
2265 END IF
2266 ELSE
2267 do_mo = .true.
2268 END IF
2269 IF (.NOT. do_mo) cycle
2270 CALL calculate_wavefunction(mo_coeff, imo, wf_r, wf_g, atomic_kind_set, &
2271 qs_kind_set, cell, dft_control, particle_set, pw_env)
2272 IF (para_env%is_source()) THEN
2273 WRITE (filename, '(A,A1,I4.4,A1,I1.1,A)') trim(filebody), "_", imo, "_", isp, ".cube"
2274 CALL open_file(filename, unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE")
2275 WRITE (title, *) "Active Orbital ", imo, " spin ", isp
2276 ELSE
2277 unit_nr = -1
2278 END IF
2279 CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=istride)
2280 IF (para_env%is_source()) THEN
2281 CALL close_file(unit_nr)
2282 END IF
2283 END DO
2284 END DO
2285
2286 CALL auxbas_pw_pool%give_back_pw(wf_r)
2287 CALL auxbas_pw_pool%give_back_pw(wf_g)
2288
2289 END SUBROUTINE print_orbital_cubes
2290
2291! **************************************************************************************************
2292!> \brief Writes a FCIDUMP file
2293!> \param active_space_env ...
2294!> \param as_input ...
2295!> \param restricted ...
2296!> \par History
2297!> 04.2016 created [JGH]
2298! **************************************************************************************************
2299 SUBROUTINE fcidump(active_space_env, as_input, restricted)
2300
2301 TYPE(active_space_type), POINTER :: active_space_env
2302 TYPE(section_vals_type), POINTER :: as_input
2303 LOGICAL, INTENT(IN) :: restricted
2304
2305 INTEGER :: i, i1, i2, i3, i4, isym, iw, m1, m2, &
2306 ms2, nmo, norb, nspins
2307 REAL(kind=dp) :: checksum, esub
2308 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: fmat
2309 TYPE(cp_logger_type), POINTER :: logger
2310 TYPE(eri_fcidump_checksum) :: eri_checksum
2311
2312 checksum = 0.0_dp
2313
2314 logger => cp_get_default_logger()
2315 iw = cp_print_key_unit_nr(logger, as_input, "FCIDUMP", &
2316 extension=".fcidump", file_status="REPLACE", file_action="WRITE", file_form="FORMATTED")
2317 !
2318 nspins = active_space_env%nspins
2319 norb = SIZE(active_space_env%active_orbitals, 1)
2320 ms2 = active_space_env%multiplicity - 1
2321 IF (nspins == 1 .OR. restricted) THEN
2322 ! Closed shell or restricted open-shell
2323 associate(nelec => active_space_env%nelec_active)
2324
2325 IF (iw > 0) THEN
2326 WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
2327 isym = 1
2328 WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb)
2329 isym = 0
2330 WRITE (iw, "(A,I1,A)") " ISYM=", isym, ","
2331 IF (restricted) WRITE (iw, "(A,I1,A)") " UHF=", 0, ","
2332 WRITE (iw, "(A)") " /"
2333 END IF
2334 !
2335 ! Print integrals: ERI
2336 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
2337 eri_fcidump_print(iw, 1, 1), 1, 1)
2338 CALL eri_checksum%set(1, 1)
2339 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
2340
2341 ! Print integrals: Fij
2342 ! replicate Fock matrix
2343 nmo = active_space_env%eri%norb
2344 ALLOCATE (fmat(nmo, nmo))
2345 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2346 IF (iw > 0) THEN
2347 i3 = 0; i4 = 0
2348 DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
2349 i1 = active_space_env%active_orbitals(m1, 1)
2350 DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
2351 i2 = active_space_env%active_orbitals(m2, 1)
2352 checksum = checksum + abs(fmat(i1, i2))
2353 WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2354 END DO
2355 END DO
2356 END IF
2357 DEALLOCATE (fmat)
2358 ! Print energy
2359 esub = active_space_env%energy_inactive
2360 i1 = 0; i2 = 0; i3 = 0; i4 = 0
2361 checksum = checksum + abs(esub)
2362 IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4
2363 END associate
2364
2365 ELSE
2366 associate(nelec => active_space_env%nelec_active)
2367
2368 IF (iw > 0) THEN
2369 WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
2370 isym = 1
2371 WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb)
2372 isym = 0
2373 WRITE (iw, "(A,I1,A)") " ISYM=", isym, ","
2374 WRITE (iw, "(A,I1,A)") " UHF=", 1, ","
2375 WRITE (iw, "(A)") " /"
2376 END IF
2377 !
2378 ! Print integrals: ERI
2379 ! alpha-alpha
2380 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
2381 eri_fcidump_print(iw, 1, 1), 1, 1)
2382 CALL eri_checksum%set(1, 1)
2383 CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
2384 ! alpha-beta
2385 CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, &
2386 eri_fcidump_print(iw, 1, norb + 1), 1, 2)
2387 CALL eri_checksum%set(1, norb + 1)
2388 CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, eri_checksum, 1, 2)
2389 ! beta-beta
2390 CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, &
2391 eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2)
2392 CALL eri_checksum%set(norb + 1, norb + 1)
2393 CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2)
2394 ! Print integrals: Fij
2395 ! alpha
2396 nmo = active_space_env%eri%norb
2397 ALLOCATE (fmat(nmo, nmo))
2398 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2399 IF (iw > 0) THEN
2400 i3 = 0; i4 = 0
2401 DO m1 = 1, norb
2402 i1 = active_space_env%active_orbitals(m1, 1)
2403 DO m2 = m1, norb
2404 i2 = active_space_env%active_orbitals(m2, 1)
2405 checksum = checksum + abs(fmat(i1, i2))
2406 WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2407 END DO
2408 END DO
2409 END IF
2410 DEALLOCATE (fmat)
2411 ! beta
2412 ALLOCATE (fmat(nmo, nmo))
2413 CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat)
2414 IF (iw > 0) THEN
2415 i3 = 0; i4 = 0
2416 DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
2417 i1 = active_space_env%active_orbitals(m1, 2)
2418 DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
2419 i2 = active_space_env%active_orbitals(m2, 2)
2420 checksum = checksum + abs(fmat(i1, i2))
2421 WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4
2422 END DO
2423 END DO
2424 END IF
2425 DEALLOCATE (fmat)
2426 ! Print energy
2427 esub = active_space_env%energy_inactive
2428 i1 = 0; i2 = 0; i3 = 0; i4 = 0
2429 checksum = checksum + abs(esub)
2430 IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4
2431 END associate
2432 END IF
2433 !
2434 CALL cp_print_key_finished_output(iw, logger, as_input, "FCIDUMP")
2435
2436 !>>
2437 iw = cp_logger_get_default_io_unit(logger)
2438 IF (iw > 0) WRITE (iw, '(T4,A,T66,F12.8)') "FCIDUMP| Checksum:", eri_checksum%checksum + checksum
2439 !<<
2440
2441 END SUBROUTINE fcidump
2442
2443! **************************************************************************************************
2444!> \brief replicate and symmetrize a matrix
2445!> \param norb the number of orbitals
2446!> \param distributed_matrix ...
2447!> \param replicated_matrix ...
2448! **************************************************************************************************
2449 SUBROUTINE replicate_and_symmetrize_matrix(norb, distributed_matrix, replicated_matrix)
2450 INTEGER, INTENT(IN) :: norb
2451 TYPE(cp_fm_type), INTENT(IN) :: distributed_matrix
2452 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: replicated_matrix
2453
2454 INTEGER :: i1, i2
2455 REAL(dp) :: mval
2456
2457 replicated_matrix(:, :) = 0.0_dp
2458 DO i1 = 1, norb
2459 DO i2 = i1, norb
2460 CALL cp_fm_get_element(distributed_matrix, i1, i2, mval)
2461 replicated_matrix(i1, i2) = mval
2462 replicated_matrix(i2, i1) = mval
2463 END DO
2464 END DO
2465 END SUBROUTINE replicate_and_symmetrize_matrix
2466
2467! **************************************************************************************************
2468!> \brief Calculates active space Fock matrix and inactive energy
2469!> \param active_space_env ...
2470!> \param restricted ...
2471!> \par History
2472!> 06.2016 created [JGH]
2473! **************************************************************************************************
2474 SUBROUTINE subspace_fock_matrix(active_space_env, restricted)
2475
2476 TYPE(active_space_type), POINTER :: active_space_env
2477 LOGICAL, INTENT(IN) :: restricted
2478
2479 INTEGER :: i1, i2, is, norb, nspins
2480 REAL(kind=dp) :: eeri, eref, esub, mval
2481 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ks_a_mat, ks_a_ref, ks_b_mat, ks_b_ref, &
2482 ks_mat, ks_ref, p_a_mat, p_b_mat, p_mat
2483 TYPE(cp_fm_type), POINTER :: matrix, mo_coef
2484 TYPE(dbcsr_csr_type), POINTER :: eri, eri_aa, eri_ab, eri_bb
2485
2486 eref = active_space_env%energy_ref
2487 nspins = active_space_env%nspins
2488
2489 IF (nspins == 1) THEN
2490 CALL get_mo_set(active_space_env%mos_active(1), nmo=norb, mo_coeff=mo_coef)
2491 !
2492 ! Loop over ERI, calculate subspace HF energy and Fock matrix
2493 !
2494 ! replicate KS, Core, and P matrices
2495 ALLOCATE (ks_mat(norb, norb), ks_ref(norb, norb), p_mat(norb, norb))
2496 ks_ref = 0.0_dp
2497
2498 ! ks_mat contains the KS/Fock matrix (of full density) projected onto the AS MO subspace (f_ref in eq. 19)
2499 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_mat)
2500 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_mat)
2501
2502 ! compute ks_ref = V_H[rho^A] + V_HFX[rho^A]
2503 eri => active_space_env%eri%eri(1)%csr_mat
2504 CALL build_subspace_fock_matrix(active_space_env%active_orbitals, eri, p_mat, ks_ref, &
2505 active_space_env%eri%comm_exchange)
2506
2507 ! compute eeri = E_H[rho^A] + E_HFX[rho^A] as
2508 ! eeri = 1/2 * (SUM_pq (V_H[rho^A] + V_HFX[rho^A])_pq * D^A_pq)
2509 eeri = 0.5_dp*sum(ks_ref*p_mat)
2510
2511 ! now calculate the inactive energy acoording to eq. 19, that is
2512 ! esub = E^I = E_ref - f_ref .* D^A + E_H[rho^A] + E_HFX[rho^A]
2513 ! where f^ref = ks_mat, which is the KS/Fock matrix in MO basis, transformed previously
2514 ! and is equal to ks_mat = h^0 + V_core + V_H[rho] + V_HFX[rho]
2515 esub = eref - sum(ks_mat(1:norb, 1:norb)*p_mat(1:norb, 1:norb)) + eeri
2516
2517 ! reuse ks_mat to store f^I = f^ref - (V_H[rho^A] + V_HFX[rho^A]) according to eq. 20
2518 ks_mat(1:norb, 1:norb) = ks_mat(1:norb, 1:norb) - ks_ref(1:norb, 1:norb)
2519 ! this is now the embedding potential for the AS calculation!
2520
2521 active_space_env%energy_inactive = esub
2522
2523 CALL cp_fm_release(active_space_env%fock_sub)
2524 ALLOCATE (active_space_env%fock_sub(nspins))
2525 DO is = 1, nspins
2526 matrix => active_space_env%ks_sub(is)
2527 CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2528 name="Active Fock operator")
2529 END DO
2530 matrix => active_space_env%fock_sub(1)
2531 DO i1 = 1, norb
2532 DO i2 = 1, norb
2533 mval = ks_mat(i1, i2)
2534 CALL cp_fm_set_element(matrix, i1, i2, mval)
2535 END DO
2536 END DO
2537 ELSE
2538
2539 CALL get_mo_set(active_space_env%mos_active(1), nmo=norb)
2540 !
2541 ! Loop over ERI, calculate subspace HF energy and Fock matrix
2542 !
2543 ! replicate KS, Core, and P matrices
2544 ALLOCATE (ks_a_mat(norb, norb), ks_b_mat(norb, norb), &
2545 & ks_a_ref(norb, norb), ks_b_ref(norb, norb), &
2546 & p_a_mat(norb, norb), p_b_mat(norb, norb))
2547 ks_a_ref(:, :) = 0.0_dp; ks_b_ref(:, :) = 0.0_dp
2548
2549 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_a_mat)
2550 CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(2), p_b_mat)
2551 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_a_mat)
2552 CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(2), ks_b_mat)
2553 !
2554 !
2555 IF (restricted) THEN
2556 ! In the restricted case, we use the same ERIs for each spin channel
2557 eri_aa => active_space_env%eri%eri(1)%csr_mat
2558 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_aa, eri_aa, p_a_mat, p_b_mat, ks_a_ref, &
2559 tr_mixed_eri=.false., comm_exchange=active_space_env%eri%comm_exchange)
2560 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_aa, eri_aa, p_b_mat, p_a_mat, ks_b_ref, &
2561 tr_mixed_eri=.true., comm_exchange=active_space_env%eri%comm_exchange)
2562 ELSE
2563 eri_aa => active_space_env%eri%eri(1)%csr_mat
2564 eri_ab => active_space_env%eri%eri(2)%csr_mat
2565 eri_bb => active_space_env%eri%eri(3)%csr_mat
2566 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, &
2567 tr_mixed_eri=.false., comm_exchange=active_space_env%eri%comm_exchange)
2568 CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_bb, eri_ab, p_b_mat, p_a_mat, ks_b_ref, &
2569 tr_mixed_eri=.true., comm_exchange=active_space_env%eri%comm_exchange)
2570 END IF
2571 !
2572 ! calculate energy
2573 eeri = 0.0_dp
2574 eeri = 0.5_dp*(sum(ks_a_ref*p_a_mat) + sum(ks_b_ref*p_b_mat))
2575 esub = eref - sum(ks_a_mat*p_a_mat) - sum(ks_b_mat*p_b_mat) + eeri
2576 ks_a_mat(:, :) = ks_a_mat(:, :) - ks_a_ref(:, :)
2577 ks_b_mat(:, :) = ks_b_mat(:, :) - ks_b_ref(:, :)
2578 !
2579 active_space_env%energy_inactive = esub
2580 !
2581 CALL cp_fm_release(active_space_env%fock_sub)
2582 ALLOCATE (active_space_env%fock_sub(nspins))
2583 DO is = 1, nspins
2584 matrix => active_space_env%ks_sub(is)
2585 CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2586 name="Active Fock operator")
2587 END DO
2588
2589 matrix => active_space_env%fock_sub(1)
2590 DO i1 = 1, norb
2591 DO i2 = 1, norb
2592 mval = ks_a_mat(i1, i2)
2593 CALL cp_fm_set_element(matrix, i1, i2, mval)
2594 END DO
2595 END DO
2596 matrix => active_space_env%fock_sub(2)
2597 DO i1 = 1, norb
2598 DO i2 = 1, norb
2599 mval = ks_b_mat(i1, i2)
2600 CALL cp_fm_set_element(matrix, i1, i2, mval)
2601 END DO
2602 END DO
2603
2604 END IF
2605
2606 END SUBROUTINE subspace_fock_matrix
2607
2608! **************************************************************************************************
2609!> \brief build subspace fockian
2610!> \param active_orbitals the active orbital indices
2611!> \param eri two electon integrals in MO
2612!> \param p_mat density matrix
2613!> \param ks_ref fockian matrix
2614!> \param comm_exchange ...
2615! **************************************************************************************************
2616 SUBROUTINE build_subspace_fock_matrix(active_orbitals, eri, p_mat, ks_ref, comm_exchange)
2617 INTEGER, DIMENSION(:, :), INTENT(IN) :: active_orbitals
2618 TYPE(dbcsr_csr_type), INTENT(IN) :: eri
2619 REAL(dp), DIMENSION(:, :), INTENT(IN) :: p_mat
2620 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: ks_ref
2621 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange
2622
2623 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_subspace_fock_matrix'
2624
2625 INTEGER :: handle, i1, i12, i12l, i2, i3, i34, &
2626 i34l, i4, irptr, m1, m2, nindex, &
2627 nmo_total, norb
2628 REAL(dp) :: erint
2629 TYPE(mp_comm_type) :: mp_group
2630
2631 CALL timeset(routinen, handle)
2632
2633 ! Nothing to do
2634 norb = SIZE(active_orbitals, 1)
2635 nmo_total = SIZE(p_mat, 1)
2636 nindex = (nmo_total*(nmo_total + 1))/2
2637 CALL mp_group%set_handle(eri%mp_group%get_handle())
2638 DO m1 = 1, norb
2639 i1 = active_orbitals(m1, 1)
2640 DO m2 = m1, norb
2641 i2 = active_orbitals(m2, 1)
2642 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2643 IF (mod(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
2644 i12l = (i12 - 1)/comm_exchange%num_pe + 1
2645 irptr = eri%rowptr_local(i12l) - 1
2646 DO i34l = 1, eri%nzerow_local(i12l)
2647 i34 = eri%colind_local(irptr + i34l)
2648 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2649 erint = eri%nzval_local%r_dp(irptr + i34l)
2650 ! Coulomb
2651 ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2652 IF (i3 /= i4) THEN
2653 ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2654 END IF
2655 IF (i12 /= i34) THEN
2656 ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2657 IF (i1 /= i2) THEN
2658 ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2659 END IF
2660 END IF
2661 ! Exchange
2662 erint = -0.5_dp*erint
2663 ks_ref(i1, i3) = ks_ref(i1, i3) + erint*p_mat(i2, i4)
2664 IF (i1 /= i2) THEN
2665 ks_ref(i2, i3) = ks_ref(i2, i3) + erint*p_mat(i1, i4)
2666 END IF
2667 IF (i3 /= i4) THEN
2668 ks_ref(i1, i4) = ks_ref(i1, i4) + erint*p_mat(i2, i3)
2669 END IF
2670 IF (i1 /= i2 .AND. i3 /= i4) THEN
2671 ks_ref(i2, i4) = ks_ref(i2, i4) + erint*p_mat(i1, i3)
2672 END IF
2673 END DO
2674 END IF
2675 END DO
2676 END DO
2677 !
2678 DO m1 = 1, norb
2679 i1 = active_orbitals(m1, 1)
2680 DO m2 = m1, norb
2681 i2 = active_orbitals(m2, 1)
2682 ks_ref(i2, i1) = ks_ref(i1, i2)
2683 END DO
2684 END DO
2685 CALL mp_group%sum(ks_ref)
2686
2687 CALL timestop(handle)
2688
2689 END SUBROUTINE build_subspace_fock_matrix
2690
2691! **************************************************************************************************
2692!> \brief build subspace fockian for unrestricted spins
2693!> \param active_orbitals the active orbital indices
2694!> \param eri_aa two electon integrals in MO with parallel spins
2695!> \param eri_ab two electon integrals in MO with anti-parallel spins
2696!> \param p_a_mat density matrix for up-spin
2697!> \param p_b_mat density matrix for down-spin
2698!> \param ks_a_ref fockian matrix for up-spin
2699!> \param tr_mixed_eri boolean to indicate Coulomb interaction alignment
2700!> \param comm_exchange ...
2701! **************************************************************************************************
2702 SUBROUTINE build_subspace_spin_fock_matrix(active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, tr_mixed_eri, &
2703 comm_exchange)
2704 INTEGER, DIMENSION(:, :), INTENT(IN) :: active_orbitals
2705 TYPE(dbcsr_csr_type), INTENT(IN) :: eri_aa, eri_ab
2706 REAL(dp), DIMENSION(:, :), INTENT(IN) :: p_a_mat, p_b_mat
2707 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: ks_a_ref
2708 LOGICAL, INTENT(IN) :: tr_mixed_eri
2709 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange
2710
2711 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_subspace_spin_fock_matrix'
2712
2713 INTEGER :: handle, i1, i12, i12l, i2, i3, i34, &
2714 i34l, i4, irptr, m1, m2, nindex, &
2715 nmo_total, norb, spin1, spin2
2716 REAL(dp) :: erint
2717 TYPE(mp_comm_type) :: mp_group
2718
2719 CALL timeset(routinen, handle)
2720
2721 norb = SIZE(active_orbitals, 1)
2722 nmo_total = SIZE(p_a_mat, 1)
2723 nindex = (nmo_total*(nmo_total + 1))/2
2724 IF (tr_mixed_eri) THEN
2725 spin1 = 2
2726 spin2 = 1
2727 ELSE
2728 spin1 = 1
2729 spin2 = 2
2730 END IF
2731 DO m1 = 1, norb
2732 i1 = active_orbitals(m1, spin1)
2733 DO m2 = m1, norb
2734 i2 = active_orbitals(m2, spin1)
2735 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2736 IF (mod(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
2737 i12l = (i12 - 1)/comm_exchange%num_pe + 1
2738 irptr = eri_aa%rowptr_local(i12l) - 1
2739 DO i34l = 1, eri_aa%nzerow_local(i12l)
2740 i34 = eri_aa%colind_local(irptr + i34l)
2741 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2742 erint = eri_aa%nzval_local%r_dp(irptr + i34l)
2743 ! Coulomb
2744 !F_ij += (ij|kl)*d_kl
2745 ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_a_mat(i3, i4)
2746 IF (i12 /= i34) THEN
2747 !F_kl += (ij|kl)*d_ij
2748 ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_a_mat(i1, i2)
2749 END IF
2750 ! Exchange
2751 erint = -1.0_dp*erint
2752 !F_ik -= (ij|kl)*d_jl
2753 ks_a_ref(i1, i3) = ks_a_ref(i1, i3) + erint*p_a_mat(i2, i4)
2754 IF (i1 /= i2) THEN
2755 !F_jk -= (ij|kl)*d_il
2756 ks_a_ref(i2, i3) = ks_a_ref(i2, i3) + erint*p_a_mat(i1, i4)
2757 END IF
2758 IF (i3 /= i4) THEN
2759 !F_il -= (ij|kl)*d_jk
2760 ks_a_ref(i1, i4) = ks_a_ref(i1, i4) + erint*p_a_mat(i2, i3)
2761 END IF
2762 IF (i1 /= i2 .AND. i3 /= i4) THEN
2763 !F_jl -= (ij|kl)*d_ik
2764 ks_a_ref(i2, i4) = ks_a_ref(i2, i4) + erint*p_a_mat(i1, i3)
2765 END IF
2766 END DO
2767 END IF
2768 END DO
2769 END DO
2770 !
2771
2772 DO m1 = 1, norb
2773 i1 = active_orbitals(m1, 1)
2774 DO m2 = m1, norb
2775 i2 = active_orbitals(m2, 1)
2776 i12 = csr_idx_to_combined(i1, i2, nmo_total)
2777 IF (mod(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
2778 i12l = (i12 - 1)/comm_exchange%num_pe + 1
2779 irptr = eri_ab%rowptr_local(i12l) - 1
2780 DO i34l = 1, eri_ab%nzerow_local(i12l)
2781 i34 = eri_ab%colind_local(irptr + i34l)
2782 CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2783 erint = eri_ab%nzval_local%r_dp(irptr + i34l)
2784 ! Coulomb
2785 IF (tr_mixed_eri) THEN
2786 !F_kl += (kl beta|ij alpha )*d_alpha_ij
2787 ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_b_mat(i1, i2)
2788 ELSE
2789 !F_ij += (ij alpha|kl beta )*d_beta_kl
2790 ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_b_mat(i3, i4)
2791 END IF
2792 END DO
2793 END IF
2794 END DO
2795 END DO
2796 !
2797 DO m1 = 1, norb
2798 i1 = active_orbitals(m1, spin1)
2799 DO m2 = m1, norb
2800 i2 = active_orbitals(m2, spin1)
2801 ks_a_ref(i2, i1) = ks_a_ref(i1, i2)
2802 END DO
2803 END DO
2804 CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
2805 CALL mp_group%sum(ks_a_ref)
2806
2807 CALL timestop(handle)
2808
2809 END SUBROUTINE build_subspace_spin_fock_matrix
2810
2811! **************************************************************************************************
2812!> \brief Creates a local basis
2813!> \param pro_basis_set ...
2814!> \param zval ...
2815!> \param ishell ...
2816!> \param nshell ...
2817!> \param lnam ...
2818!> \par History
2819!> 05.2016 created [JGH]
2820! **************************************************************************************************
2821 SUBROUTINE create_pro_basis(pro_basis_set, zval, ishell, nshell, lnam)
2822 TYPE(gto_basis_set_type), POINTER :: pro_basis_set
2823 INTEGER, INTENT(IN) :: zval, ishell
2824 INTEGER, DIMENSION(:), INTENT(IN) :: nshell
2825 CHARACTER(len=*), DIMENSION(:), INTENT(IN) :: lnam
2826
2827 CHARACTER(len=6), DIMENSION(:), POINTER :: sym
2828 INTEGER :: i, l, nj
2829 INTEGER, DIMENSION(4, 7) :: ne
2830 INTEGER, DIMENSION(:), POINTER :: lq, nq
2831 REAL(kind=dp), DIMENSION(:), POINTER :: zet
2832 TYPE(sto_basis_set_type), POINTER :: sto_basis_set
2833
2834 cpassert(.NOT. ASSOCIATED(pro_basis_set))
2835 NULLIFY (sto_basis_set)
2836
2837 ! electronic configuration
2838 ne = 0
2839 DO l = 1, 4 !lq(1)+1
2840 nj = 2*(l - 1) + 1
2841 DO i = l, 7 ! nq(1)
2842 ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
2843 ne(l, i) = max(ne(l, i), 0)
2844 ne(l, i) = min(ne(l, i), 2*nj)
2845 END DO
2846 END DO
2847 ALLOCATE (nq(ishell), lq(ishell), zet(ishell), sym(ishell))
2848 DO i = 1, ishell
2849 nq(i) = nshell(i)
2850 SELECT CASE (lnam(i))
2851 CASE ('S', 's')
2852 lq(i) = 0
2853 CASE ('P', 'p')
2854 lq(i) = 1
2855 CASE ('D', 'd')
2856 lq(i) = 2
2857 CASE ('F', 'f')
2858 lq(i) = 3
2859 CASE DEFAULT
2860 cpabort("Wrong l QN")
2861 END SELECT
2862 sym(i) = lnam(i)
2863 zet(i) = srules(zval, ne, nq(1), lq(1))
2864 END DO
2865 CALL allocate_sto_basis_set(sto_basis_set)
2866 CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zet, symbol=sym)
2867 CALL create_gto_from_sto_basis(sto_basis_set, pro_basis_set, 6)
2868 pro_basis_set%norm_type = 2
2869 CALL init_orb_basis_set(pro_basis_set)
2870 CALL deallocate_sto_basis_set(sto_basis_set)
2871
2872 END SUBROUTINE create_pro_basis
2873
2874! **************************************************************************************************
2875!> \brief Update the density matrix in AO basis with the active density contribution
2876!> \param active_space_env the active space environment
2877!> \param rho_ao the density matrix in AO basis
2878! **************************************************************************************************
2879 SUBROUTINE update_density_ao(active_space_env, rho_ao)
2880 TYPE(active_space_type), POINTER :: active_space_env
2881 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2882
2883 INTEGER :: ispin, nao, nmo, nspins
2884 TYPE(cp_fm_type) :: r, u
2885 TYPE(cp_fm_type), POINTER :: c_active, p_active_mo
2886 TYPE(dbcsr_type), POINTER :: p_inactive_ao
2887 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active
2888
2889 ! Transform the AS density matrix P_MO to the atomic orbital basis,
2890 ! this is simply C * P_MO * C^T
2891 nspins = active_space_env%nspins
2892 mos_active => active_space_env%mos_active
2893 DO ispin = 1, nspins
2894 ! size of p_inactive_ao is (nao x nao)
2895 p_inactive_ao => active_space_env%pmat_inactive(ispin)%matrix
2896
2897 ! copy p_inactive_ao to rho_ao
2898 CALL dbcsr_copy(rho_ao(ispin)%matrix, p_inactive_ao)
2899
2900 ! size of p_active_mo is (nmo x nmo)
2901 p_active_mo => active_space_env%p_active(ispin)
2902
2903 ! calculate R = p_mo
2904 CALL cp_fm_create(r, p_active_mo%matrix_struct)
2905 CALL cp_fm_to_fm(p_active_mo, r)
2906
2907 ! calculate U = C * p_mo
2908 CALL get_mo_set(mos_active(ispin), mo_coeff=c_active, nao=nao, nmo=nmo)
2909 CALL cp_fm_create(u, c_active%matrix_struct)
2910 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, c_active, r, 0.0_dp, u)
2911
2912 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(ispin)%matrix, &
2913 matrix_v=u, matrix_g=c_active, ncol=nmo, alpha=1.0_dp)
2914
2915 CALL cp_fm_release(r)
2916 CALL cp_fm_release(u)
2917 END DO
2918
2919 END SUBROUTINE update_density_ao
2920
2921! **************************************************************************************************
2922!> \brief Print each value on the master node
2923!> \param this object reference
2924!> \param i i-index
2925!> \param j j-index
2926!> \param k k-index
2927!> \param l l-index
2928!> \param val value of the integral at (i,j,k.l)
2929!> \return always true to dump all integrals
2930! **************************************************************************************************
2931 LOGICAL FUNCTION eri_fcidump_print_func(this, i, j, k, l, val) RESULT(cont)
2932 CLASS(eri_fcidump_print), INTENT(inout) :: this
2933 INTEGER, INTENT(in) :: i, j, k, l
2934 REAL(kind=dp), INTENT(in) :: val
2935
2936 ! write to the actual file only on the master
2937 IF (this%unit_nr > 0) THEN
2938 WRITE (this%unit_nr, "(ES23.16,4I4)") val, i + this%bra_start - 1, j + this%bra_start - 1, &
2939 & k + this%ket_start - 1, l + this%ket_start - 1
2940 END IF
2941
2942 cont = .true.
2943 END FUNCTION eri_fcidump_print_func
2944
2945! **************************************************************************************************
2946!> \brief checksum each value on the master node
2947!> \param this object reference
2948!> \param i i-index
2949!> \param j j-index
2950!> \param k k-index
2951!> \param l l-index
2952!> \param val value of the integral at (i,j,k.l)
2953!> \return always true to dump all integrals
2954! **************************************************************************************************
2955 LOGICAL FUNCTION eri_fcidump_checksum_func(this, i, j, k, l, val) RESULT(cont)
2956 CLASS(eri_fcidump_checksum), INTENT(inout) :: this
2957 INTEGER, INTENT(in) :: i, j, k, l
2958 REAL(kind=dp), INTENT(in) :: val
2959 mark_used(i)
2960 mark_used(j)
2961 mark_used(k)
2962 mark_used(l)
2963
2964 this%checksum = this%checksum + abs(val)
2965
2966 cont = .true.
2967 END FUNCTION eri_fcidump_checksum_func
2968
2969! **************************************************************************************************
2970!> \brief Compute and print the AS rdm and the natural orbitals occupation numbers
2971!> \param active_space_env active space environment
2972!> \param iw output unit
2973!> \author Stefano Battaglia
2974! **************************************************************************************************
2975 SUBROUTINE print_pmat_noon(active_space_env, iw)
2976 TYPE(active_space_type), POINTER :: active_space_env
2977 INTEGER :: iw
2978
2979 INTEGER :: i1, i2, ii, ispin, jm, m1, m2, &
2980 nmo_active, nspins
2981 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: noon, pmat
2982 TYPE(cp_fm_type), POINTER :: p_active
2983
2984 nspins = active_space_env%nspins
2985 nmo_active = active_space_env%nmo_active
2986
2987 ALLOCATE (noon(nmo_active, nspins))
2988 ALLOCATE (pmat(nmo_active, nmo_active))
2989
2990 DO ispin = 1, nspins
2991 p_active => active_space_env%p_active(ispin)
2992 noon(:, ispin) = 0.0_dp
2993 pmat = 0.0_dp
2994
2995 DO i1 = 1, nmo_active
2996 m1 = active_space_env%active_orbitals(i1, ispin)
2997 DO i2 = 1, nmo_active
2998 m2 = active_space_env%active_orbitals(i2, ispin)
2999 CALL cp_fm_get_element(p_active, m1, m2, pmat(i1, i2))
3000 END DO
3001 END DO
3002
3003 IF (iw > 0) THEN
3004 WRITE (iw, '(/,T3,A,I2,A)') "Active space density matrix for spin ", ispin
3005 DO i1 = 1, nmo_active
3006 DO ii = 1, nmo_active, 8
3007 jm = min(7, nmo_active - ii)
3008 WRITE (iw, '(T3,6(F9.4))') (pmat(i1, ii + i2), i2=0, jm)
3009 END DO
3010 END DO
3011 END IF
3012
3013 ! diagonalize the density matrix
3014 CALL diamat_all(pmat, noon(:, ispin))
3015
3016 IF (iw > 0) THEN
3017 WRITE (iw, '(/,T3,A,I2,A)') "Natural orbitals occupation numbers for spin ", ispin
3018 DO i1 = 1, nmo_active, 8
3019 jm = min(7, nmo_active - i1)
3020 ! noons are stored in ascending order, so reverse-print them
3021 WRITE (iw, '(T3,6(F9.4))') (noon(nmo_active - i1 - i2 + 1, ispin), i2=0, jm)
3022 END DO
3023 END IF
3024
3025 END DO
3026
3027 DEALLOCATE (noon)
3028 DEALLOCATE (pmat)
3029
3030 END SUBROUTINE print_pmat_noon
3031
3032! **************************************************************************************************
3033!> \brief Run range-separated DFT embedding with the local FCI active-space solver.
3034!> \param qs_env Quickstep environment
3035!> \param active_space_env active-space environment
3036!> \param as_input ACTIVE_SPACE input section
3037! **************************************************************************************************
3038 SUBROUTINE local_fci_embedding(qs_env, active_space_env, as_input)
3039 TYPE(qs_environment_type), POINTER :: qs_env
3040 TYPE(active_space_type), POINTER :: active_space_env
3041 TYPE(section_vals_type), POINTER :: as_input
3042
3043 CHARACTER(len=*), PARAMETER :: routinen = 'local_fci_embedding'
3044
3045 INTEGER :: handle, iter, iw, max_iter
3046 LOGICAL :: converged, do_scf_embedding
3047 REAL(kind=dp) :: delta_e, energy_corr, energy_new, &
3048 energy_old, energy_scf, eps_iter, t1, &
3049 t2
3050 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: p_act_mo_a, p_act_mo_b
3051 TYPE(cp_logger_type), POINTER :: logger
3052 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
3053 TYPE(dft_control_type), POINTER :: dft_control
3054 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active
3055 TYPE(mp_para_env_type), POINTER :: para_env
3056 TYPE(qs_energy_type), POINTER :: energy
3057 TYPE(qs_ks_env_type), POINTER :: ks_env
3058 TYPE(qs_rho_type), POINTER :: rho
3059
3060 CALL timeset(routinen, handle)
3061
3062 t1 = m_walltime()
3063 logger => cp_get_default_logger()
3064 iw = cp_logger_get_default_io_unit(logger)
3065
3066 CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
3067
3068 CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
3069 active_space_env%do_scf_embedding = do_scf_embedding
3070 CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter)
3071 IF (max_iter < 0) cpabort("Specify a non-negative number of max iterations.")
3072 CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
3073 IF (eps_iter < 0.0) cpabort("Specify a non-negative convergence threshold.")
3074
3075 CALL get_qs_env(qs_env, rho=rho, energy=energy, ks_env=ks_env)
3076 CALL qs_rho_get(rho, rho_ao=rho_ao)
3077
3078 IF (iw > 0) THEN
3079 WRITE (unit=iw, fmt="(/,T2,A,/)") &
3080 "RANGE-SEPARATED DFT EMBEDDING WITH LIBFCI SOLVER"
3081 WRITE (iw, '(T3,A,T68,I12)') "Max. iterations", max_iter
3082 WRITE (iw, '(T3,A,T68,E12.4)') "Conv. threshold", eps_iter
3083 WRITE (iw, '(T3,A,T68,A)') "Density mixer", trim(active_space_mixing_label(active_space_env))
3084 WRITE (iw, '(T3,A,T66,F14.2)') "Mixing alpha", active_space_env%alpha
3085 WRITE (unit=iw, fmt="(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") &
3086 "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", repeat("-", 78)
3087 END IF
3088
3089 iter = 0
3090 converged = .false.
3091 energy_scf = active_space_env%energy_ref
3092 energy_new = energy_scf
3093 mos_active => active_space_env%mos_active
3094
3095 DO WHILE (iter < max_iter)
3096 iter = iter + 1
3097
3098 IF (active_space_env%nspins == 2) THEN
3099 CALL solve_active_space_fci(active_space_env, para_env, p_act_mo_a, p_act_mo_b)
3100 active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3101 CALL update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
3102 DEALLOCATE (p_act_mo_a, p_act_mo_b)
3103 ELSE
3104 CALL solve_active_space_fci(active_space_env, para_env, p_act_mo_a)
3105 active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3106 CALL update_active_density(p_act_mo_a, active_space_env)
3107 DEALLOCATE (p_act_mo_a)
3108 END IF
3109
3110 energy_old = energy_new
3111 energy_new = active_space_env%energy_total
3112 energy_corr = energy_new - energy_scf
3113 delta_e = energy_new - energy_old
3114
3115 t2 = t1
3116 t1 = m_walltime()
3117 IF (iw > 0) THEN
3118 WRITE (unit=iw, &
3119 fmt="(T3,I4,T11,A,T19,F6.1,T28,F18.10,T49,F18.10,T70,ES11.2)") &
3120 iter, trim(active_space_mixing_label(active_space_env)), &
3121 t1 - t2, energy_corr, energy_new, delta_e
3122 CALL m_flush(iw)
3123 END IF
3124
3125 CALL update_density_ao(active_space_env, rho_ao)
3126 CALL qs_rho_update_rho(rho, qs_env=qs_env)
3127 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
3128 CALL evaluate_core_matrix_traces(qs_env)
3129 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., &
3130 just_energy=.false., &
3131 ext_xc_section=active_space_env%xc_section)
3132
3133 active_space_env%energy_ref = energy%total
3134 CALL calculate_operators(mos_active, qs_env, active_space_env)
3135 CALL subspace_fock_matrix(active_space_env, dft_control%roks)
3136
3137 IF (.NOT. active_space_env%do_scf_embedding) THEN
3138 IF (iw > 0) THEN
3139 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3140 "*** one-shot embedding correction finished ***"
3141 END IF
3142 converged = .true.
3143 EXIT
3144 ELSE IF (abs(delta_e) <= eps_iter) THEN
3145 IF (iw > 0) THEN
3146 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3147 "*** rs-DFT embedding run converged in ", iter, " iteration(s) ***"
3148 END IF
3149 converged = .true.
3150 EXIT
3151 END IF
3152 END DO
3153
3154 IF (.NOT. converged) THEN
3155 IF (iw > 0) THEN
3156 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3157 "*** rs-DFT embedding did not converged after ", iter, " iteration(s) ***"
3158 END IF
3159 END IF
3160
3161 energy%total = active_space_env%energy_total
3162
3163 IF (iw > 0) THEN
3164 WRITE (unit=iw, fmt="(/,T3,A)") &
3165 "Final energy contributions:"
3166 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3167 "Inactive energy:", active_space_env%energy_inactive
3168 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3169 "Active energy:", active_space_env%energy_active
3170 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3171 "Correlation energy:", energy_corr
3172 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3173 "Total rs-DFT energy:", active_space_env%energy_total
3174 END IF
3175
3176 CALL print_pmat_noon(active_space_env, iw)
3177 CALL para_env%sync()
3178 CALL timestop(handle)
3179
3180 END SUBROUTINE local_fci_embedding
3181
3182! **************************************************************************************************
3183!> \brief ...
3184!> \param qs_env ...
3185!> \param active_space_env ...
3186!> \param as_input ...
3187! **************************************************************************************************
3188 SUBROUTINE rsdft_embedding(qs_env, active_space_env, as_input)
3189 TYPE(qs_environment_type), POINTER :: qs_env
3190 TYPE(active_space_type), POINTER :: active_space_env
3191 TYPE(section_vals_type), POINTER :: as_input
3192
3193 CHARACTER(len=*), PARAMETER :: routinen = 'rsdft_embedding'
3194 INTEGER :: handle
3195
3196#ifdef __NO_SOCKETS
3197 CALL timeset(routinen, handle)
3198 cpabort("CP2K was compiled with the __NO_SOCKETS option!")
3199 mark_used(qs_env)
3200 mark_used(active_space_env)
3201 mark_used(as_input)
3202#else
3203
3204 INTEGER :: iw, client_fd, socket_fd, iter, max_iter
3205 LOGICAL :: converged, do_scf_embedding, ionode
3206 REAL(kind=dp) :: delta_e, energy_corr, energy_new, &
3207 energy_old, energy_scf, eps_iter, t1, t2
3208 TYPE(cp_logger_type), POINTER :: logger
3209 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
3210 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active
3211 TYPE(mp_para_env_type), POINTER :: para_env
3212 TYPE(qs_energy_type), POINTER :: energy
3213 TYPE(qs_ks_env_type), POINTER :: ks_env
3214 TYPE(qs_rho_type), POINTER :: rho
3215 TYPE(dft_control_type), POINTER :: dft_control
3216
3217 CALL timeset(routinen, handle)
3218
3219 t1 = m_walltime()
3220
3221 logger => cp_get_default_logger()
3222 iw = cp_logger_get_default_io_unit(logger)
3223
3224 CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
3225 ionode = para_env%is_source()
3226
3227 ! get info from the input
3228 CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
3229 active_space_env%do_scf_embedding = do_scf_embedding
3230 CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter)
3231 IF (max_iter < 0) cpabort("Specify a non-negative number of max iterations.")
3232 CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
3233 IF (eps_iter < 0.0) cpabort("Specify a non-negative convergence threshold.")
3234
3235 ! create the socket and wait for the client to connect
3236 CALL initialize_socket(socket_fd, client_fd, as_input, ionode)
3237 CALL para_env%sync()
3238
3239 ! send two-electron integrals to the client
3240 CALL send_eri_to_client(client_fd, active_space_env, para_env)
3241
3242 ! get pointer to density in ao basis
3243 CALL get_qs_env(qs_env, rho=rho, energy=energy, ks_env=ks_env)
3244 CALL qs_rho_get(rho, rho_ao=rho_ao)
3245
3246 IF (iw > 0) THEN
3247 WRITE (unit=iw, fmt="(/,T2,A,/)") &
3248 "RANGE-SEPARATED DFT EMBEDDING SELF-CONSISTENT OPTIMIZATION"
3249
3250 WRITE (iw, '(T3,A,T68,I12)') "Max. iterations", max_iter
3251 WRITE (iw, '(T3,A,T68,E12.4)') "Conv. threshold", eps_iter
3252 WRITE (iw, '(T3,A,T68,A)') "Density mixer", trim(active_space_mixing_label(active_space_env))
3253 WRITE (iw, '(T3,A,T66,F14.2)') "Mixing alpha", active_space_env%alpha
3254
3255 WRITE (unit=iw, fmt="(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") &
3256 "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", repeat("-", 78)
3257 END IF
3258 ! CALL cp_add_iter_level(logger%iter_info, "QS_SCF")
3259
3260 iter = 0
3261 converged = .false.
3262 ! store the scf energy
3263 energy_scf = active_space_env%energy_ref
3264 energy_new = energy_scf
3265 mos_active => active_space_env%mos_active
3266 ! CALL set_qs_env(qs_env, active_space=active_space_env)
3267
3268 ! start the self-consistent embedding loop
3269 DO WHILE (iter < max_iter)
3270 iter = iter + 1
3271
3272 ! send V_emb and E_ina to the active space solver and update
3273 ! the active space environment with the new active energy and density
3274 CALL send_fock_to_client(client_fd, active_space_env, para_env)
3275
3276 ! update energies
3277 energy_old = energy_new
3278 energy_new = active_space_env%energy_total
3279 energy_corr = energy_new - energy_scf
3280 delta_e = energy_new - energy_old
3281
3282 ! get timer
3283 t2 = t1
3284 t1 = m_walltime()
3285 ! print out progress
3286 IF ((iw > 0)) THEN
3287 WRITE (unit=iw, &
3288 fmt="(T3,I4,T11,A,T19,F6.1,T28,F18.10,T49,F18.10,T70,ES11.2)") &
3289 iter, trim(active_space_mixing_label(active_space_env)), &
3290 t1 - t2, energy_corr, energy_new, delta_e
3291 CALL m_flush(iw)
3292 END IF
3293
3294 ! update total density in AO basis with the AS contribution
3295 CALL update_density_ao(active_space_env, rho_ao) ! rho_ao is updated
3296
3297 ! calculate F_ks in AO basis (which contains Vxc) with the new density
3298 CALL qs_rho_update_rho(rho, qs_env=qs_env) ! updates rho_r and rho_g using rho_ao
3299 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.) ! set flags about the change
3300 ! Re-evaluate the traces between the density matrix and the core Hamiltonians
3301 CALL evaluate_core_matrix_traces(qs_env)
3302 ! the ks matrix will be rebuilt so this is fine now
3303 ! CALL set_ks_env(qs_env%ks_env, potential_changed=.FALSE.)
3304 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., &
3305 just_energy=.false., &
3306 ext_xc_section=active_space_env%xc_section)
3307
3308 ! update the reference energy
3309 active_space_env%energy_ref = energy%total
3310
3311 ! transform KS/Fock, Vxc and Hcore from AO to MO basis
3312 CALL calculate_operators(mos_active, qs_env, active_space_env)
3313
3314 ! calculate the new inactive energy and embedding potential
3315 CALL subspace_fock_matrix(active_space_env, dft_control%roks)
3316
3317 ! check if it is a one-shot correction
3318 IF (.NOT. active_space_env%do_scf_embedding) THEN
3319 IF (iw > 0) THEN
3320 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3321 "*** one-shot embedding correction finished ***"
3322 END IF
3323 converged = .true.
3324 EXIT
3325 ! check for convergence
3326 ELSE IF (abs(delta_e) <= eps_iter) THEN
3327 IF (iw > 0) THEN
3328 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3329 "*** rs-DFT embedding run converged in ", iter, " iteration(s) ***"
3330 END IF
3331 converged = .true.
3332 EXIT
3333 END IF
3334 END DO
3335
3336 IF (.NOT. converged) THEN
3337 IF (iw > 0) THEN
3338 WRITE (unit=iw, fmt="(/,T3,A,I5,A)") &
3339 "*** rs-DFT embedding did not converged after ", iter, " iteration(s) ***"
3340 END IF
3341 END IF
3342
3343 ! update qs total energy to the final rs-DFT energy
3344 energy%total = active_space_env%energy_total
3345
3346 ! print final energy contributions
3347 IF (iw > 0) THEN
3348 WRITE (unit=iw, fmt="(/,T3,A)") &
3349 "Final energy contributions:"
3350 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3351 "Inactive energy:", active_space_env%energy_inactive
3352 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3353 "Active energy:", active_space_env%energy_active
3354 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3355 "Correlation energy:", energy_corr
3356 WRITE (unit=iw, fmt="(T6,A,T56,F20.10)") &
3357 "Total rs-DFT energy:", active_space_env%energy_total
3358 END IF
3359
3360 ! print the AS rdm and the natural orbital occupation numbers
3361 CALL print_pmat_noon(active_space_env, iw)
3362
3363 CALL finalize_socket(socket_fd, client_fd, as_input, ionode)
3364 CALL para_env%sync()
3365#endif
3366
3367 CALL timestop(handle)
3368
3369 END SUBROUTINE rsdft_embedding
3370
3371#ifndef __NO_SOCKETS
3372! **************************************************************************************************
3373!> \brief Creates the socket, spawns the client and connects to it
3374!> \param socket_fd the socket file descriptor
3375!> \param client_fd the client file descriptor
3376!> \param as_input active space inpute section
3377!> \param ionode logical flag indicating if the process is the master
3378! **************************************************************************************************
3379 SUBROUTINE initialize_socket(socket_fd, client_fd, as_input, ionode)
3380 INTEGER, INTENT(OUT) :: socket_fd, client_fd
3381 TYPE(section_vals_type), INTENT(IN), POINTER :: as_input
3382 LOGICAL, INTENT(IN) :: ionode
3383
3384 CHARACTER(len=*), PARAMETER :: routinen = 'initialize_socket'
3385 INTEGER, PARAMETER :: backlog = 10
3386
3387 CHARACTER(len=default_path_length) :: hostname
3388 INTEGER :: handle, iw, port, protocol
3389 LOGICAL :: inet
3390 TYPE(cp_logger_type), POINTER :: logger
3391
3392 CALL timeset(routinen, handle)
3393
3394 logger => cp_get_default_logger()
3395 iw = cp_logger_get_default_io_unit(logger)
3396
3397 ! protocol == 0 for UNIX, protocol > 0 for INET
3398 CALL section_vals_val_get(as_input, "SOCKET%INET", l_val=inet)
3399 IF (inet) THEN
3400 protocol = 1
3401 ELSE
3402 protocol = 0
3403 END IF
3404 CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
3405 CALL section_vals_val_get(as_input, "SOCKET%PORT", i_val=port)
3406
3407 IF (ionode) THEN
3408 CALL open_bind_socket(socket_fd, protocol, port, trim(hostname)//c_null_char)
3409 WRITE (iw, '(/,T2,A,A)') "@SERVER: Created socket with address ", trim(hostname)
3410 CALL listen_socket(socket_fd, backlog)
3411
3412 ! wait until a connetion request arrives
3413 WRITE (iw, '(T2,A)') "@SERVER: Waiting for requests..."
3414 CALL accept_socket(socket_fd, client_fd)
3415 WRITE (iw, '(T2,A,I2)') "@SERVER: Accepted socket with fd ", client_fd
3416 END IF
3417
3418 CALL timestop(handle)
3419
3420 END SUBROUTINE initialize_socket
3421
3422! **************************************************************************************************
3423!> \brief Closes the connection to the socket and deletes the file
3424!> \param socket_fd the socket file descriptor
3425!> \param client_fd the client file descriptor
3426!> \param as_input active space inpute section
3427!> \param ionode logical flag indicating if the process is the master
3428! **************************************************************************************************
3429 SUBROUTINE finalize_socket(socket_fd, client_fd, as_input, ionode)
3430 INTEGER, INTENT(IN) :: socket_fd, client_fd
3431 TYPE(section_vals_type), INTENT(IN), POINTER :: as_input
3432 LOGICAL, INTENT(IN) :: ionode
3433
3434 CHARACTER(len=*), PARAMETER :: routinen = 'finalize_socket'
3435 INTEGER, PARAMETER :: header_len = 12
3436
3437 CHARACTER(len=default_path_length) :: hostname
3438 INTEGER :: handle
3439
3440 CALL timeset(routinen, handle)
3441
3442 CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
3443
3444 IF (ionode) THEN
3445 ! signal the client to quit
3446 CALL writebuffer(client_fd, "QUIT ", header_len)
3447 ! close the connection
3448 CALL close_socket(client_fd)
3449 CALL close_socket(socket_fd)
3450
3451 ! delete the socket file
3452 IF (file_exists(trim(hostname))) THEN
3453 CALL remove_socket_file(trim(hostname)//c_null_char)
3454 END IF
3455 END IF
3456
3457 CALL timestop(handle)
3458
3459 END SUBROUTINE finalize_socket
3460
3461! **************************************************************************************************
3462!> \brief Sends the two-electron integrals to the client vie the socket
3463!> \param client_fd the client file descriptor
3464!> \param active_space_env active space environment
3465!> \param para_env parallel environment
3466! **************************************************************************************************
3467 SUBROUTINE send_eri_to_client(client_fd, active_space_env, para_env)
3468 INTEGER, INTENT(IN) :: client_fd
3469 TYPE(active_space_type), INTENT(IN), POINTER :: active_space_env
3470 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
3471
3472 CHARACTER(len=*), PARAMETER :: routinen = 'send_eri_to_client'
3473 INTEGER, PARAMETER :: header_len = 12
3474
3475 CHARACTER(len=default_string_length) :: header
3476 INTEGER :: handle, iw
3477 LOGICAL :: ionode, restricted_orbitals
3478 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb, s_ab
3479 TYPE(cp_logger_type), POINTER :: logger
3480
3481 CALL timeset(routinen, handle)
3482
3483 logger => cp_get_default_logger()
3484 iw = cp_logger_get_default_io_unit(logger)
3485 ionode = para_env%is_source()
3486 restricted_orbitals = active_space_env%restricted_orbitals
3487
3488 ALLOCATE (eri_aa(active_space_env%nmo_active**4))
3489 CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
3490 IF (active_space_env%nspins == 2) THEN
3491 ALLOCATE (eri_ab(active_space_env%nmo_active**4))
3492 ALLOCATE (eri_bb(active_space_env%nmo_active**4))
3493 IF (restricted_orbitals) THEN
3494 eri_ab(:) = eri_aa
3495 eri_bb(:) = eri_aa
3496 ELSE
3497 CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
3498 CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2)
3499 END IF
3500 ! get the overlap_ab matrix into Fortran array
3501 ALLOCATE (s_ab(active_space_env%nmo_active**2))
3502 associate(act_indices_a => active_space_env%active_orbitals(:, 1), &
3503 act_indices_b => active_space_env%active_orbitals(:, 2))
3504 CALL subspace_matrix_to_array(active_space_env%sab_sub(1), s_ab, act_indices_a, act_indices_b)
3505 END associate
3506 END IF
3507
3508 ! ask the status of the client
3509 IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len)
3510 DO
3511 header = ""
3512 CALL para_env%sync()
3513 IF (ionode) THEN
3514 ! IF (iw > 0) WRITE(iw, *) "@SERVER: Waiting for messages..."
3515 CALL readbuffer(client_fd, header, header_len)
3516 END IF
3517 CALL para_env%bcast(header, para_env%source)
3518
3519 ! IF (iw > 0) WRITE(iw, *) "@SERVER: Message from client: ", TRIM(header)
3520
3521 IF (trim(header) == "READY") THEN
3522 ! if the client is ready, send the data
3523 CALL para_env%sync()
3524 IF (ionode) THEN
3525 CALL writebuffer(client_fd, "TWOBODY ", header_len)
3526 CALL writebuffer(client_fd, active_space_env%nspins)
3527 CALL writebuffer(client_fd, active_space_env%nmo_active)
3528 CALL writebuffer(client_fd, active_space_env%nelec_active)
3529 CALL writebuffer(client_fd, active_space_env%multiplicity)
3530 ! send the alpha component
3531 CALL writebuffer(client_fd, eri_aa, SIZE(eri_aa))
3532 ! send the beta part for unrestricted calculations
3533 IF (active_space_env%nspins == 2) THEN
3534 CALL writebuffer(client_fd, eri_ab, SIZE(eri_ab))
3535 CALL writebuffer(client_fd, eri_bb, SIZE(eri_bb))
3536 CALL writebuffer(client_fd, s_ab, SIZE(s_ab))
3537 END IF
3538 END IF
3539 ELSE IF (trim(header) == "RECEIVED") THEN
3540 EXIT
3541 END IF
3542 END DO
3543
3544 DEALLOCATE (eri_aa)
3545 IF (active_space_env%nspins == 2) THEN
3546 DEALLOCATE (eri_ab)
3547 DEALLOCATE (eri_bb)
3548 DEALLOCATE (s_ab)
3549 END IF
3550
3551 CALL para_env%sync()
3552
3553 CALL timestop(handle)
3554
3555 END SUBROUTINE send_eri_to_client
3556
3557! **************************************************************************************************
3558!> \brief Sends the one-electron embedding potential and the inactive energy to the client
3559!> \param client_fd the client file descriptor
3560!> \param active_space_env active space environment
3561!> \param para_env parallel environment
3562! **************************************************************************************************
3563 SUBROUTINE send_fock_to_client(client_fd, active_space_env, para_env)
3564 INTEGER, INTENT(IN) :: client_fd
3565 TYPE(active_space_type), INTENT(INOUT), POINTER :: active_space_env
3566 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
3567
3568 CHARACTER(len=*), PARAMETER :: routinen = 'send_fock_to_client'
3569 INTEGER, PARAMETER :: header_len = 12
3570
3571 CHARACTER(len=default_string_length) :: header
3572 INTEGER :: handle, iw
3573 LOGICAL :: debug, ionode
3574 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fock_a, fock_b, p_act_mo_a, p_act_mo_b
3575 TYPE(cp_logger_type), POINTER :: logger
3576
3577 CALL timeset(routinen, handle)
3578
3579 ! Set to .TRUE. to activate debug output
3580 debug = .false.
3581
3582 logger => cp_get_default_logger()
3583 iw = cp_logger_get_default_io_unit(logger)
3584 ionode = para_env%is_source()
3585
3586 ALLOCATE (p_act_mo_a(active_space_env%nmo_active**2))
3587 ALLOCATE (fock_a(active_space_env%nmo_active**2))
3588 IF (active_space_env%nspins == 2) THEN
3589 ALLOCATE (p_act_mo_b(active_space_env%nmo_active**2))
3590 ALLOCATE (fock_b(active_space_env%nmo_active**2))
3591 END IF
3592
3593 ! get the fock matrix into Fortran arrays
3594 associate(act_indices => active_space_env%active_orbitals(:, 1))
3595 CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices)
3596 END associate
3597
3598 IF (active_space_env%nspins == 2) THEN
3599 associate(act_indices => active_space_env%active_orbitals(:, 2))
3600 CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b, act_indices, act_indices)
3601 END associate
3602 END IF
3603
3604 ! ask the status of the client
3605 IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len)
3606 DO
3607 header = ""
3608
3609 CALL para_env%sync()
3610 IF (ionode) THEN
3611 IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Waiting for messages..."
3612 CALL readbuffer(client_fd, header, header_len)
3613 END IF
3614 CALL para_env%bcast(header, para_env%source)
3615
3616 IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Message from client: ", trim(header)
3617
3618 IF (trim(header) == "READY") THEN
3619 ! if the client is ready, send the data
3620 CALL para_env%sync()
3621 IF (ionode) THEN
3622 CALL writebuffer(client_fd, "ONEBODY ", header_len)
3623 CALL writebuffer(client_fd, active_space_env%energy_inactive)
3624 ! send the alpha component
3625 CALL writebuffer(client_fd, fock_a, SIZE(fock_a))
3626 ! send the beta part for unrestricted calculations
3627 IF (active_space_env%nspins == 2) THEN
3628 CALL writebuffer(client_fd, fock_b, SIZE(fock_b))
3629 END IF
3630 END IF
3631
3632 ELSE IF (trim(header) == "HAVEDATA") THEN
3633 ! qiskit has data to transfer, let them know we want it and wait for it
3634 CALL para_env%sync()
3635 IF (ionode) THEN
3636 IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Qiskit has data to transfer"
3637 CALL writebuffer(client_fd, "GETDENSITY ", header_len)
3638
3639 ! read the active energy and density
3640 CALL readbuffer(client_fd, active_space_env%energy_active)
3641 CALL readbuffer(client_fd, p_act_mo_a, SIZE(p_act_mo_a))
3642 IF (active_space_env%nspins == 2) THEN
3643 CALL readbuffer(client_fd, p_act_mo_b, SIZE(p_act_mo_b))
3644 END IF
3645 END IF
3646
3647 ! broadcast the data to all processors
3648 CALL para_env%bcast(active_space_env%energy_active, para_env%source)
3649 CALL para_env%bcast(p_act_mo_a, para_env%source)
3650 IF (active_space_env%nspins == 2) THEN
3651 CALL para_env%bcast(p_act_mo_b, para_env%source)
3652 END IF
3653
3654 ! update total and reference energies in active space enviornment
3655 active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3656
3657 ! update the active density matrix in the active space environment
3658 IF (active_space_env%nspins == 2) THEN
3659 CALL update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
3660 ELSE
3661 CALL update_active_density(p_act_mo_a, active_space_env)
3662 END IF
3663
3664 ! the non-iterative part is done, we can continue
3665 EXIT
3666 END IF
3667
3668 END DO
3669
3670 DEALLOCATE (p_act_mo_a)
3671 DEALLOCATE (fock_a)
3672 IF (active_space_env%nspins == 2) THEN
3673 DEALLOCATE (p_act_mo_b)
3674 DEALLOCATE (fock_b)
3675 END IF
3676
3677 CALL para_env%sync()
3678
3679 CALL timestop(handle)
3680
3681 END SUBROUTINE send_fock_to_client
3682#endif
3683
3684END MODULE qs_active_space_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
group_dist_types::get_group_dist
Definition group_dist_types.F:47
group_dist_types::release_group_dist
Definition group_dist_types.F:53
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_multiply_with
Definition pw_methods.F:203
pw_methods::pw_multiply
Definition pw_methods.F:184
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_rebuild
Definition pw_poisson_methods.F:73
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
qs_density_matrices::calculate_density_matrix
Definition qs_density_matrices.F:58
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:68
sockets_interface::accept_socket
Definition sockets_interface.F:79
sockets_interface::close_socket
Definition sockets_interface.F:85
sockets_interface::listen_socket
Definition sockets_interface.F:73
sockets_interface::open_bind_socket
Definition sockets_interface.F:66
sockets_interface::readbuffer
Definition sockets_interface.F:43
sockets_interface::remove_socket_file
Definition sockets_interface.F:91
sockets_interface::writebuffer
Definition sockets_interface.F:36
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:599
admm_types::admm_env_release
subroutine, public admm_env_release(admm_env)
releases the ADMM environment, cleans up all types
Definition admm_types.F:429
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
basis_set_types
Definition basis_set_types.F:15
basis_set_types::srules
pure real(dp) function, public srules(z, ne, n, l)
...
Definition basis_set_types.F:2707
basis_set_types::deallocate_sto_basis_set
subroutine, public deallocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:2136
basis_set_types::allocate_sto_basis_set
subroutine, public allocate_sto_basis_set(sto_basis_set)
...
Definition basis_set_types.F:2122
basis_set_types::create_gto_from_sto_basis
subroutine, public create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
...
Definition basis_set_types.F:2426
basis_set_types::set_sto_basis_set
subroutine, public set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
...
Definition basis_set_types.F:2213
basis_set_types::init_orb_basis_set
subroutine, public init_orb_basis_set(gto_basis_set)
Initialise a Gaussian-type orbital (GTO) basis set data set.
Definition basis_set_types.F:1138
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::write_cell_low
subroutine, public write_cell_low(cell, unit_str, output_unit, label)
Write the cell parameters to the output unit.
Definition cell_methods.F:1396
cell_methods::set_cell_param
subroutine, public set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma) using the convention: a parall...
Definition cell_methods.F:672
cell_methods::init_cell
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Definition cell_methods.F:185
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::use_perd_none
integer, parameter, public use_perd_none
Definition cell_types.F:42
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::blacs_grid_square
integer, parameter, public blacs_grid_square
Definition cp_blacs_env.F:32
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:893
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:693
cp_dbcsr_operations::cp_dbcsr_m_by_n_from_template
subroutine, public cp_dbcsr_m_by_n_from_template(matrix, template, m, n, sym)
Utility function to create an arbitrary shaped dbcsr matrix with the same processor grid as the templ...
Definition cp_dbcsr_operations.F:892
cp_dbcsr_output
DBCSR output in CP2K.
Definition cp_dbcsr_output.F:17
cp_dbcsr_output::cp_dbcsr_write_sparse_matrix
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers, cartesian_basis)
...
Definition cp_dbcsr_output.F:165
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:504
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:627
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_set_element
subroutine, public cp_fm_set_element(matrix, irow_global, icol_global, alpha)
sets an element of a matrix
Definition cp_fm_types.F:684
cp_fm_types::cp_fm_init_random
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
Definition cp_fm_types.F:452
cp_fm_types::cp_fm_write_formatted
subroutine, public cp_fm_write_formatted(fm, unit, header, value_format)
Write out a full matrix in plain text.
Definition cp_fm_types.F:2236
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::debug_print_level
integer, parameter, public debug_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::high_print_level
integer, parameter, public high_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_output_handling::silent_print_level
integer, parameter, public silent_print_level
Definition cp_output_handling.F:73
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
erf_complex
Module to compute the error function of a complex argument.
Definition erf_complex.F:50
erf_complex::erfz_fast
elemental complex(kind=dp) function, public erfz_fast(z)
Computes the error function of a complex argument using the Poppe and Wijers algorithm.
Definition erf_complex.F:475
group_dist_types
Types to describe group distributions.
Definition group_dist_types.F:14
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
hfx_types::hfx_create
subroutine, public hfx_create(x_data, para_env, hfx_section, atomic_kind_set, qs_kind_set, particle_set, dft_control, cell, orb_basis, ri_basis, nelectron_total, nkp_grid)
This routine allocates and initializes all types in hfx_data
Definition hfx_types.F:607
hfx_types::hfx_release
subroutine, public hfx_release(x_data)
This routine deallocates all data structures
Definition hfx_types.F:1971
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::eri_operator_erf
integer, parameter, public eri_operator_erf
Definition input_constants.F:1290
input_constants::qiskit_solver
integer, parameter, public qiskit_solver
Definition input_constants.F:1181
input_constants::no_solver
integer, parameter, public no_solver
Definition input_constants.F:1181
input_constants::wannier_projection
integer, parameter, public wannier_projection
Definition input_constants.F:1284
input_constants::mao_projection
integer, parameter, public mao_projection
Definition input_constants.F:1284
input_constants::eri_method_full_gpw
integer, parameter, public eri_method_full_gpw
Definition input_constants.F:1288
input_constants::casci_canonical
integer, parameter, public casci_canonical
Definition input_constants.F:1284
input_constants::manual_selection
integer, parameter, public manual_selection
Definition input_constants.F:1284
input_constants::eri_operator_gaussian
integer, parameter, public eri_operator_gaussian
Definition input_constants.F:1290
input_constants::high_spin_roks
integer, parameter, public high_spin_roks
Definition input_constants.F:813
input_constants::eri_method_gpw_ht
integer, parameter, public eri_method_gpw_ht
Definition input_constants.F:1288
input_constants::eri_operator_erfc
integer, parameter, public eri_operator_erfc
Definition input_constants.F:1290
input_constants::eri_poisson_mt
integer, parameter, public eri_poisson_mt
Definition input_constants.F:1298
input_constants::eri_poisson_analytic
integer, parameter, public eri_poisson_analytic
Definition input_constants.F:1298
input_constants::eri_operator_trunc
integer, parameter, public eri_operator_trunc
Definition input_constants.F:1290
input_constants::eri_poisson_periodic
integer, parameter, public eri_poisson_periodic
Definition input_constants.F:1298
input_constants::eri_operator_coulomb
integer, parameter, public eri_operator_coulomb
Definition input_constants.F:1290
input_constants::eri_operator_yukawa
integer, parameter, public eri_operator_yukawa
Definition input_constants.F:1290
input_constants::eri_operator_lr_trunc
integer, parameter, public eri_operator_lr_trunc
Definition input_constants.F:1290
input_constants::fci_solver
integer, parameter, public fci_solver
Definition input_constants.F:1181
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1233
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_set_subs_vals
subroutine, public section_vals_set_subs_vals(section_vals, subsection_name, new_section_vals, i_rep_section)
replaces of the requested subsection with the one given
Definition input_section_types.F:1723
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:529
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_para_env_release
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Definition message_passing.F:2029
mp2_gpw
Calls routines to get RI integrals and calculate total energies.
Definition mp2_gpw.F:14
mp2_gpw::grep_rows_in_subgroups
subroutine, public grep_rows_in_subgroups(para_env, para_env_sub, mo_coeff, gd_array, c)
...
Definition mp2_gpw.F:758
mp2_gpw::build_dbcsr_from_rows
subroutine, public build_dbcsr_from_rows(para_env_sub, mo_coeff_to_build, cread, mat_munu, gd_array, eps_filter)
Encapsulate the building of dbcsr_matrices mo_coeff_(v,o,all)
Definition mp2_gpw.F:886
mp2_gpw::create_mat_munu
subroutine, public create_mat_munu(mat_munu, qs_env, eps_grid, blacs_env_sub, do_ri_aux_basis, do_mixed_basis, group_size_prim, do_alloc_blocks_from_nbl, do_kpoints, sab_orb_sub, dbcsr_sym_type)
Encapsulate the building of dbcsr_matrix mat_munu.
Definition mp2_gpw.F:989
mt_util
Definition mt_util.F:9
mt_util::mt0d
integer, parameter, public mt0d
Definition mt_util.F:33
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
preconditioner_types
types of preconditioners
Definition preconditioner_types.F:14
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition preconditioner.F:15
pw_env_methods
methods of pw_env that have dependence on qs_env
Definition pw_env_methods.F:17
pw_env_methods::pw_env_rebuild
subroutine, public pw_env_rebuild(pw_env, qs_env, external_para_env)
rebuilds the pw_env data (necessary if cell or cutoffs change)
Definition pw_env_methods.F:164
pw_env_methods::pw_env_create
subroutine, public pw_env_create(pw_env)
creates a pw_env, if qs_env is given calls pw_env_rebuild
Definition pw_env_methods.F:133
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_release
subroutine, public pw_env_release(pw_env, para_env)
releases the given pw_env (see doc/ReferenceCounting.html)
Definition pw_env_types.F:176
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_methods::pw_gauss_damp
subroutine, public pw_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
Definition pw_methods.F:9943
pw_methods::pw_compl_gauss_damp
subroutine, public pw_compl_gauss_damp(pw, omega)
Multiply all data points with a Gaussian damping factor Needed for longrange Coulomb potential V(\vec...
Definition pw_methods.F:10016
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::pw_poisson_periodic
integer, parameter, public pw_poisson_periodic
Definition pw_poisson_types.F:67
pw_poisson_types::periodic3d
integer, parameter, public periodic3d
Definition pw_poisson_types.F:81
pw_poisson_types::analytic0d
integer, parameter, public analytic0d
Definition pw_poisson_types.F:81
pw_poisson_types::pw_poisson_analytic
integer, parameter, public pw_poisson_analytic
Definition pw_poisson_types.F:67
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qcschema
The module to read/write QCSchema HDF5 files for interfacing CP2K with other programs.
Definition qcschema.F:14
qcschema::qcschema_env_release
subroutine, public qcschema_env_release(qcschema_env)
Releases the allocated memory of a qcschema environment.
Definition qcschema.F:611
qcschema::qcschema_env_create
subroutine, public qcschema_env_create(qcschema_env, qs_env)
Create and initialize a qcschema object from a quickstep environment.
Definition qcschema.F:272
qcschema::qcschema_to_hdf5
subroutine, public qcschema_to_hdf5(qcschema_env, filename)
Writes a qcschema object to an hdf5 file.
Definition qcschema.F:741
qs_active_space_fci
Local FCI solver interface for active-space embedding.
Definition qs_active_space_fci.F:11
qs_active_space_fci::solve_active_space_fci
subroutine, public solve_active_space_fci(active_space_env, para_env, p_act_mo_a, p_act_mo_b)
Solve the current active-space Hamiltonian with the local FCI kernel.
Definition qs_active_space_fci.F:64
qs_active_space_methods
Determine active space Hamiltonian.
Definition qs_active_space_methods.F:14
qs_active_space_methods::eri_fcidump_set
subroutine eri_fcidump_set(this, bra_start, ket_start)
Sets the starting indices of the bra and ket.
Definition qs_active_space_methods.F:206
qs_active_space_methods::active_space_main
subroutine, public active_space_main(qs_env)
Main method for determining the active space Hamiltonian.
Definition qs_active_space_methods.F:217
qs_active_space_mixing
Dense density mixing for active-space embedding.
Definition qs_active_space_mixing.F:11
qs_active_space_mixing::active_space_mixing_label
character(len=15) function, public active_space_mixing_label(active_space_env)
Return the current active-space mixer label for iteration output.
Definition qs_active_space_mixing.F:123
qs_active_space_mixing::update_active_density
subroutine, public update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
Update active space density matrix from Fortran arrays.
Definition qs_active_space_mixing.F:502
qs_active_space_mixing::initialize_active_space_mixing
subroutine, public initialize_active_space_mixing(active_space_env, as_input)
Initialize the density mixer used in the self-consistent active-space embedding loop.
Definition qs_active_space_mixing.F:43
qs_active_space_types
The types needed for the calculation of active space Hamiltonians.
Definition qs_active_space_types.F:14
qs_active_space_types::csr_idx_from_combined
subroutine, public csr_idx_from_combined(ij, n, i, j)
extracts indices i and j from combined index ij
Definition qs_active_space_types.F:304
qs_active_space_types::csr_idx_to_combined
integer function, public csr_idx_to_combined(i, j, n)
calculates combined index (ij)
Definition qs_active_space_types.F:282
qs_active_space_types::create_active_space_type
subroutine, public create_active_space_type(active_space_env)
Creates an active space environment type, nullifying all quantities.
Definition qs_active_space_types.F:159
qs_active_space_utils
Contains utility routines for the active space module.
Definition qs_active_space_utils.F:14
qs_active_space_utils::eri_to_array
subroutine, public eri_to_array(eri_env, array, active_orbitals, spin1, spin2)
Copy the eri tensor for spins isp1 and isp2 to a standard 1D Fortran array.
Definition qs_active_space_utils.F:85
qs_active_space_utils::subspace_matrix_to_array
subroutine, public subspace_matrix_to_array(source_matrix, target_array, row_index, col_index)
Copy a (square portion) of a cp_fm_type matrix to a standard 1D Fortran array.
Definition qs_active_space_utils.F:45
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_wavefunction
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type)
maps a given wavefunction on the grid
Definition qs_collocate_density.F:2869
qs_density_matrices
collects routines that calculate density matrices
Definition qs_density_matrices.F:14
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1120
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::evaluate_core_matrix_traces
subroutine, public evaluate_core_matrix_traces(qs_env, rho_ao_ext)
Calculates the traces of the core matrices and the density matrix.
Definition qs_ks_methods.F:1611
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition qs_ks_methods.F:1540
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix, ext_xc_section)
routine where the real calculations are made: the KS matrix is calculated
Definition qs_ks_methods.F:201
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition qs_ks_types.F:977
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::write_mo_set_to_output_unit
subroutine, public write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set, qs_env)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
Definition qs_mo_io.F:1023
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name, counter)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:247
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:207
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1096
qs_ot_eigensolver
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
Definition qs_ot_eigensolver.F:13
qs_ot_eigensolver::ot_eigensolver
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Definition qs_ot_eigensolver.F:74
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition qs_rho_methods.F:380
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_post_scf
Utility routines for qs_scf.
Definition qs_scf_post_scf.F:11
qs_scf_post_scf::qs_scf_compute_properties
subroutine, public qs_scf_compute_properties(qs_env, wf_type, do_mp2)
computes properties for a given hamilonian using the current wfn
Definition qs_scf_post_scf.F:45
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:137
scf_control_types
parameters that control an scf iteration
Definition scf_control_types.F:17
sockets_interface
Implements UNIX and INET sockets.
Definition sockets_interface.F:16
task_list_methods
generate the tasks lists used by collocate and integrate routines
Definition task_list_methods.F:14
task_list_methods::generate_qs_task_list
subroutine, public generate_qs_task_list(ks_env, task_list, basis_type, reorder_rs_grid_ranks, skip_load_balance_distributed, pw_env_external, sab_orb_external, ext_kpoints)
...
Definition task_list_methods.F:119
task_list_types
types for task lists
Definition task_list_types.F:14
task_list_types::deallocate_task_list
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
Definition task_list_types.F:145
task_list_types::allocate_task_list
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
Definition task_list_types.F:97
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
open_bind_socket
void open_bind_socket(int *psockfd, int *inet, int *port, char *host)
Opens and binds a socket.
Definition sockets.c:139
writebuffer
void writebuffer(int *psockfd, char *data, int *plen)
Writes to a socket.
Definition sockets.c:203
accept_socket
void accept_socket(int *psockfd, int *pclientfd)
Listens to a socket.
Definition sockets.c:258
close_socket
void close_socket(int *psockfd)
Closes a socket.
Definition sockets.c:269
listen_socket
void listen_socket(int *psockfd, int *backlog)
Listens to a socket.
Definition sockets.c:245
readbuffer
void readbuffer(int *psockfd, char *data, int *plen)
Reads from a socket.
Definition sockets.c:221
remove_socket_file
void remove_socket_file(char *host)
Removes a socket file.
Definition sockets.c:275
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
basis_set_types::sto_basis_set_type
Definition basis_set_types.F:101
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_control_types::qs_control_type
Definition cp_control_types.F:499
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_comm_type
Definition message_passing.F:161
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
preconditioner_types::preconditioner_type
Definition preconditioner_types.F:42
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::greens_fn_type
contains all the informations needed by the fft based poisson solvers
Definition pw_poisson_types.F:146
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qcschema::qcschema_type
The full QCSchema output type. For more information refer to: https://molssi-qc-schema....
Definition qcschema.F:253
qs_active_space_types::active_space_type
Definition qs_active_space_types.F:87
qs_active_space_types::eri_type_eri_element_func
Abstract function object for the eri_type_eri_foreach method.
Definition qs_active_space_types.F:82
qs_active_space_types::eri_type
Definition qs_active_space_types.F:59
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56
scf_control_types::scf_control_type
Definition scf_control_types.F:126
task_list_types::task_list_type
Definition task_list_types.F:58