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