Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Determine active space Hamiltonian
10 : !> \par History
11 : !> 04.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE 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 :
206 : CONTAINS
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 128 : SUBROUTINE eri_fcidump_set(this, bra_start, ket_start)
215 : CLASS(eri_fcidump_checksum) :: this
216 : INTEGER, INTENT(IN) :: bra_start, ket_start
217 128 : this%bra_start = bra_start
218 128 : this%ket_start = ket_start
219 128 : END SUBROUTINE eri_fcidump_set
220 :
221 : ! **************************************************************************************************
222 : !> \brief Main method for determining the active space Hamiltonian
223 : !> \param qs_env ...
224 : ! **************************************************************************************************
225 29031 : 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 29031 : 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 29031 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
246 29031 : REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virtual
247 : TYPE(active_space_type), POINTER :: active_space_env
248 29031 : 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 58062 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, rho_ao, s_matrix
258 58062 : 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 29031 : 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 29031 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
266 : TYPE(preconditioner_type), POINTER :: local_preconditioner
267 116124 : TYPE(qcschema_type) :: qcschema_env
268 : TYPE(qs_energy_type), POINTER :: energy
269 29031 : 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 29031 : CALL get_qs_env(qs_env, input=input)
280 29031 : as_input => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE")
281 29031 : CALL section_vals_get(as_input, explicit=explicit)
282 29031 : IF (.NOT. explicit) RETURN
283 108 : CALL timeset(routineN, handle)
284 :
285 108 : logger => cp_get_default_logger()
286 108 : iw = cp_logger_get_default_io_unit(logger)
287 :
288 108 : IF (iw > 0) THEN
289 : WRITE (iw, '(/,T2,A)') &
290 54 : '!-----------------------------------------------------------------------------!'
291 54 : WRITE (iw, '(T26,A)') "Active Space Embedding Module"
292 : WRITE (iw, '(T2,A)') &
293 54 : '!-----------------------------------------------------------------------------!'
294 : END IF
295 :
296 : ! k-points?
297 108 : NULLIFY (kpoints)
298 108 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, dft_control=dft_control, kpoints=kpoints)
299 108 : IF (do_kpoints) THEN
300 2 : IF (.NOT. ASSOCIATED(kpoints)) THEN
301 0 : CALL cp_abort(__LOCATION__, "Missing Gamma-point environment for active space module")
302 : END IF
303 2 : CALL get_kpoint_info(kpoints, kp_scheme=kp_scheme, nkp=nkp, use_real_wfn=use_real_wfn)
304 2 : IF (TRIM(kp_scheme) /= "GAMMA" .OR. nkp /= 1 .OR. .NOT. use_real_wfn) THEN
305 : CALL cp_abort(__LOCATION__, &
306 0 : "Only Gamma-point DFT%KPOINTS are supported in the active space module")
307 : END IF
308 2 : IF (.NOT. ASSOCIATED(kpoints%kp_env)) THEN
309 0 : CALL cp_abort(__LOCATION__, "Missing Gamma-point environment for active space module")
310 : END IF
311 2 : IF (.NOT. ASSOCIATED(kpoints%kp_env(1)%kpoint_env)) THEN
312 0 : CALL cp_abort(__LOCATION__, "Missing Gamma-point environment for active space module")
313 : END IF
314 2 : IF (.NOT. ASSOCIATED(kpoints%kp_env(1)%kpoint_env%mos)) THEN
315 0 : CALL cp_abort(__LOCATION__, "Missing Gamma-point MOs for active space module")
316 : END IF
317 : END IF
318 :
319 : ! adiabatic rescaling?
320 108 : adiabatic_rescaling => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
321 108 : CALL section_vals_get(adiabatic_rescaling, explicit=explicit)
322 108 : IF (explicit) THEN
323 0 : 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 108 : qs_env%cp_ddapc_ewald%do_solvation
331 : IF (do_ddapc) THEN
332 0 : CALL cp_abort(__LOCATION__, "DDAPC charges are not supported in the active space module")
333 : END IF
334 108 : IF (dft_control%do_sccs) THEN
335 0 : CALL cp_abort(__LOCATION__, "SCCS is not supported in the active space module")
336 : END IF
337 108 : IF (dft_control%correct_surf_dip) THEN
338 0 : IF (dft_control%surf_dip_correct_switch) THEN
339 0 : CALL cp_abort(__LOCATION__, "Surface dipole correction not supported in the AS module")
340 : END IF
341 : END IF
342 108 : IF (dft_control%smeagol_control%smeagol_enabled) THEN
343 0 : CALL cp_abort(__LOCATION__, "SMEAGOL is not supported in the active space module")
344 : END IF
345 108 : IF (dft_control%qs_control%do_kg) THEN
346 0 : CALL cp_abort(__LOCATION__, "KG correction not supported in the active space module")
347 : END IF
348 :
349 108 : NULLIFY (active_space_env)
350 108 : CALL create_active_space_type(active_space_env)
351 108 : active_space_env%energy_total = 0.0_dp
352 108 : active_space_env%energy_ref = 0.0_dp
353 108 : active_space_env%energy_inactive = 0.0_dp
354 108 : active_space_env%energy_active = 0.0_dp
355 :
356 : ! input options
357 :
358 108 : CALL section_vals_val_get(as_input, "EPS_FCI", r_val=active_space_env%eps_fci)
359 108 : IF (.NOT. (active_space_env%eps_fci > 0.0_dp)) THEN
360 0 : CPABORT("ACTIVE_SPACE%EPS_FCI must be positive.")
361 : END IF
362 :
363 : ! figure out what needs to be printed/stored
364 108 : IF (BTEST(cp_print_key_should_output(logger%iter_info, as_input, "FCIDUMP"), cp_p_file)) THEN
365 104 : active_space_env%fcidump = .TRUE.
366 : END IF
367 :
368 108 : CALL section_vals_val_get(as_input, "QCSCHEMA", c_val=qcschema_filename, explicit=explicit)
369 108 : IF (explicit) THEN
370 6 : active_space_env%qcschema = .TRUE.
371 6 : active_space_env%qcschema_filename = qcschema_filename
372 : END IF
373 :
374 108 : CALL section_vals_val_get(as_input, "ACTIVE_ELECTRONS", i_val=nelec_active)
375 108 : CALL get_qs_env(qs_env, nelectron_total=nelec_total)
376 :
377 108 : IF (nelec_active <= 0) CPABORT("Specify a positive number of active electrons.")
378 108 : IF (nelec_active > nelec_total) CPABORT("More active electrons than total electrons.")
379 :
380 108 : nelec_inactive = nelec_total - nelec_active
381 108 : IF (MOD(nelec_inactive, 2) /= 0) THEN
382 0 : CPABORT("The remaining number of inactive electrons has to be even.")
383 : END IF
384 :
385 108 : IF (iw > 0) THEN
386 54 : WRITE (iw, '(T3,A,T70,I10)') "Total number of electrons", nelec_total
387 54 : WRITE (iw, '(T3,A,T70,I10)') "Number of inactive electrons", nelec_inactive
388 54 : WRITE (iw, '(T3,A,T70,I10)') "Number of active electrons", nelec_active
389 : END IF
390 :
391 108 : CALL get_qs_env(qs_env, dft_control=dft_control)
392 108 : nspins = dft_control%nspins
393 :
394 108 : active_space_env%nelec_active = nelec_active
395 108 : active_space_env%nelec_inactive = nelec_inactive
396 108 : active_space_env%nelec_total = nelec_total
397 108 : active_space_env%nspins = nspins
398 108 : active_space_env%multiplicity = dft_control%multiplicity
399 108 : active_space_env%restricted_orbitals = dft_control%roks
400 :
401 : ! define the active/inactive space orbitals
402 108 : CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit, i_val=nmo_active)
403 108 : IF (.NOT. explicit) THEN
404 0 : CALL cp_abort(__LOCATION__, "Number of Active Orbitals has to be specified.")
405 : END IF
406 108 : active_space_env%nmo_active = nmo_active
407 : ! this is safe because nelec_inactive is always even
408 108 : nmo_inactive = nelec_inactive/2
409 108 : active_space_env%nmo_inactive = nmo_inactive
410 :
411 108 : CALL initialize_active_space_mixing(active_space_env, as_input)
412 :
413 108 : CALL section_vals_val_get(as_input, "ORBITAL_SELECTION", i_val=mselect)
414 108 : IF (iw > 0) THEN
415 0 : SELECT CASE (mselect)
416 : CASE DEFAULT
417 0 : CPABORT("Unknown orbital selection method")
418 : CASE (casci_canonical)
419 : WRITE (iw, '(/,T3,A)') &
420 39 : "Active space orbitals selected using energy ordered canonical orbitals"
421 : CASE (wannier_projection)
422 : WRITE (iw, '(/,T3,A)') &
423 0 : "Active space orbitals selected using projected Wannier orbitals"
424 : CASE (mao_projection)
425 : WRITE (iw, '(/,T3,A)') &
426 0 : "Active space orbitals selected using modified atomic orbitals (MAO)"
427 : CASE (manual_selection)
428 : WRITE (iw, '(/,T3,A)') &
429 54 : "Active space orbitals selected manually"
430 : END SELECT
431 :
432 54 : WRITE (iw, '(T3,A,T70,I10)') "Number of inactive orbitals", nmo_inactive
433 54 : WRITE (iw, '(T3,A,T70,I10)') "Number of active orbitals", nmo_active
434 : END IF
435 :
436 : ! get projection spaces
437 108 : CALL section_vals_val_get(as_input, "SUBSPACE_ATOM", i_val=iatom, explicit=explicit)
438 108 : IF (explicit) THEN
439 0 : CALL get_qs_env(qs_env, natom=natom)
440 0 : IF (iatom <= 0 .OR. iatom > natom) THEN
441 0 : IF (iw > 0) THEN
442 0 : WRITE (iw, '(/,T3,A,I3)') "ERROR: SUBSPACE_ATOM number is not valid", iatom
443 : END IF
444 0 : CPABORT("Select a valid SUBSPACE_ATOM")
445 : END IF
446 : END IF
447 108 : CALL section_vals_val_get(as_input, "SUBSPACE_SHELL", c_val=cshell, explicit=explicit)
448 108 : nshell = 0
449 648 : lnam = ""
450 108 : IF (explicit) THEN
451 0 : cshell = ADJUSTL(cshell)
452 0 : n1 = 1
453 0 : DO i = 1, 5
454 0 : ishell = i
455 0 : IF (cshell(n1:n1) == " ") THEN
456 108 : ishell = ishell - 1
457 : EXIT
458 : END IF
459 0 : READ (cshell(n1:), "(I1,A1)") nshell(i), lnam(i)
460 0 : n1 = n1 + 2
461 : END DO
462 : END IF
463 :
464 : ! generate orbitals
465 0 : SELECT CASE (mselect)
466 : CASE DEFAULT
467 0 : CPABORT("Unknown orbital selection method")
468 : CASE (casci_canonical)
469 78 : IF (do_kpoints) THEN
470 2 : mos => kpoints%kp_env(1)%kpoint_env%mos(:)
471 : ELSE
472 76 : 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 78 : nmo_occ = nmo_inactive + nmo_active
477 :
478 : ! set inactive orbital indices, these are trivially 1...nmo_inactive
479 260 : ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
480 166 : DO ispin = 1, nspins
481 238 : DO i = 1, nmo_inactive
482 160 : 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 312 : ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
488 166 : DO ispin = 1, nspins
489 432 : DO i = 1, nmo_active
490 354 : 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 78 : maxocc = 2.0_dp
499 78 : IF (nspins > 1) maxocc = 1.0_dp
500 322 : ALLOCATE (active_space_env%mos_active(nspins))
501 244 : ALLOCATE (active_space_env%mos_inactive(nspins))
502 166 : DO ispin = 1, nspins
503 88 : CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
504 88 : 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 88 : 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 88 : nrow_global=nrow_global, ncol_global=nmo_occ)
509 88 : CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
510 88 : CALL cp_fm_struct_release(fm_struct_tmp)
511 88 : IF (nspins == 2) THEN
512 20 : nel = nelec_inactive/2
513 : ELSE
514 68 : nel = nelec_inactive
515 : END IF
516 : CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo_occ, nel, &
517 88 : REAL(nel, KIND=dp), maxocc, 0.0_dp)
518 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
519 88 : nrow_global=nrow_global, ncol_global=nmo_occ)
520 88 : CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO")
521 254 : CALL cp_fm_struct_release(fm_struct_tmp)
522 : END DO
523 :
524 : ! create canonical orbitals
525 78 : CALL get_qs_env(qs_env, scf_control=scf_control)
526 78 : 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 0 : "general ROKS MO definitions are not implemented.")
530 : ELSE
531 78 : IF (dft_control%do_admm) THEN
532 0 : IF (dft_control%do_admm_mo) THEN
533 0 : CPABORT("ADMM currently possible only with purification none_dm")
534 : END IF
535 : END IF
536 :
537 312 : ALLOCATE (eigenvalues(nmo_occ, nspins))
538 78 : eigenvalues = 0.0_dp
539 78 : IF (do_kpoints) THEN
540 : CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix_kp, matrix_s_kp=s_matrix_kp, &
541 2 : scf_control=scf_control)
542 2 : ks_matrix => ks_matrix_kp(:, 1)
543 2 : s_matrix => s_matrix_kp(:, 1)
544 : ELSE
545 76 : 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 78 : IF (iw > 0) THEN
550 39 : WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
551 : END IF
552 166 : 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 88 : .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 88 : 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 88 : nmo_virtual = nmo_occ - nmo_available
563 88 : nmo_virtual = MAX(nmo_virtual, 0)
564 :
565 : NULLIFY (evals_virtual)
566 178 : ALLOCATE (evals_virtual(nmo_virtual))
567 :
568 : CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
569 88 : nrow_global=nrow_global)
570 :
571 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
572 88 : nrow_global=nrow_global, ncol_global=nmo_virtual)
573 88 : CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
574 88 : CALL cp_fm_struct_release(fm_struct_tmp)
575 88 : 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 88 : 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 88 : size_ortho_space=nmo_available)
586 88 : IF (ASSOCIATED(local_preconditioner)) THEN
587 2 : CALL destroy_preconditioner(local_preconditioner)
588 2 : DEALLOCATE (local_preconditioner)
589 : END IF
590 :
591 : ! get the eigenvalues
592 88 : 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 88 : eigenvalues(:, ispin), max_offdiag)
597 :
598 : ! copy inactive orbitals
599 88 : mo_set => active_space_env%mos_inactive(ispin)
600 88 : CALL get_mo_set(mo_set, mo_coeff=mo_target)
601 160 : DO i = 1, SIZE(active_space_env%inactive_orbitals, 1)
602 72 : m = active_space_env%inactive_orbitals(i, ispin)
603 72 : CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
604 72 : mo_set%eigenvalues(m) = eigenvalues(m, ispin)
605 160 : IF (nspins > 1) THEN
606 28 : mo_set%occupation_numbers(m) = 1.0_dp
607 : ELSE
608 44 : mo_set%occupation_numbers(m) = 2.0_dp
609 : END IF
610 : END DO
611 :
612 : ! copy active orbitals
613 88 : mo_set => active_space_env%mos_active(ispin)
614 88 : CALL get_mo_set(mo_set, mo_coeff=mo_target)
615 : ! for mult > 1, put the polarized electrons in the alpha channel
616 88 : IF (nspins == 2) THEN
617 20 : IF (ispin == 1) THEN
618 10 : nel = (nelec_active + active_space_env%multiplicity - 1)/2
619 : ELSE
620 10 : nel = (nelec_active - active_space_env%multiplicity + 1)/2
621 : END IF
622 : ELSE
623 68 : nel = nelec_active
624 : END IF
625 88 : mo_set%nelectron = nel
626 88 : mo_set%n_el_f = REAL(nel, KIND=dp)
627 354 : DO i = 1, nmo_active
628 266 : m = active_space_env%active_orbitals(i, ispin)
629 266 : IF (m > nmo_available) THEN
630 8 : CALL cp_fm_to_fm(mo_virtual, mo_target, 1, m - nmo_available, m)
631 8 : eigenvalues(m, ispin) = evals_virtual(m - nmo_available)
632 8 : mo_set%occupation_numbers(m) = 0.0_dp
633 : ELSE
634 258 : CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
635 258 : mo_set%occupation_numbers(m) = mos(ispin)%occupation_numbers(m)
636 : END IF
637 354 : mo_set%eigenvalues(m) = eigenvalues(m, ispin)
638 : END DO
639 : ! Release
640 88 : DEALLOCATE (evals_virtual)
641 430 : CALL cp_fm_release(mo_virtual)
642 : END DO
643 :
644 78 : IF (iw > 0) THEN
645 83 : DO ispin = 1, nspins
646 44 : WRITE (iw, '(/,T3,A,I3,T66,A)') "Canonical Orbital Selection for spin", ispin, &
647 88 : "[atomic units]"
648 60 : DO i = 1, nmo_inactive, 4
649 16 : jm = MIN(3, nmo_inactive - i)
650 96 : WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [I]", j=0, jm)
651 : END DO
652 93 : DO i = nmo_inactive + 1, nmo_inactive + nmo_active, 4
653 49 : jm = MIN(3, nmo_inactive + nmo_active - i)
654 226 : WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [A]", j=0, jm)
655 : END DO
656 44 : WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
657 132 : DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
658 49 : jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i)
659 226 : 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 78 : DEALLOCATE (eigenvalues)
664 : END IF
665 :
666 : CASE (manual_selection)
667 : ! create canonical orbitals
668 30 : IF (dft_control%roks) THEN
669 : CALL cp_abort(__LOCATION__, &
670 : "Manual ACTIVE_SPACE orbital selection is not supported for ROKS; "// &
671 0 : "use canonical high-spin ROKS.")
672 : ELSE
673 30 : 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 0 : IF (dft_control%do_admm_mo) THEN
679 0 : CPABORT("ADMM currently possible only with purification none_dm")
680 : END IF
681 : END IF
682 :
683 30 : CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", explicit=explicit, i_vals=invals)
684 30 : IF (.NOT. explicit) THEN
685 : CALL cp_abort(__LOCATION__, "Manual orbital selection requires to explicitly "// &
686 0 : "set the active orbital indices via ACTIVE_ORBITAL_INDICES")
687 : END IF
688 :
689 30 : IF (nspins == 1) THEN
690 22 : CPASSERT(SIZE(invals) == nmo_active)
691 : ELSE
692 8 : CPASSERT(SIZE(invals) == 2*nmo_active)
693 : END IF
694 104 : ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
695 120 : ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
696 :
697 68 : DO ispin = 1, nspins
698 180 : DO i = 1, nmo_active
699 150 : active_space_env%active_orbitals(i, ispin) = invals(i + (ispin - 1)*nmo_active)
700 : END DO
701 : END DO
702 :
703 30 : IF (do_kpoints) THEN
704 0 : mos => kpoints%kp_env(1)%kpoint_env%mos(:)
705 : ELSE
706 30 : CALL get_qs_env(qs_env, mos=mos)
707 : END IF
708 :
709 : ! include MOs up to the largest index in the list
710 142 : max_orb_ind = MAXVAL(invals)
711 30 : maxocc = 2.0_dp
712 30 : IF (nspins > 1) maxocc = 1.0_dp
713 128 : ALLOCATE (active_space_env%mos_active(nspins))
714 98 : ALLOCATE (active_space_env%mos_inactive(nspins))
715 68 : DO ispin = 1, nspins
716 : ! init active orbitals
717 38 : CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
718 38 : CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
719 38 : 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 38 : nrow_global=nrow_global, ncol_global=max_orb_ind)
722 38 : CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
723 38 : CALL cp_fm_struct_release(fm_struct_tmp)
724 :
725 : ! init inactive orbitals
726 38 : IF (nspins == 2) THEN
727 16 : nel = nelec_inactive/2
728 : ELSE
729 22 : nel = nelec_inactive
730 : END IF
731 38 : 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 38 : nrow_global=nrow_global, ncol_global=max_orb_ind)
734 38 : 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 186 : active_space_env%mos_inactive(ispin)%occupation_numbers = 0.0_dp
737 106 : CALL cp_fm_struct_release(fm_struct_tmp)
738 : END DO
739 :
740 120 : ALLOCATE (eigenvalues(max_orb_ind, nspins))
741 30 : eigenvalues = 0.0_dp
742 30 : IF (do_kpoints) THEN
743 : CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix_kp, matrix_s_kp=s_matrix_kp, &
744 0 : scf_control=scf_control)
745 0 : ks_matrix => ks_matrix_kp(:, 1)
746 0 : s_matrix => s_matrix_kp(:, 1)
747 : ELSE
748 30 : 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 30 : IF (iw > 0) THEN
753 15 : WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
754 : END IF
755 68 : 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 38 : .TRUE., do_kpoints, iw)
759 38 : CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
760 38 : nmo_virtual = max_orb_ind - nmo_available
761 38 : nmo_virtual = MAX(nmo_virtual, 0)
762 :
763 : NULLIFY (evals_virtual)
764 78 : ALLOCATE (evals_virtual(nmo_virtual))
765 :
766 : CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
767 38 : nrow_global=nrow_global)
768 :
769 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
770 38 : nrow_global=nrow_global, ncol_global=nmo_virtual)
771 38 : CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
772 38 : CALL cp_fm_struct_release(fm_struct_tmp)
773 38 : CALL cp_fm_init_random(mo_virtual, nmo_virtual)
774 :
775 : CALL virtual_orbital_preconditioner(local_preconditioner, mos(ispin), ks_matrix(ispin)%matrix, &
776 38 : 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 38 : size_ortho_space=nmo_available)
783 38 : IF (ASSOCIATED(local_preconditioner)) THEN
784 2 : CALL destroy_preconditioner(local_preconditioner)
785 2 : DEALLOCATE (local_preconditioner)
786 : END IF
787 :
788 : CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, &
789 38 : 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 38 : eigenvalues(:, ispin), max_offdiag)
794 :
795 38 : mo_set_active => active_space_env%mos_active(ispin)
796 38 : CALL get_mo_set(mo_set_active, mo_coeff=fm_target_active)
797 38 : mo_set_inactive => active_space_env%mos_inactive(ispin)
798 38 : CALL get_mo_set(mo_set_inactive, mo_coeff=fm_target_inactive)
799 :
800 : ! copy orbitals
801 38 : nmo_inactive_remaining = nmo_inactive
802 186 : DO i = 1, max_orb_ind
803 : ! case for i being an active orbital
804 454 : IF (ANY(active_space_env%active_orbitals(:, ispin) == i)) THEN
805 112 : IF (i > nmo_available) THEN
806 8 : CALL cp_fm_to_fm(mo_virtual, fm_target_active, 1, i - nmo_available, i)
807 8 : eigenvalues(i, ispin) = evals_virtual(i - nmo_available)
808 8 : mo_set_active%occupation_numbers(i) = 0.0_dp
809 : ELSE
810 104 : CALL cp_fm_to_fm(mo_ref, fm_target_active, 1, i, i)
811 104 : mo_set_active%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
812 : END IF
813 112 : mo_set_active%eigenvalues(i) = eigenvalues(i, ispin)
814 : ! if it was not an active orbital, check whether it is an inactive orbital
815 36 : ELSE IF (nmo_inactive_remaining > 0) THEN
816 20 : CALL cp_fm_to_fm(mo_ref, fm_target_inactive, 1, i, i)
817 : ! store on the fly the mapping of inactive orbitals
818 20 : active_space_env%inactive_orbitals(nmo_inactive - nmo_inactive_remaining + 1, ispin) = i
819 20 : mo_set_inactive%eigenvalues(i) = eigenvalues(i, ispin)
820 20 : mo_set_inactive%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
821 : ! hack: set homo and lumo manually
822 20 : IF (nmo_inactive_remaining == 1) THEN
823 8 : mo_set_inactive%homo = i
824 8 : mo_set_inactive%lfomo = i + 1
825 : END IF
826 20 : nmo_inactive_remaining = nmo_inactive_remaining - 1
827 : ELSE
828 16 : CYCLE
829 : END IF
830 : END DO
831 :
832 : ! Release
833 38 : DEALLOCATE (evals_virtual)
834 182 : CALL cp_fm_release(mo_virtual)
835 : END DO
836 :
837 30 : IF (iw > 0) THEN
838 34 : DO ispin = 1, nspins
839 19 : WRITE (iw, '(/,T3,A,I3,T66,A)') "Orbital Energies and Selection for spin", ispin, "[atomic units]"
840 :
841 42 : DO i = 1, max_orb_ind, 4
842 23 : jm = MIN(3, max_orb_ind - i)
843 23 : WRITE (iw, '(T4)', advance="no")
844 97 : DO j = 0, jm
845 231 : IF (ANY(active_space_env%active_orbitals(:, ispin) == i + j)) THEN
846 56 : WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [A]"
847 23 : ELSE IF (ANY(active_space_env%inactive_orbitals(:, ispin) == i + j)) THEN
848 10 : WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [I]"
849 : ELSE
850 8 : WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [V]"
851 : END IF
852 : END DO
853 42 : WRITE (iw, *)
854 : END DO
855 19 : WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
856 56 : DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
857 22 : jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i)
858 97 : 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 60 : DEALLOCATE (eigenvalues)
863 : END IF
864 :
865 : CASE (wannier_projection)
866 0 : NULLIFY (loc_section, loc_print)
867 0 : loc_section => section_vals_get_subs_vals(as_input, "LOCALIZE")
868 0 : CPASSERT(ASSOCIATED(loc_section))
869 0 : loc_print => section_vals_get_subs_vals(as_input, "LOCALIZE%PRINT")
870 : !
871 0 : CPABORT("not yet available")
872 : !
873 : CASE (mao_projection)
874 : !
875 108 : CPABORT("not yet available")
876 : !
877 : END SELECT
878 :
879 : ! Warn if the active space splits a set of degenerate orbitals
880 108 : CALL check_active_space_degeneracy(active_space_env, mos)
881 :
882 : ! Print orbitals on Cube files
883 108 : print_orb => section_vals_get_subs_vals(as_input, "PRINT_ORBITAL_CUBES")
884 108 : CALL section_vals_get(print_orb, explicit=explicit)
885 108 : CALL section_vals_val_get(print_orb, "STOP_AFTER_CUBES", l_val=stop_after_print)
886 108 : IF (explicit) THEN
887 : !
888 4 : CALL print_orbital_cubes(print_orb, qs_env, active_space_env%mos_active)
889 : !
890 4 : IF (stop_after_print) THEN
891 :
892 0 : IF (iw > 0) THEN
893 : WRITE (iw, '(/,T2,A)') &
894 0 : '!----------------- Early End of Active Space Interface -----------------------!'
895 : END IF
896 :
897 0 : CALL timestop(handle)
898 :
899 0 : RETURN
900 : END IF
901 : END IF
902 :
903 : ! calculate inactive density matrix
904 108 : CALL get_qs_env(qs_env, rho=rho)
905 108 : IF (do_kpoints) THEN
906 2 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
907 2 : rho_ao => rho_ao_kp(:, 1)
908 : ELSE
909 106 : CALL qs_rho_get(rho, rho_ao=rho_ao)
910 : END IF
911 108 : CPASSERT(ASSOCIATED(rho_ao))
912 108 : CALL dbcsr_allocate_matrix_set(active_space_env%pmat_inactive, nspins)
913 234 : DO ispin = 1, nspins
914 126 : ALLOCATE (denmat)
915 126 : CALL dbcsr_copy(denmat, rho_ao(ispin)%matrix)
916 126 : mo_set => active_space_env%mos_inactive(ispin)
917 126 : CALL calculate_density_matrix(mo_set, denmat)
918 234 : active_space_env%pmat_inactive(ispin)%matrix => denmat
919 : END DO
920 :
921 : ! read in ERI parameters
922 108 : CALL section_vals_val_get(as_input, "ERI%METHOD", i_val=eri_method)
923 108 : active_space_env%eri%method = eri_method
924 108 : CALL section_vals_val_get(as_input, "ERI%OPERATOR", i_val=eri_operator, explicit=ex_operator)
925 108 : active_space_env%eri%operator = eri_operator
926 108 : CALL section_vals_val_get(as_input, "ERI%OMEGA", r_val=eri_op_omega, explicit=ex_omega)
927 108 : active_space_env%eri%omega = eri_op_omega
928 108 : CALL section_vals_val_get(as_input, "ERI%CUTOFF_RADIUS", r_val=eri_rcut, explicit=ex_rcut)
929 108 : active_space_env%eri%cutoff_radius = eri_rcut ! this is already converted to bohr!
930 108 : CALL section_vals_val_get(as_input, "ERI%PERIODICITY", i_vals=invals, explicit=ex_perd)
931 108 : CALL section_vals_val_get(as_input, "ERI%EPS_INTEGRAL", r_val=eri_eps_int)
932 108 : 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 108 : IF (ex_perd) THEN
935 96 : IF (SIZE(invals) == 1) THEN
936 0 : active_space_env%eri%periodicity(1:3) = invals(1)
937 : ELSE
938 672 : active_space_env%eri%periodicity(1:3) = invals(1:3)
939 : END IF
940 : ELSE
941 12 : CALL get_qs_env(qs_env, cell=cell)
942 84 : active_space_env%eri%periodicity(1:3) = cell%perd(1:3)
943 : END IF
944 108 : IF (iw > 0) THEN
945 54 : WRITE (iw, '(/,T3,A)') "Calculation of Electron Repulsion Integrals"
946 :
947 46 : SELECT CASE (eri_method)
948 : CASE (eri_method_full_gpw)
949 46 : WRITE (iw, '(T3,A,T50,A)') "Integration method", "GPW Fourier transform over MOs"
950 : CASE (eri_method_gpw_ht)
951 8 : WRITE (iw, '(T3,A,T44,A)') "Integration method", "Half transformed integrals from GPW"
952 : CASE DEFAULT
953 54 : CPABORT("Unknown ERI method")
954 : END SELECT
955 :
956 44 : SELECT CASE (eri_operator)
957 : CASE (eri_operator_coulomb)
958 44 : WRITE (iw, '(T3,A,T73,A)') "ERI operator", "Coulomb"
959 :
960 : CASE (eri_operator_yukawa)
961 0 : WRITE (iw, '(T3,A,T74,A)') "ERI operator", "Yukawa"
962 0 : IF (.NOT. ex_omega) CALL cp_abort(__LOCATION__, &
963 0 : "Yukawa operator requires OMEGA to be explicitly set")
964 0 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
965 :
966 : CASE (eri_operator_erf)
967 8 : WRITE (iw, '(T3,A,T63,A)') "ERI operator", "Longrange Coulomb"
968 8 : IF (.NOT. ex_omega) CALL cp_abort(__LOCATION__, &
969 0 : "Longrange operator requires OMEGA to be explicitly set")
970 8 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
971 :
972 : CASE (eri_operator_erfc)
973 0 : WRITE (iw, '(T3,A,T62,A)') "ERI operator", "Shortrange Coulomb"
974 0 : IF (.NOT. ex_omega) CALL cp_abort(__LOCATION__, &
975 0 : "Shortrange operator requires OMEGA to be explicitly set")
976 0 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
977 :
978 : CASE (eri_operator_trunc)
979 0 : WRITE (iw, '(T3,A,T63,A)') "ERI operator", "Truncated Coulomb"
980 0 : IF (.NOT. ex_rcut) CALL cp_abort(__LOCATION__, &
981 0 : "Cutoff radius not specified for trunc. Coulomb operator")
982 0 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator cutoff radius (au)", eri_rcut
983 :
984 : CASE (eri_operator_lr_trunc)
985 2 : WRITE (iw, '(T3,A,T53,A)') "ERI operator", "Longrange truncated Coulomb"
986 2 : IF (.NOT. ex_rcut) CALL cp_abort(__LOCATION__, &
987 0 : "Cutoff radius not specified for trunc. longrange operator")
988 2 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator cutoff radius (au)", eri_rcut
989 2 : IF (.NOT. ex_omega) CALL cp_abort(__LOCATION__, &
990 0 : "LR truncated operator requires OMEGA to be explicitly set")
991 2 : WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter OMEGA", eri_op_omega
992 2 : IF (eri_op_omega < 0.01_dp) THEN
993 0 : CPABORT("LR truncated operator requires OMEGA >= 0.01 to be stable")
994 : END IF
995 :
996 : CASE DEFAULT
997 54 : CPABORT("Unknown ERI operator")
998 :
999 : END SELECT
1000 :
1001 54 : WRITE (iw, '(T3,A,T68,E12.4)') "Accuracy of ERIs", eri_eps_int
1002 216 : 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 54 : IF (nspins < 2) THEN
1006 45 : WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI", (nmo_active**4)/8
1007 : ELSE
1008 9 : WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|aa)", (nmo_active**4)/8
1009 9 : WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (bb|bb)", (nmo_active**4)/8
1010 9 : 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 108 : CALL get_qs_env(qs_env, para_env=para_env)
1016 108 : m = (nspins*(nspins + 1))/2
1017 : ! With ROHF/ROKS, we need ERIs from only a single set of orbitals
1018 108 : IF (dft_control%roks) m = 1
1019 460 : ALLOCATE (active_space_env%eri%eri(m))
1020 244 : DO i = 1, m
1021 136 : CALL get_mo_set(active_space_env%mos_active(1), nmo=nmo)
1022 136 : ALLOCATE (active_space_env%eri%eri(i)%csr_mat)
1023 136 : eri_mat => active_space_env%eri%eri(i)%csr_mat
1024 136 : IF (i == 1) THEN
1025 108 : n1 = nmo
1026 108 : n2 = nmo
1027 28 : ELSE IF (i == 2) THEN
1028 14 : n1 = nmo
1029 14 : n2 = nmo
1030 : ELSE
1031 14 : n1 = nmo
1032 14 : n2 = nmo
1033 : END IF
1034 136 : nn1 = (n1*(n1 + 1))/2
1035 136 : nn2 = (n2*(n2 + 1))/2
1036 136 : CALL dbcsr_csr_create(eri_mat, nn1, nn2, 0_int_8, 0, 0, para_env%get_handle())
1037 380 : active_space_env%eri%norb = nmo
1038 : END DO
1039 :
1040 108 : SELECT CASE (eri_method)
1041 : CASE (eri_method_full_gpw, eri_method_gpw_ht)
1042 108 : CALL section_vals_val_get(as_input, "ERI_GPW%EPS_GRID", r_val=eri_eps_grid)
1043 108 : active_space_env%eri%eri_gpw%eps_grid = eri_eps_grid
1044 108 : CALL section_vals_val_get(as_input, "ERI_GPW%EPS_FILTER", r_val=eri_eps_filter)
1045 108 : active_space_env%eri%eri_gpw%eps_filter = eri_eps_filter
1046 108 : CALL section_vals_val_get(as_input, "ERI_GPW%CUTOFF", r_val=eri_gpw_cutoff)
1047 108 : active_space_env%eri%eri_gpw%cutoff = eri_gpw_cutoff
1048 108 : CALL section_vals_val_get(as_input, "ERI_GPW%REL_CUTOFF", r_val=eri_rel_cutoff)
1049 108 : active_space_env%eri%eri_gpw%rel_cutoff = eri_rel_cutoff
1050 108 : CALL section_vals_val_get(as_input, "ERI_GPW%PRINT_LEVEL", i_val=eri_print)
1051 108 : active_space_env%eri%eri_gpw%print_level = eri_print
1052 108 : CALL section_vals_val_get(as_input, "ERI_GPW%STORE_WFN", l_val=store_wfn)
1053 108 : active_space_env%eri%eri_gpw%store_wfn = store_wfn
1054 108 : CALL section_vals_val_get(as_input, "ERI_GPW%GROUP_SIZE", i_val=group_size)
1055 108 : active_space_env%eri%eri_gpw%group_size = group_size
1056 : ! Always redo Poisson solver for now
1057 108 : active_space_env%eri%eri_gpw%redo_poisson = .TRUE.
1058 : ! active_space_env%eri%eri_gpw%redo_poisson = (ex_operator .OR. ex_perd)
1059 108 : IF (iw > 0) THEN
1060 54 : WRITE (iw, '(/,T2,A,T71,F10.1)') "ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff
1061 54 : 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 108 : dft_control%roks)
1066 : !
1067 : CASE DEFAULT
1068 108 : CPABORT("Unknown ERI method")
1069 : END SELECT
1070 108 : IF (iw > 0) THEN
1071 122 : DO isp = 1, SIZE(active_space_env%eri%eri)
1072 68 : eri_mat => active_space_env%eri%eri(isp)%csr_mat
1073 : nze_percentage = 100.0_dp*(REAL(eri_mat%nze_total, KIND=dp) &
1074 68 : /REAL(eri_mat%nrows_total, KIND=dp))/REAL(eri_mat%ncols_total, KIND=dp)
1075 68 : WRITE (iw, '(/,T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1076 136 : "Number of CSR non-zero elements:", eri_mat%nze_total
1077 68 : WRITE (iw, '(T2,A,I2,T30,A,T68,F12.4)') "ERI_GPW| Spinmatrix:", isp, &
1078 136 : "Percentage CSR non-zero elements:", nze_percentage
1079 68 : WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1080 136 : "nrows_total", eri_mat%nrows_total
1081 68 : WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1082 136 : "ncols_total", eri_mat%ncols_total
1083 68 : WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
1084 190 : "nrows_local", eri_mat%nrows_local
1085 : END DO
1086 54 : CALL m_flush(iw)
1087 : END IF
1088 108 : CALL para_env%sync()
1089 :
1090 : ! set the reference active space density matrix
1091 108 : nspins = active_space_env%nspins
1092 450 : ALLOCATE (active_space_env%p_active(nspins))
1093 234 : DO isp = 1, nspins
1094 126 : mo_set => active_space_env%mos_active(isp)
1095 126 : CALL get_mo_set(mo_set, mo_coeff=mo_coeff, nmo=nmo)
1096 234 : CALL create_subspace_matrix(mo_coeff, active_space_env%p_active(isp), nmo)
1097 : END DO
1098 0 : SELECT CASE (mselect)
1099 : CASE DEFAULT
1100 0 : CPABORT("Unknown orbital selection method")
1101 : CASE (casci_canonical, manual_selection)
1102 108 : focc = 2.0_dp
1103 108 : IF (nspins == 2) focc = 1.0_dp
1104 234 : DO isp = 1, nspins
1105 126 : fmat => active_space_env%p_active(isp)
1106 126 : CALL cp_fm_set_all(fmat, alpha=0.0_dp)
1107 126 : IF (nspins == 2) THEN
1108 36 : IF (isp == 1) THEN
1109 18 : nel = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2
1110 : ELSE
1111 18 : nel = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2
1112 : END IF
1113 : ELSE
1114 90 : nel = active_space_env%nelec_active
1115 : END IF
1116 612 : DO i = 1, nmo_active
1117 378 : m = active_space_env%active_orbitals(i, isp)
1118 378 : fel = MIN(focc, REAL(nel, KIND=dp))
1119 378 : CALL cp_fm_set_element(fmat, m, m, fel)
1120 378 : nel = nel - NINT(fel)
1121 504 : nel = MAX(nel, 0)
1122 : END DO
1123 : END DO
1124 : CASE (wannier_projection)
1125 0 : CPABORT("NOT IMPLEMENTED")
1126 : CASE (mao_projection)
1127 108 : CPABORT("NOT IMPLEMENTED")
1128 : END SELECT
1129 :
1130 : ! compute alpha-beta overlap matrix in case of spin-polarized calculation
1131 108 : 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 108 : xc_section => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE%XC")
1135 108 : explicit = .FALSE.
1136 108 : IF (ASSOCIATED(xc_section)) CALL section_vals_get(xc_section, explicit=explicit)
1137 :
1138 : ! rebuild KS matrix if needed
1139 108 : IF (explicit) THEN
1140 : ! release the hfx data if it was part of the SCF functional
1141 2 : 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 2 : 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 2 : particle_set=particle_set, cell=cell, ks_env=ks_env)
1147 2 : IF (dft_control%do_admm) THEN
1148 0 : basis_type = 'AUX_FIT'
1149 : ELSE
1150 2 : basis_type = 'ORB'
1151 : END IF
1152 2 : 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 2 : nelectron_total=nelec_total)
1156 :
1157 2 : 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 2 : 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 2 : ext_xc_section=xc_section)
1164 : ! we need to reset it to false
1165 2 : CALL set_ks_env(ks_env, s_mstruct_changed=.FALSE.)
1166 : ELSE
1167 106 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1168 : END IF
1169 : ! set the xc_section
1170 108 : active_space_env%xc_section => xc_section
1171 :
1172 108 : CALL get_qs_env(qs_env, energy=energy)
1173 : ! transform KS/Fock, Vxc and Hcore to AS MO basis
1174 108 : CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env)
1175 : ! set the reference energy in the active space
1176 108 : active_space_env%energy_ref = energy%total
1177 : ! calculate inactive energy and embedding potential
1178 108 : CALL subspace_fock_matrix(active_space_env, dft_control%roks)
1179 :
1180 : ! associate the active space environment with the qs environment
1181 108 : 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 108 : CALL section_vals_val_get(as_input, "AS_SOLVER", i_val=as_solver)
1185 102 : SELECT CASE (as_solver)
1186 : CASE (no_solver)
1187 102 : IF (iw > 0) THEN
1188 51 : WRITE (iw, '(/,T3,A)') "No active space solver specified, skipping embedding calculation"
1189 51 : CALL m_flush(iw)
1190 : END IF
1191 102 : CALL para_env%sync()
1192 : CASE (qiskit_solver)
1193 0 : CALL rsdft_embedding(qs_env, active_space_env, as_input)
1194 0 : CALL qs_scf_compute_properties(qs_env, wf_type="MC-DFT", do_mp2=.FALSE.)
1195 : CASE (fci_solver)
1196 6 : CALL local_fci_embedding(qs_env, active_space_env, as_input)
1197 6 : CALL qs_scf_compute_properties(qs_env, wf_type="MC-DFT", do_mp2=.FALSE.)
1198 : CASE DEFAULT
1199 108 : CPABORT("Unknown active space solver")
1200 : END SELECT
1201 :
1202 : ! Output a FCIDUMP file if requested
1203 108 : 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 108 : IF (BTEST(cp_print_key_should_output(logger%iter_info, as_input, "MOLDEN"), cp_p_file)) THEN
1207 2 : CALL write_active_space_molden(qs_env, active_space_env, as_input)
1208 : END IF
1209 :
1210 : ! Output a QCSchema file if requested
1211 108 : IF (active_space_env%qcschema) THEN
1212 6 : CALL qcschema_env_create(qcschema_env, qs_env)
1213 6 : CALL qcschema_to_hdf5(qcschema_env, active_space_env%qcschema_filename)
1214 6 : CALL qcschema_env_release(qcschema_env)
1215 : END IF
1216 :
1217 108 : IF (iw > 0) THEN
1218 : WRITE (iw, '(/,T2,A)') &
1219 54 : '!-------------------- End of Active Space Interface --------------------------!'
1220 54 : CALL m_flush(iw)
1221 : END IF
1222 108 : CALL para_env%sync()
1223 :
1224 108 : CALL timestop(handle)
1225 :
1226 117312 : 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 126 : 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 126 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: energies
1255 126 : 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 116 : IF (do_kpoints) RETURN
1264 124 : 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 124 : mo_derivs=mo_derivs)
1267 124 : IF (dft_control%roks) RETURN
1268 :
1269 : ! a converged diagonalization SCF leaves off-diagonal elements of the order of EPS_SCF
1270 116 : eps_canonical = MAX(1.0E-4_dp, 100.0_dp*scf_control%eps_scf)
1271 116 : IF (mo_set%use_mo_coeff_b) CALL copy_dbcsr_to_fm(mo_set%mo_coeff_b, mo_set%mo_coeff)
1272 348 : ALLOCATE (energies(mo_set%nmo))
1273 116 : CALL calculate_orbital_fock_elements(mo_set%mo_coeff, ks_matrix, energies, max_offdiag)
1274 116 : DEALLOCATE (energies)
1275 116 : IF (max_offdiag <= eps_canonical) THEN
1276 156 : IF (iw > 0) WRITE (iw, '(/,T3,A,I0,A,ES9.2,A)') "SCF orbitals of spin ", ispin, &
1277 104 : " are canonical (max |<i|F|j>| = ", max_offdiag, ")"
1278 104 : RETURN
1279 : END IF
1280 :
1281 12 : on_the_fly_loc = .FALSE.
1282 12 : 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 12 : IF (on_the_fly_loc) THEN
1284 2 : 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 0 : "to select localized orbitals, or switch off ON_THE_FLY_LOC.")
1289 : END IF
1290 3 : IF (iw > 0) WRITE (iw, '(/,T3,A,I0,A)') "SCF orbitals of spin ", ispin, &
1291 2 : " are localized (ON_THE_FLY_LOC): ACTIVE_ORBITAL_INDICES refer to the localized orbitals"
1292 2 : RETURN
1293 : END IF
1294 :
1295 10 : IF (ASSOCIATED(mo_derivs)) THEN
1296 10 : IF (SIZE(mo_derivs) >= ispin) mo_deriv => mo_derivs(ispin)%matrix
1297 : END IF
1298 10 : CALL make_mo_set_canonical(mo_set, ks_matrix, eps_canonical, rotated, max_offdiag, mo_deriv)
1299 10 : IF (iw > 0) THEN
1300 5 : WRITE (iw, '(/,T3,A,I0,A,ES9.2,A)') "SCF orbitals of spin ", ispin, &
1301 10 : " are not canonical (max |<i|F|j>| = ", max_offdiag, ")"
1302 5 : WRITE (iw, '(T3,A)') "Rotated into canonical orbitals, separately for occupied and virtual orbitals"
1303 : END IF
1304 :
1305 366 : 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 126 : 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 126 : NULLIFY (preconditioner)
1328 126 : IF (nmo_virtual <= 0) RETURN
1329 16 : ALLOCATE (preconditioner)
1330 4 : 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 4 : 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 108 : 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 108 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues
1355 :
1356 234 : DO ispin = 1, SIZE(active_space_env%active_orbitals, 2)
1357 126 : mspin = MIN(ispin, SIZE(mos))
1358 126 : CALL get_mo_set(mos(mspin), nmo=nmo, eigenvalues=eigenvalues)
1359 : ! no orbital energies available
1360 126 : IF (.NOT. ASSOCIATED(eigenvalues)) CYCLE
1361 134 : 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 124 : eps_degenerate, npairs, first_active, first_other)
1364 358 : IF (npairs > 0) THEN
1365 : WRITE (message, '(A,I0,A,I0,A,I0,A,I0,A)') &
1366 0 : "The active space splits a set of degenerate orbitals: active orbital ", first_active, &
1367 0 : " and inactive/virtual orbital ", first_other, " (spin ", ispin, ", ", npairs, &
1368 0 : " degenerate pairs in total)."
1369 0 : CPWARN(TRIM(message))
1370 : END IF
1371 : END DO
1372 :
1373 108 : 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 2 : 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 2 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_out
1395 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1396 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1397 : TYPE(section_vals_type), POINTER :: molden_section
1398 :
1399 2 : NULLIFY (cell, fm_struct, mo_coeff_active, mo_coeff_inactive, particle_set, qs_kind_set)
1400 2 : CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set)
1401 2 : molden_section => section_vals_get_subs_vals(as_input, "MOLDEN")
1402 :
1403 2 : nspins = SIZE(active_space_env%mos_active)
1404 2 : nmo_active = SIZE(active_space_env%active_orbitals, 1)
1405 2 : nmo_inactive = SIZE(active_space_env%inactive_orbitals, 1)
1406 2 : nmo_out = nmo_inactive + nmo_active
1407 8 : ALLOCATE (mos_out(nspins))
1408 4 : DO ispin = 1, nspins
1409 : CALL get_mo_set(active_space_env%mos_active(ispin), nao=nao, maxocc=maxocc, &
1410 2 : mo_coeff=mo_coeff_active)
1411 2 : CALL get_mo_set(active_space_env%mos_inactive(ispin), mo_coeff=mo_coeff_inactive)
1412 2 : 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 2 : template_fmstruct=mo_coeff_active%matrix_struct)
1415 2 : CALL init_mo_set(mos_out(ispin), fm_struct=fm_struct, name="active space orbitals")
1416 2 : 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 6 : DO i = 1, nmo_inactive
1419 4 : m = active_space_env%inactive_orbitals(i, ispin)
1420 4 : CALL cp_fm_to_fm(mo_coeff_inactive, mos_out(ispin)%mo_coeff, 1, m, i)
1421 4 : mos_out(ispin)%eigenvalues(i) = active_space_env%mos_inactive(ispin)%eigenvalues(m)
1422 6 : mos_out(ispin)%occupation_numbers(i) = active_space_env%mos_inactive(ispin)%occupation_numbers(m)
1423 : END DO
1424 18 : DO i = 1, nmo_active
1425 12 : m = active_space_env%active_orbitals(i, ispin)
1426 12 : iorb = nmo_inactive + i
1427 12 : CALL cp_fm_to_fm(mo_coeff_active, mos_out(ispin)%mo_coeff, 1, m, iorb)
1428 12 : mos_out(ispin)%eigenvalues(iorb) = active_space_env%mos_active(ispin)%eigenvalues(m)
1429 14 : mos_out(ispin)%occupation_numbers(iorb) = active_space_env%mos_active(ispin)%occupation_numbers(m)
1430 : END DO
1431 : END DO
1432 :
1433 2 : CALL write_mos_molden(mos_out, qs_kind_set, particle_set, molden_section, cell=cell)
1434 :
1435 4 : DO ispin = 1, nspins
1436 4 : CALL deallocate_mo_set(mos_out(ispin))
1437 : END DO
1438 2 : DEALLOCATE (mos_out)
1439 :
1440 2 : 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 108 : 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 108 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_matrix
1462 108 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: s_matrix_kp
1463 :
1464 108 : CALL timeset(routineN, handle)
1465 :
1466 108 : nspins = active_space_env%nspins
1467 :
1468 : ! overlap in AO
1469 108 : IF (nspins > 1) THEN
1470 18 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
1471 18 : IF (do_kpoints) THEN
1472 0 : CALL get_qs_env(qs_env, matrix_s_kp=s_matrix_kp)
1473 0 : s_matrix => s_matrix_kp(:, 1)
1474 : ELSE
1475 18 : CALL get_qs_env(qs_env, matrix_s=s_matrix)
1476 : END IF
1477 36 : ALLOCATE (active_space_env%sab_sub(1))
1478 :
1479 18 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff_a, nmo=nmo)
1480 18 : CALL get_mo_set(mo_set=mos(2), mo_coeff=mo_coeff_b, nmo=nmo)
1481 18 : 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 108 : CALL timestop(handle)
1485 :
1486 108 : 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 130 : 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 130 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: h_matrix, ks_matrix
1507 :
1508 130 : CALL timeset(routineN, handle)
1509 :
1510 130 : nspins = active_space_env%nspins
1511 :
1512 : ! Kohn-Sham / Fock operator
1513 130 : CALL cp_fm_release(active_space_env%ks_sub)
1514 130 : CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix)
1515 552 : ALLOCATE (active_space_env%ks_sub(nspins))
1516 292 : DO ispin = 1, nspins
1517 162 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1518 292 : 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 130 : CALL cp_fm_release(active_space_env%h_sub)
1523 :
1524 130 : NULLIFY (h_matrix)
1525 130 : CALL get_qs_env(qs_env=qs_env, matrix_h_kp=h_matrix)
1526 422 : ALLOCATE (active_space_env%h_sub(nspins))
1527 292 : DO ispin = 1, nspins
1528 162 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1529 292 : CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(ispin))
1530 : END DO
1531 :
1532 130 : CALL timestop(handle)
1533 :
1534 130 : 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 684 : 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 342 : CALL timeset(routineN, handle)
1560 :
1561 342 : CALL cp_fm_get_info(matrix=mo_coeff, ncol_global=ncol, nrow_global=nrow)
1562 342 : CPASSERT(nmo <= ncol)
1563 :
1564 342 : IF (nmo > 0) THEN
1565 342 : CALL cp_fm_create(vectors, mo_coeff%matrix_struct, "vectors")
1566 342 : CALL create_subspace_matrix(mo_coeff, op_sub, nmo)
1567 :
1568 342 : IF (PRESENT(mo_coeff_b)) THEN
1569 : ! if beta orbitals are present, compute the cross alpha_beta term
1570 18 : CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff_b, vectors, nmo)
1571 : ELSE
1572 : ! otherwise the same spin, whatever that is
1573 324 : CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff, vectors, nmo)
1574 : END IF
1575 :
1576 342 : CALL parallel_gemm('T', 'N', nmo, nmo, nrow, 1.0_dp, mo_coeff, vectors, 0.0_dp, op_sub)
1577 342 : CALL cp_fm_release(vectors)
1578 : END IF
1579 :
1580 342 : CALL timestop(handle)
1581 :
1582 342 : 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 468 : 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 468 : IF (n > 0) THEN
1601 :
1602 468 : 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 468 : context=orbitals%matrix_struct%context)
1606 468 : CALL cp_fm_create(op_sub, fm_struct, name="Subspace operator")
1607 468 : CALL cp_fm_struct_release(fm_struct)
1608 :
1609 : END IF
1610 :
1611 468 : 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 108 : 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 108 : INTEGER, ALLOCATABLE, DIMENSION(:) :: eri_index
1640 108 : 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 108 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri
1646 108 : 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 108 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_matrix_pq_rnu, fm_matrix_pq_rs, &
1651 108 : fm_mo_coeff_as
1652 : TYPE(cp_fm_type), POINTER :: mo_coeff
1653 : TYPE(dbcsr_p_type) :: mat_munu
1654 108 : 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 108 : POINTER :: sab_orb_sub
1659 108 : 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 108 : 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 108 : 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 108 : CALL timeset(routineN, handle)
1674 :
1675 108 : IF (iw > 0) t1 = m_walltime()
1676 :
1677 : ! print levels
1678 204 : SELECT CASE (eri_env%eri_gpw%print_level)
1679 : CASE (silent_print_level)
1680 96 : print1 = .FALSE.
1681 96 : print2 = .FALSE.
1682 : CASE (low_print_level)
1683 6 : print1 = .FALSE.
1684 6 : print2 = .FALSE.
1685 : CASE (medium_print_level)
1686 6 : print1 = .TRUE.
1687 6 : print2 = .FALSE.
1688 : CASE (high_print_level)
1689 0 : print1 = .TRUE.
1690 0 : print2 = .TRUE.
1691 : CASE (debug_print_level)
1692 0 : print1 = .TRUE.
1693 108 : print2 = .TRUE.
1694 : CASE DEFAULT
1695 : ! do nothing
1696 : END SELECT
1697 :
1698 : ! Check the input group
1699 108 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
1700 108 : IF (eri_env%eri_gpw%group_size < 1) eri_env%eri_gpw%group_size = para_env%num_pe
1701 108 : IF (MOD(para_env%num_pe, eri_env%eri_gpw%group_size) /= 0) THEN
1702 0 : 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 108 : IF (eri_env%eri_gpw%group_size == para_env%num_pe) THEN
1706 102 : eri_env%para_env_sub => para_env
1707 102 : CALL eri_env%para_env_sub%retain()
1708 102 : blacs_env_sub => blacs_env
1709 102 : CALL blacs_env_sub%retain()
1710 102 : number_of_subgroups = 1
1711 102 : color = 0
1712 : ELSE
1713 6 : number_of_subgroups = para_env%num_pe/eri_env%eri_gpw%group_size
1714 6 : color = para_env%mepos/eri_env%eri_gpw%group_size
1715 6 : ALLOCATE (eri_env%para_env_sub)
1716 6 : CALL eri_env%para_env_sub%from_split(para_env, color)
1717 6 : NULLIFY (blacs_env_sub)
1718 6 : CALL cp_blacs_env_create(blacs_env_sub, eri_env%para_env_sub, BLACS_GRID_SQUARE, .TRUE.)
1719 : END IF
1720 108 : 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 108 : CALL get_qs_env(qs_env, dft_control=dft_control)
1724 432 : ALLOCATE (qs_control)
1725 108 : qs_control_old => dft_control%qs_control
1726 108 : qs_control = qs_control_old
1727 108 : dft_control%qs_control => qs_control
1728 108 : progression_factor = qs_control%progression_factor
1729 108 : n_multigrid = SIZE(qs_control%e_cutoff)
1730 108 : 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 108 : IF (restricted) nspins = 1
1734 : ! Allocate new cutoffs (just in private qs_control, not in qs_control_old)
1735 324 : ALLOCATE (qs_control%e_cutoff(n_multigrid))
1736 :
1737 108 : qs_control%cutoff = eri_env%eri_gpw%cutoff*0.5_dp
1738 108 : qs_control%e_cutoff(1) = qs_control%cutoff
1739 432 : DO i_multigrid = 2, n_multigrid
1740 : qs_control%e_cutoff(i_multigrid) = qs_control%e_cutoff(i_multigrid - 1) &
1741 432 : /progression_factor
1742 : END DO
1743 108 : 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 108 : 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 108 : do_alloc_blocks_from_nbl=.TRUE., dbcsr_sym_type=dbcsr_type_symmetric)
1749 108 : CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
1750 :
1751 : ! Generate the appropriate pw_env
1752 108 : NULLIFY (pw_env_sub)
1753 108 : CALL pw_env_create(pw_env_sub)
1754 108 : CALL pw_env_rebuild(pw_env_sub, qs_env, external_para_env=eri_env%para_env_sub)
1755 108 : 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 108 : IF (eri_env%eri_gpw%redo_poisson) THEN
1759 : ! We need to rebuild the Poisson solver on the fly
1760 432 : IF (SUM(eri_env%periodicity) /= 0) THEN
1761 10 : poisson_env%parameters%solver = pw_poisson_periodic
1762 : ELSE
1763 98 : poisson_env%parameters%solver = pw_poisson_analytic
1764 : END IF
1765 432 : 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 108 : CALL pw_poisson_rebuild(poisson_env)
1771 :
1772 : ! set the cutoff radius for the Greens function in case we use ANALYTIC Poisson solver
1773 108 : CALL get_qs_env(qs_env, cell=cell)
1774 108 : rc = cell%hmat(1, 1)
1775 432 : DO iwa1 = 1, 3
1776 : ! TODO: I think this is not the largest possible radius inscribed in the cell
1777 432 : rc = MIN(rc, 0.5_dp*cell%hmat(iwa1, iwa1))
1778 : END DO
1779 108 : poisson_env%green_fft%radius = rc
1780 :
1781 : ! Overwrite the Greens function with the one we want
1782 108 : CALL pw_eri_green_create(poisson_env%green_fft, eri_env)
1783 :
1784 108 : IF (iw > 0) THEN
1785 54 : CALL get_qs_env(qs_env, cell=cell)
1786 432 : IF (SUM(cell%perd) /= SUM(eri_env%periodicity)) THEN
1787 0 : IF (SUM(eri_env%periodicity) /= 0) THEN
1788 : WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") &
1789 0 : "ERI_GPW| Switching Poisson solver to", "PERIODIC"
1790 : ELSE
1791 : WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") &
1792 0 : "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 59 : SELECT CASE (poisson_env%green_fft%method)
1797 : CASE (PERIODIC3D)
1798 : WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
1799 5 : "ERI_GPW| Poisson Greens function", "PERIODIC"
1800 : CASE (ANALYTIC0D)
1801 : WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
1802 49 : "ERI_GPW| Poisson Greens function", "ANALYTIC"
1803 49 : WRITE (UNIT=iw, FMT="(T2,A,T71,F10.4)") "ERI_GPW| Poisson cutoff radius", &
1804 98 : poisson_env%green_fft%radius*angstrom
1805 : CASE DEFAULT
1806 54 : CPABORT("Wrong Greens function setup")
1807 : END SELECT
1808 : END IF
1809 : END IF
1810 :
1811 784 : ALLOCATE (mo_coeff_as(nspins), fm_mo_coeff_as(nspins))
1812 230 : DO ispin = 1, nspins
1813 366 : BLOCK
1814 122 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: C, C_active
1815 : INTEGER :: nmo
1816 122 : TYPE(group_dist_d1_type) :: gd_array
1817 : TYPE(cp_fm_type), POINTER :: mo_coeff
1818 122 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1819 122 : CALL grep_rows_in_subgroups(para_env, eri_env%para_env_sub, mo_coeff, gd_array, C)
1820 :
1821 488 : ALLOCATE (C_active(SIZE(C, 1), SIZE(orbitals, 1)))
1822 488 : DO i1 = 1, SIZE(orbitals, 1)
1823 2386 : 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 122 : C_active, mat_munu%matrix, gd_array, eri_env%eri_gpw%eps_filter)
1827 122 : CALL release_group_dist(gd_array)
1828 488 : DEALLOCATE (C, C_active)
1829 : END BLOCK
1830 :
1831 122 : CALL dbcsr_get_info(mo_coeff_as(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
1832 :
1833 122 : NULLIFY (fm_struct)
1834 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1835 122 : nrow_global=nrow_global, ncol_global=ncol_global)
1836 122 : CALL cp_fm_create(fm_mo_coeff_as(ispin), fm_struct)
1837 122 : CALL cp_fm_struct_release(fm_struct)
1838 :
1839 352 : CALL copy_dbcsr_to_fm(mo_coeff_as(ispin), fm_mo_coeff_as(ispin))
1840 : END DO
1841 :
1842 108 : IF (eri_env%method == eri_method_gpw_ht) THEN
1843 : ! We need a task list
1844 16 : NULLIFY (task_list_sub)
1845 16 : skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
1846 16 : 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 16 : 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 112 : ALLOCATE (matrix_pq_rnu(nspins), fm_matrix_pq_rnu(nspins), fm_matrix_pq_rs(nspins))
1856 32 : DO ispin = 1, nspins
1857 16 : CALL dbcsr_create(matrix_pq_rnu(ispin), template=mo_coeff_as(ispin))
1858 16 : CALL dbcsr_set(matrix_pq_rnu(ispin), 0.0_dp)
1859 :
1860 16 : CALL dbcsr_get_info(matrix_pq_rnu(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
1861 :
1862 16 : NULLIFY (fm_struct)
1863 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1864 16 : nrow_global=nrow_global, ncol_global=ncol_global)
1865 16 : CALL cp_fm_create(fm_matrix_pq_rnu(ispin), fm_struct)
1866 16 : CALL cp_fm_struct_release(fm_struct)
1867 :
1868 16 : NULLIFY (fm_struct)
1869 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, para_env=eri_env%para_env_sub, &
1870 16 : nrow_global=ncol_global, ncol_global=ncol_global)
1871 16 : CALL cp_fm_create(fm_matrix_pq_rs(ispin), fm_struct)
1872 48 : 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 108 : CALL auxbas_pw_pool%create_pw(wfn_r)
1879 108 : CALL auxbas_pw_pool%create_pw(rho_g)
1880 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, &
1881 108 : particle_set=particle_set, atomic_kind_set=atomic_kind_set)
1882 :
1883 : ! pre-calculate wavefunctions on reals space grid
1884 108 : 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 108 : IF (eri_env%eri_gpw%store_wfn) THEN
1889 : ! pre-calculate wavefunctions on reals space grid
1890 96 : rsize = 0.0_dp
1891 96 : nmo = 0
1892 206 : DO ispin = 1, nspins
1893 110 : CALL get_mo_set(mo_set=mos(ispin), nmo=nx)
1894 110 : nmo = MAX(nmo, nx)
1895 536 : rsize = REAL(SIZE(wfn_r%array), KIND=dp)*nx
1896 : END DO
1897 96 : IF (print1 .AND. iw > 0) THEN
1898 3 : rsize = rsize*8._dp/1000000._dp
1899 3 : WRITE (iw, "(T2,'ERI_GPW|',' Store active orbitals on real space grid ',T66,F12.3,' MB')") rsize
1900 : END IF
1901 944 : ALLOCATE (wfn_a(nmo, nspins))
1902 206 : DO ispin = 1, nspins
1903 110 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1904 548 : DO i1 = 1, SIZE(orbitals, 1)
1905 342 : iwfn = orbitals(i1, ispin)
1906 342 : 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 342 : qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1909 452 : IF (print2 .AND. iw > 0) THEN
1910 0 : 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 12 : ALLOCATE (wfn1, wfn2)
1917 12 : CALL auxbas_pw_pool%create_pw(wfn1)
1918 12 : CALL auxbas_pw_pool%create_pw(wfn2)
1919 12 : IF (eri_env%method /= eri_method_gpw_ht) THEN
1920 6 : ALLOCATE (wfn3, wfn4)
1921 6 : CALL auxbas_pw_pool%create_pw(wfn3)
1922 6 : CALL auxbas_pw_pool%create_pw(wfn4)
1923 : END IF
1924 : END IF
1925 :
1926 : ! get some of the grids ready
1927 108 : CALL auxbas_pw_pool%create_pw(rho_r)
1928 108 : 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 108 : CALL pw_zero(rho_r)
1932 108 : CALL pw_transfer(rho_r, rho_g)
1933 108 : dvol = rho_r%pw_grid%dvol
1934 :
1935 108 : IF (iw > 0) THEN
1936 54 : CALL m_flush(iw)
1937 : END IF
1938 : ! calculate the integrals
1939 108 : stored_integrals = 0
1940 230 : DO isp1 = 1, nspins
1941 122 : CALL get_mo_set(mo_set=mos(isp1), nmo=nmo1)
1942 122 : nmm = (nmo1*(nmo1 + 1))/2
1943 596 : DO i1 = 1, SIZE(orbitals, 1)
1944 366 : iwa1 = orbitals(i1, isp1)
1945 366 : IF (eri_env%eri_gpw%store_wfn) THEN
1946 342 : wfn1 => wfn_a(iwa1, isp1)
1947 : ELSE
1948 : CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwa1, wfn1, rho_g, atomic_kind_set, &
1949 24 : qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1950 : END IF
1951 1334 : DO i2 = i1, SIZE(orbitals, 1)
1952 846 : iwa2 = orbitals(i2, isp1)
1953 846 : iwa12 = csr_idx_to_combined(iwa1, iwa2, nmo1)
1954 : ! Skip calculation directly if the pair is not part of our subgroup
1955 846 : IF (MOD(iwa12 - 1, eri_env%comm_exchange%num_pe) /= eri_env%comm_exchange%mepos) CYCLE
1956 837 : iwa12 = (iwa12 - 1)/eri_env%comm_exchange%num_pe + 1
1957 837 : IF (eri_env%eri_gpw%store_wfn) THEN
1958 807 : wfn2 => wfn_a(iwa2, isp1)
1959 : ELSE
1960 : CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwa2, wfn2, rho_g, atomic_kind_set, &
1961 30 : qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
1962 : END IF
1963 : ! calculate charge distribution and potential
1964 837 : CALL pw_zero(rho_r)
1965 837 : CALL pw_multiply(rho_r, wfn1, wfn2)
1966 837 : CALL pw_transfer(rho_r, rho_g)
1967 837 : CALL pw_poisson_solve(poisson_env, rho_g, pair_int, pot_g)
1968 :
1969 : ! screening using pair_int
1970 837 : IF (pair_int < eri_env%eps_integral) CYCLE
1971 837 : CALL pw_transfer(pot_g, rho_r)
1972 : !
1973 2040 : IF (eri_env%method == eri_method_gpw_ht) THEN
1974 56 : CALL pw_scale(rho_r, dvol)
1975 112 : DO isp2 = isp1, nspins
1976 56 : CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2)
1977 56 : nx = (nmo2*(nmo2 + 1))/2
1978 280 : ALLOCATE (eri(nx), eri_index(nx))
1979 56 : 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 56 : 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 56 : 0.0_dp, matrix_pq_rnu(isp2), filter_eps=eri_env%eri_gpw%eps_filter)
1986 56 : CALL copy_dbcsr_to_fm(matrix_pq_rnu(isp2), fm_matrix_pq_rnu(isp2))
1987 :
1988 56 : 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 56 : 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 56 : 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 56 : col_indices=col_indices, row_indices=row_indices)
1999 :
2000 56 : icount2 = 0
2001 208 : DO col_local = 1, ncol_local
2002 152 : iwb1 = orbitals(col_indices(col_local), isp2)
2003 152 : IF (isp1 == isp2 .AND. iwb1 < iwa1) CYCLE
2004 366 : DO row_local = 1, nrow_local
2005 190 : iwb2 = orbitals(row_indices(row_local), isp2)
2006 190 : IF (iwb2 < iwb1) CYCLE
2007 114 : IF (isp1 == isp2 .AND. iwa1 == iwb1 .AND. iwb2 < iwa2) CYCLE
2008 :
2009 97 : iwb12 = csr_idx_to_combined(iwb1, iwb2, nmo2)
2010 97 : erint = fm_matrix_pq_rs(isp2)%local_data(row_local, col_local)
2011 249 : IF (ABS(erint) > eri_env%eps_integral) THEN
2012 90 : icount2 = icount2 + 1
2013 90 : eri(icount2) = erint
2014 90 : eri_index(icount2) = iwb12
2015 : END IF
2016 : END DO
2017 : END DO
2018 56 : stored_integrals = stored_integrals + icount2
2019 : !
2020 56 : isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
2021 56 : CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
2022 : !
2023 280 : DEALLOCATE (eri, eri_index)
2024 : END DO
2025 781 : ELSE IF (eri_env%method == eri_method_full_gpw) THEN
2026 1636 : DO isp2 = isp1, nspins
2027 855 : CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2)
2028 855 : nx = (nmo2*(nmo2 + 1))/2
2029 4275 : ALLOCATE (eri(nx), eri_index(nx))
2030 855 : icount2 = 0
2031 855 : iwbs = 1
2032 855 : IF (isp1 == isp2) iwbs = i1
2033 855 : isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
2034 3578 : DO i3 = iwbs, SIZE(orbitals, 1)
2035 2723 : iwb1 = orbitals(i3, isp2)
2036 2723 : IF (eri_env%eri_gpw%store_wfn) THEN
2037 2698 : wfn3 => wfn_a(iwb1, isp2)
2038 : ELSE
2039 : CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwb1, wfn3, rho_g, atomic_kind_set, &
2040 25 : qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
2041 : END IF
2042 2723 : CALL pw_zero(wfn_r)
2043 2723 : CALL pw_multiply(wfn_r, rho_r, wfn3)
2044 2723 : iwbt = i3
2045 2723 : IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
2046 9392 : DO i4 = iwbt, SIZE(orbitals, 1)
2047 5814 : iwb2 = orbitals(i4, isp2)
2048 5814 : IF (eri_env%eri_gpw%store_wfn) THEN
2049 5784 : wfn4 => wfn_a(iwb2, isp2)
2050 : ELSE
2051 : CALL calculate_wavefunction(fm_mo_coeff_as(isp1), iwb2, wfn4, rho_g, atomic_kind_set, &
2052 30 : 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 5814 : erint = pw_integral_ab(wfn_r, wfn4, local_only=.TRUE.)
2056 5814 : icount2 = icount2 + 1
2057 5814 : eri(icount2) = erint
2058 8537 : eri_index(icount2) = csr_idx_to_combined(iwb1, iwb2, nmo2)
2059 : END DO
2060 : END DO
2061 : ! Now, we sum the integrals globally
2062 855 : CALL eri_env%para_env_sub%sum(eri)
2063 : ! and we reorder the integrals to prevent storing too small integrals
2064 855 : intcount = 0
2065 855 : icount2 = 0
2066 : iwbs = 1
2067 : IF (isp1 == isp2) iwbs = i1
2068 855 : isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
2069 3578 : DO i3 = iwbs, SIZE(orbitals, 1)
2070 2723 : iwb1 = orbitals(i3, isp2)
2071 2723 : iwbt = i3
2072 2723 : IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
2073 9392 : DO i4 = iwbt, SIZE(orbitals, 1)
2074 5814 : iwb2 = orbitals(i4, isp2)
2075 5814 : intcount = intcount + 1
2076 5814 : erint = eri(intcount)
2077 8537 : IF (ABS(erint) > eri_env%eps_integral) THEN
2078 5138 : IF (MOD(intcount, eri_env%para_env_sub%num_pe) == eri_env%para_env_sub%mepos) THEN
2079 2571 : icount2 = icount2 + 1
2080 2571 : eri(icount2) = erint
2081 2571 : eri_index(icount2) = eri_index(intcount)
2082 : END IF
2083 : END IF
2084 : END DO
2085 : END DO
2086 855 : stored_integrals = stored_integrals + icount2
2087 : !
2088 855 : CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
2089 : !
2090 2491 : DEALLOCATE (eri, eri_index)
2091 : END DO
2092 : ELSE
2093 0 : CPABORT("Unknown option")
2094 : END IF
2095 : END DO
2096 : END DO
2097 : END DO
2098 :
2099 108 : IF (print1 .AND. iw > 0) THEN
2100 3 : WRITE (iw, "(T2,'ERI_GPW|',' Number of Integrals stored locally',T71,I10)") stored_integrals
2101 : END IF
2102 :
2103 108 : IF (eri_env%eri_gpw%store_wfn) THEN
2104 206 : DO ispin = 1, nspins
2105 548 : DO i1 = 1, SIZE(orbitals, 1)
2106 342 : iwfn = orbitals(i1, ispin)
2107 452 : CALL wfn_a(iwfn, ispin)%release()
2108 : END DO
2109 : END DO
2110 96 : DEALLOCATE (wfn_a)
2111 : ELSE
2112 12 : CALL wfn1%release()
2113 12 : CALL wfn2%release()
2114 12 : DEALLOCATE (wfn1, wfn2)
2115 12 : IF (eri_env%method /= eri_method_gpw_ht) THEN
2116 6 : CALL wfn3%release()
2117 6 : CALL wfn4%release()
2118 6 : DEALLOCATE (wfn3, wfn4)
2119 : END IF
2120 : END IF
2121 108 : CALL auxbas_pw_pool%give_back_pw(wfn_r)
2122 108 : CALL auxbas_pw_pool%give_back_pw(rho_g)
2123 108 : CALL auxbas_pw_pool%give_back_pw(rho_r)
2124 108 : CALL auxbas_pw_pool%give_back_pw(pot_g)
2125 :
2126 108 : IF (eri_env%method == eri_method_gpw_ht) THEN
2127 32 : DO ispin = 1, nspins
2128 16 : CALL dbcsr_release(mo_coeff_as(ispin))
2129 16 : CALL dbcsr_release(matrix_pq_rnu(ispin))
2130 16 : CALL cp_fm_release(fm_matrix_pq_rnu(ispin))
2131 32 : CALL cp_fm_release(fm_matrix_pq_rs(ispin))
2132 : END DO
2133 16 : DEALLOCATE (matrix_pq_rnu, fm_matrix_pq_rnu, fm_matrix_pq_rs)
2134 16 : CALL deallocate_task_list(task_list_sub)
2135 : END IF
2136 230 : DO ispin = 1, nspins
2137 122 : CALL dbcsr_release(mo_coeff_as(ispin))
2138 230 : CALL cp_fm_release(fm_mo_coeff_as(ispin))
2139 : END DO
2140 108 : DEALLOCATE (mo_coeff_as, fm_mo_coeff_as)
2141 108 : CALL release_neighbor_list_sets(sab_orb_sub)
2142 108 : CALL cp_blacs_env_release(blacs_env_sub)
2143 108 : CALL dbcsr_release(mat_munu%matrix)
2144 108 : DEALLOCATE (mat_munu%matrix)
2145 108 : CALL pw_env_release(pw_env_sub)
2146 : ! Return to the old qs_control
2147 108 : dft_control%qs_control => qs_control_old
2148 108 : DEALLOCATE (qs_control%e_cutoff)
2149 108 : DEALLOCATE (qs_control)
2150 :
2151 : ! print out progress
2152 108 : IF (iw > 0) THEN
2153 54 : t2 = m_walltime()
2154 54 : WRITE (iw, '(/,T2,A,T66,F14.2)') "ERI_GPW| ERI calculation took (sec)", t2 - t1
2155 54 : CALL m_flush(iw)
2156 : END IF
2157 :
2158 108 : CALL timestop(handle)
2159 :
2160 216 : 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 108 : 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 118 : SELECT CASE (green%method)
2184 : CASE (PERIODIC3D)
2185 :
2186 114 : SELECT CASE (eri_env%operator)
2187 : CASE (eri_operator_coulomb)
2188 786435 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2189 786429 : g2 = grid%gsq(ig)
2190 786435 : gf%array(ig) = fourpi/g2
2191 : END DO
2192 6 : IF (grid%have_g0) gf%array(1) = 0.0_dp
2193 :
2194 : CASE (eri_operator_yukawa)
2195 0 : CALL cp_warn(__LOCATION__, "Yukawa operator has not been tested")
2196 0 : omega2 = eri_env%omega**2
2197 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2198 0 : g2 = grid%gsq(ig)
2199 0 : gf%array(ig) = fourpi/(omega2 + g2)
2200 : END DO
2201 0 : IF (grid%have_g0) gf%array(1) = fourpi/omega2
2202 :
2203 : CASE (eri_operator_erf)
2204 0 : omega2 = eri_env%omega**2
2205 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2206 0 : g2 = grid%gsq(ig)
2207 0 : gf%array(ig) = fourpi/g2*EXP(-0.25_dp*g2/omega2)
2208 : END DO
2209 0 : IF (grid%have_g0) gf%array(1) = 0.0_dp
2210 :
2211 : CASE (eri_operator_erfc)
2212 0 : omega2 = eri_env%omega**2
2213 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2214 0 : g2 = grid%gsq(ig)
2215 0 : gf%array(ig) = fourpi/g2*(1.0_dp - EXP(-0.25_dp*g2/omega2))
2216 : END DO
2217 0 : IF (grid%have_g0) gf%array(1) = pi/omega2
2218 :
2219 : CASE (eri_operator_trunc)
2220 0 : Rc = eri_env%cutoff_radius
2221 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2222 0 : g2 = grid%gsq(ig)
2223 0 : g = SQRT(g2)
2224 : ! Taylor expansion around zero
2225 0 : IF (g*Rc >= 0.005_dp) THEN
2226 0 : gf%array(ig) = fourpi/g2*(1.0_dp - COS(g*Rc))
2227 : ELSE
2228 0 : 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 0 : IF (grid%have_g0) gf%array(1) = twopi*Rc**2
2232 :
2233 : CASE (eri_operator_lr_trunc)
2234 4 : omega = eri_env%omega
2235 4 : omega2 = omega**2
2236 4 : Rc = eri_env%cutoff_radius
2237 4 : Rc2 = Rc**2
2238 4 : G0 = 0.001_dp ! threshold for the G=0 case
2239 4 : Ginf = 20.0_dp ! threshold for the Taylor exapnsion arounf G=∞
2240 843752 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2241 843748 : g2 = grid%gsq(ig)
2242 843748 : g = SQRT(g2)
2243 843752 : IF (g <= 2.0_dp*G0) THEN
2244 : gf%array(ig) = -pi/omega2*erf(omega*Rc) &
2245 : + twopi*Rc2*erf(omega*Rc) &
2246 0 : + 2*rootpi*Rc*EXP(-omega2*Rc2)/omega
2247 843748 : ELSE IF (g >= 2.0_dp*Ginf*omega) THEN
2248 : ! exponential prefactor
2249 1488 : exp_prefac = EXP(-omega2*Rc2)/(rootpi*(omega2*Rc2 + 0.25_dp*g2/omega2))
2250 : ! cos sin factor
2251 1488 : cossin_fac = omega*Rc*COS(g*Rc) - 0.5_dp*g/omega*SIN(g*Rc)
2252 : ! real erf term with cosine
2253 1488 : erfcos_fac = ERF(omega*Rc)*COS(g*Rc)
2254 : ! Combine terms
2255 1488 : gf%array(ig) = fourpi/g2*(-exp_prefac*cossin_fac - erfcos_fac)
2256 : ELSE
2257 : ! exponential prefactor
2258 842260 : exp_prefac = twopi/g2*EXP(-0.25_dp*g2/omega2)
2259 : ! Compute complex arguments for erf
2260 842260 : z_p = CMPLX(omega*Rc, 0.5_dp*g/omega, kind=dp)
2261 : ! Evaluate complex error functions
2262 842260 : erf_fac_p = 2.0_dp*REAL(erfz_fast(z_p))
2263 : ! Real erf term with cosine
2264 842260 : erfcos_fac = fourpi/g2*ERF(omega*Rc)*COS(g*Rc)
2265 : ! Combine terms
2266 842260 : gf%array(ig) = exp_prefac*erf_fac_p - erfcos_fac
2267 : END IF
2268 : END DO
2269 4 : IF (grid%have_g0) THEN
2270 : gf%array(1) = -pi/omega2*ERF(omega*Rc) &
2271 : + twopi*Rc2*ERF(omega*Rc) &
2272 2 : + 2*rootpi*Rc*EXP(-omega2*Rc2)/omega
2273 : END IF
2274 :
2275 : CASE DEFAULT
2276 10 : 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 180 : 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 82 : IF (eri_env%operator == eri_operator_coulomb) THEN
2290 82 : Rc = green%radius
2291 : ELSE
2292 0 : Rc = eri_env%cutoff_radius
2293 : END IF
2294 24981734 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2295 24981652 : g2 = grid%gsq(ig)
2296 24981652 : g = SQRT(g2)
2297 : ! Taylor expansion around zero
2298 24981734 : IF (g*Rc >= 0.005_dp) THEN
2299 24981652 : gf%array(ig) = fourpi/g2*(1.0_dp - COS(g*Rc))
2300 : ELSE
2301 0 : 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 82 : IF (grid%have_g0) gf%array(1) = twopi*Rc**2
2305 :
2306 : ! Not tested
2307 : CASE (eri_operator_yukawa)
2308 0 : CALL cp_warn(__LOCATION__, "Yukawa operator has not been tested")
2309 0 : Rc = green%radius
2310 0 : omega = eri_env%omega
2311 0 : ea = EXP(-omega*Rc)
2312 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2313 0 : g2 = grid%gsq(ig)
2314 0 : g = SQRT(g2)
2315 0 : g3d = fourpi/(omega**2 + g2)
2316 0 : gf%array(ig) = g3d*(1.0_dp - ea*(COS(g*Rc) + omega/g*SIN(g*Rc)))
2317 : END DO
2318 0 : 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 16 : IF (eri_env%operator == eri_operator_erf) THEN
2324 16 : Rc = green%radius
2325 : ELSE
2326 0 : Rc = eri_env%cutoff_radius
2327 : END IF
2328 16 : omega2 = eri_env%omega**2
2329 1728008 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2330 1727992 : g2 = grid%gsq(ig)
2331 1727992 : g = SQRT(g2)
2332 1727992 : ga = -0.25_dp*g2/omega2
2333 1728008 : gf%array(ig) = fourpi/g2*EXP(ga)*(1.0_dp - COS(g*Rc))
2334 : END DO
2335 16 : 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 0 : "Short-range Coulomb operator may be incorrect with ANALYTIC0D Poisson solver")
2342 0 : Rc = green%radius
2343 0 : omega2 = eri_env%omega**2
2344 0 : DO ig = grid%first_gne0, grid%ngpts_cut_local
2345 0 : g2 = grid%gsq(ig)
2346 0 : g = SQRT(g2)
2347 0 : ga = -0.25_dp*g2/omega2
2348 0 : gf%array(ig) = fourpi/g2*(1.0_dp - EXP(ga))*(1.0_dp - COS(g*Rc))
2349 : END DO
2350 0 : IF (grid%have_g0) gf%array(1) = pi/omega2
2351 :
2352 : CASE DEFAULT
2353 98 : CPABORT("Unsupported operator")
2354 : END SELECT
2355 :
2356 : CASE DEFAULT
2357 108 : CPABORT("Unsupported Poisson solver")
2358 : END SELECT
2359 : END ASSOCIATE
2360 :
2361 108 : 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 911 : 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 911 : IF (irow /= 0) THEN
2384 911 : nze = csr_mat%nze_local
2385 911 : nze_new = nze + nnz
2386 : ! values
2387 911 : CALL reallocate(csr_mat%nzval_local%r_dp, 1, nze_new)
2388 3572 : csr_mat%nzval_local%r_dp(nze + 1:nze_new) = rdat(1:nnz)
2389 : ! col indices
2390 911 : CALL reallocate(csr_mat%colind_local, 1, nze_new)
2391 3572 : csr_mat%colind_local(nze + 1:nze_new) = rind(1:nnz)
2392 : ! rows
2393 911 : nrow = csr_mat%nrows_local
2394 911 : CALL reallocate(csr_mat%rowptr_local, 1, irow + 1)
2395 2426 : csr_mat%rowptr_local(nrow + 1:irow) = nze + 1
2396 911 : csr_mat%rowptr_local(irow + 1) = nze_new + 1
2397 : ! nzerow
2398 911 : CALL reallocate(csr_mat%nzerow_local, 1, irow)
2399 2426 : DO k = nrow + 1, irow
2400 2426 : csr_mat%nzerow_local(k) = csr_mat%rowptr_local(k + 1) - csr_mat%rowptr_local(k)
2401 : END DO
2402 911 : csr_mat%nrows_local = irow
2403 911 : csr_mat%nze_local = csr_mat%nze_local + nnz
2404 : END IF
2405 911 : csr_mat%nze_total = csr_mat%nze_total + nnz
2406 911 : csr_mat%has_indices = .TRUE.
2407 :
2408 911 : 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 4 : 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 4 : INTEGER, DIMENSION(:), POINTER :: alist, blist, istride
2424 : LOGICAL :: do_mo, explicit_a, explicit_b
2425 4 : 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 4 : 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 4 : 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 4 : CALL section_vals_val_get(input, "FILENAME", c_val=filebody)
2441 4 : CALL section_vals_val_get(input, "STRIDE", i_vals=istride)
2442 4 : IF (SIZE(istride) == 1) THEN
2443 16 : str(1:3) = istride(1)
2444 0 : ELSE IF (SIZE(istride) == 3) THEN
2445 0 : str(1:3) = istride(1:3)
2446 : ELSE
2447 0 : CPABORT("STRIDE arguments inconsistent")
2448 : END IF
2449 4 : CALL section_vals_val_get(input, "ALIST", i_vals=alist, explicit=explicit_a)
2450 4 : 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 4 : input=scf_input)
2462 :
2463 4 : CALL qs_subsys_get(subsys, particles=particles)
2464 : !
2465 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
2466 4 : CALL auxbas_pw_pool%create_pw(wf_r)
2467 4 : CALL auxbas_pw_pool%create_pw(wf_g)
2468 : !
2469 4 : dft_section => section_vals_get_subs_vals(scf_input, "DFT")
2470 : !
2471 8 : DO isp = 1, SIZE(mos)
2472 4 : CALL get_mo_set(mo_set=mos(isp), mo_coeff=mo_coeff, nmo=nmo)
2473 :
2474 4 : IF (SIZE(mos) > 1) THEN
2475 0 : SELECT CASE (isp)
2476 : CASE (1)
2477 : CALL write_mo_set_to_output_unit(mos(isp), qs_kind_set, particle_set, &
2478 0 : 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 0 : dft_section, 4, 0, final_mos=.TRUE., spin="BETA")
2482 : CASE DEFAULT
2483 0 : CPABORT("Invalid spin")
2484 : END SELECT
2485 : ELSE
2486 : CALL write_mo_set_to_output_unit(mos(isp), qs_kind_set, particle_set, &
2487 4 : dft_section, 4, 0, final_mos=.TRUE.)
2488 : END IF
2489 :
2490 44 : DO imo = 1, nmo
2491 32 : IF (isp == 1 .AND. explicit_a) THEN
2492 32 : IF (alist(1) == -1) THEN
2493 : do_mo = .TRUE.
2494 : ELSE
2495 32 : do_mo = .FALSE.
2496 128 : DO i = 1, SIZE(alist)
2497 128 : IF (imo == alist(i)) do_mo = .TRUE.
2498 : END DO
2499 : END IF
2500 0 : ELSE IF (isp == 2 .AND. explicit_b) THEN
2501 0 : IF (blist(1) == -1) THEN
2502 : do_mo = .TRUE.
2503 : ELSE
2504 0 : do_mo = .FALSE.
2505 0 : DO i = 1, SIZE(blist)
2506 0 : IF (imo == blist(i)) do_mo = .TRUE.
2507 : END DO
2508 : END IF
2509 : ELSE
2510 : do_mo = .TRUE.
2511 : END IF
2512 32 : IF (.NOT. do_mo) CYCLE
2513 : CALL calculate_wavefunction(mo_coeff, imo, wf_r, wf_g, atomic_kind_set, &
2514 12 : qs_kind_set, cell, dft_control, particle_set, pw_env)
2515 12 : IF (para_env%is_source()) THEN
2516 6 : WRITE (filename, '(A,A1,I4.4,A1,I1.1,A)') TRIM(filebody), "_", imo, "_", isp, ".cube"
2517 6 : CALL open_file(filename, unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE")
2518 6 : WRITE (title, *) "Active Orbital ", imo, " spin ", isp
2519 : ELSE
2520 6 : unit_nr = -1
2521 : END IF
2522 12 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=istride)
2523 16 : IF (para_env%is_source()) THEN
2524 26 : CALL close_file(unit_nr)
2525 : END IF
2526 : END DO
2527 : END DO
2528 :
2529 4 : CALL auxbas_pw_pool%give_back_pw(wf_r)
2530 4 : CALL auxbas_pw_pool%give_back_pw(wf_g)
2531 :
2532 4 : 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 104 : 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 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fmat
2552 : TYPE(cp_logger_type), POINTER :: logger
2553 : TYPE(eri_fcidump_checksum) :: eri_checksum
2554 :
2555 104 : checksum = 0.0_dp
2556 :
2557 208 : logger => cp_get_default_logger()
2558 : iw = cp_print_key_unit_nr(logger, as_input, "FCIDUMP", &
2559 104 : extension=".fcidump", file_status="REPLACE", file_action="WRITE", file_form="FORMATTED")
2560 : !
2561 104 : nspins = active_space_env%nspins
2562 104 : norb = SIZE(active_space_env%active_orbitals, 1)
2563 104 : ms2 = active_space_env%multiplicity - 1
2564 104 : IF (nspins == 1 .OR. restricted) THEN
2565 : ! Closed shell or restricted open-shell
2566 : ASSOCIATE (nelec => active_space_env%nelec_active)
2567 :
2568 92 : IF (iw > 0) THEN
2569 46 : WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
2570 46 : isym = 1
2571 187 : WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb)
2572 46 : isym = 0
2573 46 : WRITE (iw, "(A,I1,A)") " ISYM=", isym, ","
2574 46 : IF (restricted) WRITE (iw, "(A,I1,A)") " UHF=", 0, ","
2575 46 : 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 92 : eri_fcidump_print(iw, 1, 1), 1, 1)
2581 92 : CALL eri_checksum%set(1, 1)
2582 92 : 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 92 : nmo = active_space_env%eri%norb
2587 368 : ALLOCATE (fmat(nmo, nmo))
2588 92 : CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2589 92 : IF (iw > 0) THEN
2590 46 : i3 = 0; i4 = 0
2591 187 : DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
2592 141 : i1 = active_space_env%active_orbitals(m1, 1)
2593 526 : DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
2594 339 : i2 = active_space_env%active_orbitals(m2, 1)
2595 339 : checksum = checksum + ABS(fmat(i1, i2))
2596 480 : WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2597 : END DO
2598 : END DO
2599 : END IF
2600 92 : DEALLOCATE (fmat)
2601 : ! Print energy
2602 92 : esub = active_space_env%energy_inactive
2603 92 : i1 = 0; i2 = 0; i3 = 0; i4 = 0
2604 92 : checksum = checksum + ABS(esub)
2605 184 : 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 12 : IF (iw > 0) THEN
2612 6 : WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
2613 6 : isym = 1
2614 21 : WRITE (iw, "(A,1000(I1,','))") " ORBSYM=", (isym, i=1, norb)
2615 6 : isym = 0
2616 6 : WRITE (iw, "(A,I1,A)") " ISYM=", isym, ","
2617 6 : WRITE (iw, "(A,I1,A)") " UHF=", 1, ","
2618 6 : 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 12 : eri_fcidump_print(iw, 1, 1), 1, 1)
2625 12 : CALL eri_checksum%set(1, 1)
2626 12 : 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 12 : eri_fcidump_print(iw, 1, norb + 1), 1, 2)
2630 12 : CALL eri_checksum%set(1, norb + 1)
2631 12 : 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 12 : eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2)
2635 12 : CALL eri_checksum%set(norb + 1, norb + 1)
2636 12 : CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2)
2637 : ! Print integrals: Fij
2638 : ! alpha
2639 12 : nmo = active_space_env%eri%norb
2640 48 : ALLOCATE (fmat(nmo, nmo))
2641 12 : CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
2642 12 : IF (iw > 0) THEN
2643 6 : i3 = 0; i4 = 0
2644 21 : DO m1 = 1, norb
2645 15 : i1 = active_space_env%active_orbitals(m1, 1)
2646 48 : DO m2 = m1, norb
2647 27 : i2 = active_space_env%active_orbitals(m2, 1)
2648 27 : checksum = checksum + ABS(fmat(i1, i2))
2649 42 : WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
2650 : END DO
2651 : END DO
2652 : END IF
2653 12 : DEALLOCATE (fmat)
2654 : ! beta
2655 36 : ALLOCATE (fmat(nmo, nmo))
2656 12 : CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat)
2657 12 : IF (iw > 0) THEN
2658 6 : i3 = 0; i4 = 0
2659 21 : DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
2660 15 : i1 = active_space_env%active_orbitals(m1, 2)
2661 48 : DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
2662 27 : i2 = active_space_env%active_orbitals(m2, 2)
2663 27 : checksum = checksum + ABS(fmat(i1, i2))
2664 42 : WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4
2665 : END DO
2666 : END DO
2667 : END IF
2668 12 : DEALLOCATE (fmat)
2669 : ! Print energy
2670 12 : esub = active_space_env%energy_inactive
2671 12 : i1 = 0; i2 = 0; i3 = 0; i4 = 0
2672 12 : checksum = checksum + ABS(esub)
2673 24 : IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4
2674 : END ASSOCIATE
2675 : END IF
2676 : !
2677 104 : CALL cp_print_key_finished_output(iw, logger, as_input, "FCIDUMP")
2678 :
2679 : !>>
2680 104 : iw = cp_logger_get_default_io_unit(logger)
2681 104 : IF (iw > 0) WRITE (iw, '(T4,A,T66,F12.8)') "FCIDUMP| Checksum:", eri_checksum%checksum + checksum
2682 : !<<
2683 :
2684 208 : 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 440 : 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 10520 : replicated_matrix(:, :) = 0.0_dp
2701 2182 : DO i1 = 1, norb
2702 7222 : DO i2 = i1, norb
2703 5040 : CALL cp_fm_get_element(distributed_matrix, i1, i2, mval)
2704 5040 : replicated_matrix(i1, i2) = mval
2705 6782 : replicated_matrix(i2, i1) = mval
2706 : END DO
2707 : END DO
2708 440 : 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 130 : 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 130 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ks_a_mat, ks_a_ref, ks_b_mat, ks_b_ref, &
2725 130 : 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 130 : eref = active_space_env%energy_ref
2730 130 : nspins = active_space_env%nspins
2731 :
2732 130 : IF (nspins == 1) THEN
2733 98 : 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 784 : ALLOCATE (ks_mat(norb, norb), ks_ref(norb, norb), p_mat(norb, norb))
2739 98 : 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 98 : CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_mat)
2743 98 : 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 98 : 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 98 : 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 2542 : 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 2542 : 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 2542 : 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 98 : active_space_env%energy_inactive = esub
2765 :
2766 98 : CALL cp_fm_release(active_space_env%fock_sub)
2767 392 : ALLOCATE (active_space_env%fock_sub(nspins))
2768 196 : DO is = 1, nspins
2769 98 : matrix => active_space_env%ks_sub(is)
2770 : CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2771 196 : name="Active Fock operator")
2772 : END DO
2773 98 : matrix => active_space_env%fock_sub(1)
2774 594 : DO i1 = 1, norb
2775 2542 : DO i2 = 1, norb
2776 2046 : mval = ks_mat(i1, i2)
2777 2444 : CALL cp_fm_set_element(matrix, i1, i2, mval)
2778 : END DO
2779 : END DO
2780 : ELSE
2781 :
2782 32 : 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 448 : & p_a_mat(norb, norb), p_b_mat(norb, norb))
2790 32 : ks_a_ref(:, :) = 0.0_dp; ks_b_ref(:, :) = 0.0_dp
2791 :
2792 32 : CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_a_mat)
2793 32 : CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(2), p_b_mat)
2794 32 : CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_a_mat)
2795 32 : CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(2), ks_b_mat)
2796 : !
2797 : !
2798 32 : IF (restricted) THEN
2799 : ! In the restricted case, we use the same ERIs for each spin channel
2800 12 : 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 12 : 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 12 : restricted=restricted)
2807 : ELSE
2808 20 : eri_aa => active_space_env%eri%eri(1)%csr_mat
2809 20 : eri_ab => active_space_env%eri%eri(2)%csr_mat
2810 20 : 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 20 : 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 20 : restricted=restricted)
2817 : END IF
2818 : !
2819 : ! calculate energy
2820 32 : eeri = 0.0_dp
2821 1336 : eeri = 0.5_dp*(SUM(ks_a_ref*p_a_mat) + SUM(ks_b_ref*p_b_mat))
2822 1336 : esub = eref - SUM(ks_a_mat*p_a_mat) - SUM(ks_b_mat*p_b_mat) + eeri
2823 668 : ks_a_mat(:, :) = ks_a_mat(:, :) - ks_a_ref(:, :)
2824 668 : ks_b_mat(:, :) = ks_b_mat(:, :) - ks_b_ref(:, :)
2825 : !
2826 32 : active_space_env%energy_inactive = esub
2827 : !
2828 32 : CALL cp_fm_release(active_space_env%fock_sub)
2829 160 : ALLOCATE (active_space_env%fock_sub(nspins))
2830 96 : DO is = 1, nspins
2831 64 : matrix => active_space_env%ks_sub(is)
2832 : CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
2833 96 : name="Active Fock operator")
2834 : END DO
2835 :
2836 32 : matrix => active_space_env%fock_sub(1)
2837 156 : DO i1 = 1, norb
2838 668 : DO i2 = 1, norb
2839 512 : mval = ks_a_mat(i1, i2)
2840 636 : CALL cp_fm_set_element(matrix, i1, i2, mval)
2841 : END DO
2842 : END DO
2843 32 : matrix => active_space_env%fock_sub(2)
2844 188 : DO i1 = 1, norb
2845 668 : DO i2 = 1, norb
2846 512 : mval = ks_b_mat(i1, i2)
2847 636 : CALL cp_fm_set_element(matrix, i1, i2, mval)
2848 : END DO
2849 : END DO
2850 :
2851 : END IF
2852 :
2853 130 : 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 98 : 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 98 : CALL timeset(routineN, handle)
2879 :
2880 : ! Nothing to do
2881 98 : norb = SIZE(active_orbitals, 1)
2882 98 : nmo_total = SIZE(p_mat, 1)
2883 98 : nindex = (nmo_total*(nmo_total + 1))/2
2884 98 : CALL mp_group%set_handle(eri%mp_group%get_handle())
2885 424 : DO m1 = 1, norb
2886 326 : i1 = active_orbitals(m1, 1)
2887 1264 : DO m2 = m1, norb
2888 840 : i2 = active_orbitals(m2, 1)
2889 840 : i12 = csr_idx_to_combined(i1, i2, nmo_total)
2890 1166 : IF (MOD(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
2891 831 : i12l = (i12 - 1)/comm_exchange%num_pe + 1
2892 831 : irptr = eri%rowptr_local(i12l) - 1
2893 3486 : DO i34l = 1, eri%nzerow_local(i12l)
2894 2655 : i34 = eri%colind_local(irptr + i34l)
2895 2655 : CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
2896 2655 : erint = eri%nzval_local%r_dp(irptr + i34l)
2897 : ! Coulomb
2898 2655 : ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2899 2655 : IF (i3 /= i4) THEN
2900 1614 : ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
2901 : END IF
2902 2655 : IF (i12 /= i34) THEN
2903 2235 : ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2904 2235 : IF (i1 /= i2) THEN
2905 1617 : ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
2906 : END IF
2907 : END IF
2908 : ! Exchange
2909 2655 : erint = -0.5_dp*erint
2910 2655 : ks_ref(i1, i3) = ks_ref(i1, i3) + erint*p_mat(i2, i4)
2911 2655 : IF (i1 /= i2) THEN
2912 1874 : ks_ref(i2, i3) = ks_ref(i2, i3) + erint*p_mat(i1, i4)
2913 : END IF
2914 2655 : IF (i3 /= i4) THEN
2915 1614 : ks_ref(i1, i4) = ks_ref(i1, i4) + erint*p_mat(i2, i3)
2916 : END IF
2917 2655 : IF (i1 /= i2 .AND. i3 /= i4) THEN
2918 1253 : 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 6141 : IF (i12 /= i34) THEN
2922 2235 : ks_ref(i3, i1) = ks_ref(i3, i1) + erint*p_mat(i4, i2)
2923 2235 : IF (i1 /= i2) THEN
2924 1617 : ks_ref(i3, i2) = ks_ref(i3, i2) + erint*p_mat(i4, i1)
2925 : END IF
2926 2235 : IF (i3 /= i4) THEN
2927 1357 : ks_ref(i4, i1) = ks_ref(i4, i1) + erint*p_mat(i3, i2)
2928 : END IF
2929 2235 : IF (i1 /= i2 .AND. i3 /= i4) THEN
2930 996 : 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 424 : DO m1 = 1, norb
2939 326 : i1 = active_orbitals(m1, 1)
2940 1264 : DO m2 = m1, norb
2941 840 : i2 = active_orbitals(m2, 1)
2942 1166 : ks_ref(i2, i1) = ks_ref(i1, i2)
2943 : END DO
2944 : END DO
2945 4986 : CALL mp_group%sum(ks_ref)
2946 :
2947 98 : CALL timestop(handle)
2948 :
2949 98 : 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 64 : 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 64 : CALL timeset(routineN, handle)
2982 :
2983 64 : norb = SIZE(active_orbitals, 1)
2984 64 : nmo_total = SIZE(p_a_mat, 1)
2985 64 : nindex = (nmo_total*(nmo_total + 1))/2
2986 64 : IF (tr_mixed_eri) THEN
2987 : spin1 = 2
2988 64 : spin2 = 1
2989 : ELSE
2990 32 : spin1 = 1
2991 32 : spin2 = 2
2992 : END IF
2993 276 : DO m1 = 1, norb
2994 212 : i1 = active_orbitals(m1, spin1)
2995 756 : DO m2 = m1, norb
2996 480 : i2 = active_orbitals(m2, spin1)
2997 480 : i12 = csr_idx_to_combined(i1, i2, nmo_total)
2998 692 : IF (MOD(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
2999 480 : i12l = (i12 - 1)/comm_exchange%num_pe + 1
3000 480 : irptr = eri_aa%rowptr_local(i12l) - 1
3001 1631 : DO i34l = 1, eri_aa%nzerow_local(i12l)
3002 1151 : i34 = eri_aa%colind_local(irptr + i34l)
3003 1151 : CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
3004 1151 : erint = eri_aa%nzval_local%r_dp(irptr + i34l)
3005 : ! Coulomb
3006 1151 : p12 = p_a_mat(i1, i2)
3007 1151 : p34 = p_a_mat(i3, i4)
3008 1151 : IF (restricted) THEN
3009 560 : p12 = p12 + p_b_mat(i1, i2)
3010 560 : p34 = p34 + p_b_mat(i3, i4)
3011 : END IF
3012 1151 : IF (i1 /= i2) p12 = 2.0_dp*p12
3013 1151 : IF (i3 /= i4) p34 = 2.0_dp*p34
3014 : !F_ij += (ij|kl)*d_kl
3015 1151 : ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p34
3016 1151 : IF (i12 /= i34) THEN
3017 : !F_kl += (ij|kl)*d_ij
3018 911 : ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p12
3019 : END IF
3020 : ! Exchange
3021 1151 : erint = -1.0_dp*erint
3022 : !F_ik -= (ij|kl)*d_jl
3023 1151 : ks_a_ref(i1, i3) = ks_a_ref(i1, i3) + erint*p_a_mat(i2, i4)
3024 1151 : IF (i1 /= i2) THEN
3025 : !F_jk -= (ij|kl)*d_il
3026 705 : ks_a_ref(i2, i3) = ks_a_ref(i2, i3) + erint*p_a_mat(i1, i4)
3027 : END IF
3028 1151 : IF (i3 /= i4) THEN
3029 : !F_il -= (ij|kl)*d_jk
3030 628 : ks_a_ref(i1, i4) = ks_a_ref(i1, i4) + erint*p_a_mat(i2, i3)
3031 : END IF
3032 1151 : IF (i1 /= i2 .AND. i3 /= i4) THEN
3033 : !F_jl -= (ij|kl)*d_ik
3034 422 : ks_a_ref(i2, i4) = ks_a_ref(i2, i4) + erint*p_a_mat(i1, i3)
3035 : END IF
3036 2782 : IF (i12 /= i34) THEN
3037 911 : ks_a_ref(i3, i1) = ks_a_ref(i3, i1) + erint*p_a_mat(i4, i2)
3038 911 : IF (i1 /= i2) THEN
3039 571 : ks_a_ref(i3, i2) = ks_a_ref(i3, i2) + erint*p_a_mat(i4, i1)
3040 : END IF
3041 911 : IF (i3 /= i4) THEN
3042 494 : ks_a_ref(i4, i1) = ks_a_ref(i4, i1) + erint*p_a_mat(i3, i2)
3043 : END IF
3044 911 : IF (i1 /= i2 .AND. i3 /= i4) THEN
3045 288 : 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 64 : IF (.NOT. restricted) THEN
3056 164 : DO m1 = 1, norb
3057 124 : i1 = active_orbitals(m1, 1)
3058 432 : DO m2 = m1, norb
3059 268 : i2 = active_orbitals(m2, 1)
3060 268 : i12 = csr_idx_to_combined(i1, i2, nmo_total)
3061 392 : IF (MOD(i12 - 1, comm_exchange%num_pe) == comm_exchange%mepos) THEN
3062 268 : i12l = (i12 - 1)/comm_exchange%num_pe + 1
3063 268 : irptr = eri_ab%rowptr_local(i12l) - 1
3064 1310 : DO i34l = 1, eri_ab%nzerow_local(i12l)
3065 1042 : i34 = eri_ab%colind_local(irptr + i34l)
3066 1042 : CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
3067 1042 : erint = eri_ab%nzval_local%r_dp(irptr + i34l)
3068 : ! Coulomb
3069 1310 : IF (tr_mixed_eri) THEN
3070 : !F_kl += (kl beta|ij alpha )*d_alpha_ij
3071 521 : ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_b_mat(i1, i2)
3072 521 : 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 521 : ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_b_mat(i3, i4)
3076 521 : 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 276 : DO m1 = 1, norb
3085 212 : i1 = active_orbitals(m1, spin1)
3086 756 : DO m2 = m1, norb
3087 480 : i2 = active_orbitals(m2, spin1)
3088 692 : ks_a_ref(i2, i1) = ks_a_ref(i1, i2)
3089 : END DO
3090 : END DO
3091 64 : CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
3092 2608 : CALL mp_group%sum(ks_a_ref)
3093 :
3094 64 : CALL timestop(handle)
3095 :
3096 64 : 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 0 : 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 0 : CHARACTER(len=6), DIMENSION(:), POINTER :: sym
3115 : INTEGER :: i, l, nj
3116 : INTEGER, DIMENSION(4, 7) :: ne
3117 0 : INTEGER, DIMENSION(:), POINTER :: lq, nq
3118 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: zet
3119 : TYPE(sto_basis_set_type), POINTER :: sto_basis_set
3120 :
3121 0 : CPASSERT(.NOT. ASSOCIATED(pro_basis_set))
3122 0 : NULLIFY (sto_basis_set)
3123 :
3124 : ! electronic configuration
3125 0 : ne = 0
3126 0 : DO l = 1, 4 !lq(1)+1
3127 0 : nj = 2*(l - 1) + 1
3128 0 : DO i = l, 7 ! nq(1)
3129 0 : ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
3130 0 : ne(l, i) = MAX(ne(l, i), 0)
3131 0 : ne(l, i) = MIN(ne(l, i), 2*nj)
3132 : END DO
3133 : END DO
3134 0 : ALLOCATE (nq(ishell), lq(ishell), zet(ishell), sym(ishell))
3135 0 : DO i = 1, ishell
3136 0 : nq(i) = nshell(i)
3137 0 : SELECT CASE (lnam(i))
3138 : CASE ('S', 's')
3139 0 : lq(i) = 0
3140 : CASE ('P', 'p')
3141 0 : lq(i) = 1
3142 : CASE ('D', 'd')
3143 0 : lq(i) = 2
3144 : CASE ('F', 'f')
3145 0 : lq(i) = 3
3146 : CASE DEFAULT
3147 0 : CPABORT("Wrong l QN")
3148 : END SELECT
3149 0 : sym(i) = lnam(i)
3150 0 : zet(i) = srules(zval, ne, nq(1), lq(1))
3151 : END DO
3152 0 : CALL allocate_sto_basis_set(sto_basis_set)
3153 0 : CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zet, symbol=sym)
3154 0 : CALL create_gto_from_sto_basis(sto_basis_set, pro_basis_set, 6)
3155 0 : pro_basis_set%norm_type = 2
3156 0 : CALL init_orb_basis_set(pro_basis_set)
3157 0 : CALL deallocate_sto_basis_set(sto_basis_set)
3158 :
3159 0 : 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 22 : 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 22 : 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 22 : nspins = active_space_env%nspins
3179 22 : mos_active => active_space_env%mos_active
3180 58 : DO ispin = 1, nspins
3181 : ! size of p_inactive_ao is (nao x nao)
3182 36 : p_inactive_ao => active_space_env%pmat_inactive(ispin)%matrix
3183 :
3184 : ! copy p_inactive_ao to rho_ao
3185 36 : CALL dbcsr_copy(rho_ao(ispin)%matrix, p_inactive_ao)
3186 :
3187 : ! size of p_active_mo is (nmo x nmo)
3188 36 : p_active_mo => active_space_env%p_active(ispin)
3189 :
3190 : ! calculate R = p_mo
3191 36 : CALL cp_fm_create(R, p_active_mo%matrix_struct)
3192 36 : CALL cp_fm_to_fm(p_active_mo, R)
3193 :
3194 : ! calculate U = C * p_mo
3195 36 : CALL get_mo_set(mos_active(ispin), mo_coeff=C_active, nao=nao, nmo=nmo)
3196 36 : CALL cp_fm_create(U, C_active%matrix_struct)
3197 36 : 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 36 : matrix_v=U, matrix_g=C_active, ncol=nmo, alpha=1.0_dp)
3201 :
3202 36 : CALL cp_fm_release(R)
3203 130 : CALL cp_fm_release(U)
3204 : END DO
3205 :
3206 22 : 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 4792 : 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 4792 : IF (this%unit_nr > 0) THEN
3225 2396 : WRITE (this%unit_nr, "(ES23.16,4I4)") val, i + this%bra_start - 1, j + this%bra_start - 1, &
3226 4792 : & k + this%ket_start - 1, l + this%ket_start - 1
3227 : END IF
3228 :
3229 4792 : cont = .TRUE.
3230 4792 : 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 4792 : 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 4792 : this%checksum = this%checksum + ABS(val)
3252 :
3253 4792 : cont = .TRUE.
3254 4792 : 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 6 : 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 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: noon, pmat
3269 : TYPE(cp_fm_type), POINTER :: p_active
3270 :
3271 6 : nspins = active_space_env%nspins
3272 6 : nmo_active = active_space_env%nmo_active
3273 :
3274 24 : ALLOCATE (noon(nmo_active, nspins))
3275 24 : ALLOCATE (pmat(nmo_active, nmo_active))
3276 :
3277 16 : DO ispin = 1, nspins
3278 10 : p_active => active_space_env%p_active(ispin)
3279 54 : noon(:, ispin) = 0.0_dp
3280 10 : pmat = 0.0_dp
3281 :
3282 54 : DO i1 = 1, nmo_active
3283 44 : m1 = active_space_env%active_orbitals(i1, ispin)
3284 254 : DO i2 = 1, nmo_active
3285 200 : m2 = active_space_env%active_orbitals(i2, ispin)
3286 244 : CALL cp_fm_get_element(p_active, m1, m2, pmat(i1, i2))
3287 : END DO
3288 : END DO
3289 :
3290 10 : IF (iw > 0) THEN
3291 5 : WRITE (iw, '(/,T3,A,I2,A)') "Active space density matrix for spin ", ispin
3292 27 : DO i1 = 1, nmo_active
3293 49 : DO ii = 1, nmo_active, 8
3294 22 : jm = MIN(7, nmo_active - ii)
3295 144 : 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 10 : CALL diamat_all(pmat, noon(:, ispin))
3302 :
3303 16 : IF (iw > 0) THEN
3304 5 : WRITE (iw, '(/,T3,A,I2,A)') "Natural orbitals occupation numbers for spin ", ispin
3305 10 : DO i1 = 1, nmo_active, 8
3306 5 : jm = MIN(7, nmo_active - i1)
3307 : ! noons are stored in ascending order, so reverse-print them
3308 32 : 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 6 : DEALLOCATE (noon)
3315 6 : DEALLOCATE (pmat)
3316 :
3317 6 : 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 6 : 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 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p_act_mo_a, p_act_mo_b
3338 : TYPE(cp_logger_type), POINTER :: logger
3339 6 : 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 6 : CALL timeset(routineN, handle)
3348 :
3349 6 : t1 = m_walltime()
3350 6 : logger => cp_get_default_logger()
3351 6 : iw = cp_logger_get_default_io_unit(logger)
3352 :
3353 6 : CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
3354 :
3355 6 : CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
3356 6 : active_space_env%do_scf_embedding = do_scf_embedding
3357 6 : CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter)
3358 6 : IF (max_iter < 0) CPABORT("Specify a non-negative number of max iterations.")
3359 6 : CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
3360 6 : IF (eps_iter < 0.0_dp) CPABORT("Specify a non-negative convergence threshold.")
3361 :
3362 6 : CALL get_qs_env(qs_env, rho=rho, energy=energy, ks_env=ks_env)
3363 6 : CALL qs_rho_get(rho, rho_ao=rho_ao)
3364 :
3365 6 : IF (iw > 0) THEN
3366 : WRITE (UNIT=iw, FMT="(/,T2,A,/)") &
3367 3 : "RANGE-SEPARATED DFT EMBEDDING WITH LIBFCI SOLVER"
3368 3 : WRITE (iw, '(T3,A,T68,I12)') "Max. iterations", max_iter
3369 3 : WRITE (iw, '(T3,A,T68,E12.4)') "Conv. threshold", eps_iter
3370 3 : WRITE (iw, '(T3,A,T68,A)') "Density mixer", TRIM(active_space_mixing_label(active_space_env))
3371 3 : 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 3 : "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", REPEAT("-", 78)
3374 : END IF
3375 :
3376 6 : iter = 0
3377 6 : converged = .FALSE.
3378 6 : energy_scf = active_space_env%energy_ref
3379 6 : energy_new = energy_scf
3380 6 : mos_active => active_space_env%mos_active
3381 :
3382 22 : DO WHILE (iter < max_iter)
3383 22 : iter = iter + 1
3384 :
3385 22 : IF (active_space_env%nspins == 2) THEN
3386 14 : CALL solve_active_space_fci(active_space_env, para_env, p_act_mo_a, p_act_mo_b)
3387 14 : active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3388 14 : CALL update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
3389 14 : DEALLOCATE (p_act_mo_a, p_act_mo_b)
3390 : ELSE
3391 8 : CALL solve_active_space_fci(active_space_env, para_env, p_act_mo_a)
3392 8 : active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
3393 8 : CALL update_active_density(p_act_mo_a, active_space_env)
3394 8 : DEALLOCATE (p_act_mo_a)
3395 : END IF
3396 :
3397 22 : energy_old = energy_new
3398 22 : energy_new = active_space_env%energy_total
3399 22 : energy_corr = energy_new - energy_scf
3400 22 : delta_E = energy_new - energy_old
3401 :
3402 22 : t2 = t1
3403 22 : t1 = m_walltime()
3404 22 : 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 11 : iter, TRIM(active_space_mixing_label(active_space_env)), &
3408 22 : t1 - t2, energy_corr, energy_new, delta_E
3409 11 : CALL m_flush(iw)
3410 : END IF
3411 :
3412 22 : CALL update_density_ao(active_space_env, rho_ao)
3413 22 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
3414 22 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
3415 22 : 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 22 : ext_xc_section=active_space_env%xc_section)
3419 :
3420 22 : active_space_env%energy_ref = energy%total
3421 22 : CALL calculate_operators(mos_active, qs_env, active_space_env)
3422 22 : CALL subspace_fock_matrix(active_space_env, dft_control%roks)
3423 :
3424 22 : IF (.NOT. active_space_env%do_scf_embedding) THEN
3425 0 : IF (iw > 0) THEN
3426 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3427 0 : "*** one-shot embedding correction finished ***"
3428 : END IF
3429 : converged = .TRUE.
3430 : EXIT
3431 22 : ELSE IF (ABS(delta_E) <= eps_iter) THEN
3432 6 : IF (iw > 0) THEN
3433 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3434 3 : "*** 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 0 : IF (iw > 0) THEN
3443 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3444 0 : "*** rs-DFT embedding did not converged after ", iter, " iteration(s) ***"
3445 : END IF
3446 : END IF
3447 :
3448 6 : energy%total = active_space_env%energy_total
3449 :
3450 6 : IF (iw > 0) THEN
3451 : WRITE (UNIT=iw, FMT="(/,T3,A)") &
3452 3 : "Final energy contributions:"
3453 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3454 3 : "Inactive energy:", active_space_env%energy_inactive
3455 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3456 3 : "Active energy:", active_space_env%energy_active
3457 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3458 3 : "Correlation energy:", energy_corr
3459 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3460 3 : "Total rs-DFT energy:", active_space_env%energy_total
3461 : END IF
3462 :
3463 6 : CALL print_pmat_noon(active_space_env, iw)
3464 6 : CALL para_env%sync()
3465 6 : CALL timestop(handle)
3466 :
3467 12 : END SUBROUTINE local_fci_embedding
3468 :
3469 : ! **************************************************************************************************
3470 : !> \brief ...
3471 : !> \param qs_env ...
3472 : !> \param active_space_env ...
3473 : !> \param as_input ...
3474 : ! **************************************************************************************************
3475 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
3505 :
3506 0 : t1 = m_walltime()
3507 :
3508 0 : logger => cp_get_default_logger()
3509 0 : iw = cp_logger_get_default_io_unit(logger)
3510 :
3511 0 : CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
3512 0 : ionode = para_env%is_source()
3513 :
3514 : ! get info from the input
3515 0 : CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
3516 0 : active_space_env%do_scf_embedding = do_scf_embedding
3517 0 : CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter)
3518 0 : IF (max_iter < 0) CPABORT("Specify a non-negative number of max iterations.")
3519 0 : CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
3520 0 : 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 0 : CALL initialize_socket(socket_fd, client_fd, as_input, ionode)
3524 0 : CALL para_env%sync()
3525 :
3526 : ! send two-electron integrals to the client
3527 0 : CALL send_eri_to_client(client_fd, active_space_env, para_env)
3528 :
3529 : ! get pointer to density in ao basis
3530 0 : CALL get_qs_env(qs_env, rho=rho, energy=energy, ks_env=ks_env)
3531 0 : CALL qs_rho_get(rho, rho_ao=rho_ao)
3532 :
3533 0 : IF (iw > 0) THEN
3534 : WRITE (UNIT=iw, FMT="(/,T2,A,/)") &
3535 0 : "RANGE-SEPARATED DFT EMBEDDING SELF-CONSISTENT OPTIMIZATION"
3536 :
3537 0 : WRITE (iw, '(T3,A,T68,I12)') "Max. iterations", max_iter
3538 0 : WRITE (iw, '(T3,A,T68,E12.4)') "Conv. threshold", eps_iter
3539 0 : WRITE (iw, '(T3,A,T68,A)') "Density mixer", TRIM(active_space_mixing_label(active_space_env))
3540 0 : 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 0 : "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", REPEAT("-", 78)
3544 : END IF
3545 : ! TODO: consider using an iteration level
3546 :
3547 0 : iter = 0
3548 0 : converged = .FALSE.
3549 : ! store the scf energy
3550 0 : energy_scf = active_space_env%energy_ref
3551 0 : energy_new = energy_scf
3552 0 : 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 0 : DO WHILE (iter < max_iter)
3557 0 : 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 0 : CALL send_fock_to_client(client_fd, active_space_env, para_env)
3562 :
3563 : ! update energies
3564 0 : energy_old = energy_new
3565 0 : energy_new = active_space_env%energy_total
3566 0 : energy_corr = energy_new - energy_scf
3567 0 : delta_E = energy_new - energy_old
3568 :
3569 : ! get timer
3570 0 : t2 = t1
3571 0 : t1 = m_walltime()
3572 : ! print out progress
3573 0 : 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 0 : iter, TRIM(active_space_mixing_label(active_space_env)), &
3577 0 : t1 - t2, energy_corr, energy_new, delta_E
3578 0 : CALL m_flush(iw)
3579 : END IF
3580 :
3581 : ! update total density in AO basis with the AS contribution
3582 0 : 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 0 : CALL qs_rho_update_rho(rho, qs_env=qs_env) ! updates rho_r and rho_g using rho_ao
3586 0 : 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 0 : 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 0 : ext_xc_section=active_space_env%xc_section)
3594 :
3595 : ! update the reference energy
3596 0 : active_space_env%energy_ref = energy%total
3597 :
3598 : ! transform KS/Fock, Vxc and Hcore from AO to MO basis
3599 0 : CALL calculate_operators(mos_active, qs_env, active_space_env)
3600 :
3601 : ! calculate the new inactive energy and embedding potential
3602 0 : CALL subspace_fock_matrix(active_space_env, dft_control%roks)
3603 :
3604 : ! check if it is a one-shot correction
3605 0 : IF (.NOT. active_space_env%do_scf_embedding) THEN
3606 0 : IF (iw > 0) THEN
3607 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3608 0 : "*** one-shot embedding correction finished ***"
3609 : END IF
3610 : converged = .TRUE.
3611 : EXIT
3612 : ! check for convergence
3613 0 : ELSE IF (ABS(delta_E) <= eps_iter) THEN
3614 0 : IF (iw > 0) THEN
3615 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3616 0 : "*** 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 0 : IF (iw > 0) THEN
3625 : WRITE (UNIT=iw, FMT="(/,T3,A,I5,A)") &
3626 0 : "*** 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 0 : energy%total = active_space_env%energy_total
3632 :
3633 : ! print final energy contributions
3634 0 : IF (iw > 0) THEN
3635 : WRITE (UNIT=iw, FMT="(/,T3,A)") &
3636 0 : "Final energy contributions:"
3637 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3638 0 : "Inactive energy:", active_space_env%energy_inactive
3639 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3640 0 : "Active energy:", active_space_env%energy_active
3641 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3642 0 : "Correlation energy:", energy_corr
3643 : WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
3644 0 : "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 0 : CALL print_pmat_noon(active_space_env, iw)
3649 :
3650 0 : CALL finalize_socket(socket_fd, client_fd, as_input, ionode)
3651 0 : CALL para_env%sync()
3652 : #endif
3653 :
3654 0 : CALL timestop(handle)
3655 :
3656 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
3680 :
3681 0 : logger => cp_get_default_logger()
3682 0 : iw = cp_logger_get_default_io_unit(logger)
3683 :
3684 : ! protocol == 0 for UNIX, protocol > 0 for INET
3685 0 : CALL section_vals_val_get(as_input, "SOCKET%INET", l_val=inet)
3686 0 : IF (inet) THEN
3687 0 : protocol = 1
3688 : ELSE
3689 0 : protocol = 0
3690 : END IF
3691 0 : CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
3692 0 : CALL section_vals_val_get(as_input, "SOCKET%PORT", i_val=port)
3693 :
3694 0 : IF (ionode) THEN
3695 0 : CALL open_bind_socket(socket_fd, protocol, port, TRIM(hostname)//C_NULL_CHAR)
3696 0 : WRITE (iw, '(/,T2,A,A)') "@SERVER: Created socket with address ", TRIM(hostname)
3697 0 : CALL listen_socket(socket_fd, backlog)
3698 :
3699 : ! wait until a connetion request arrives
3700 0 : WRITE (iw, '(T2,A)') "@SERVER: Waiting for requests..."
3701 0 : CALL accept_socket(socket_fd, client_fd)
3702 0 : WRITE (iw, '(T2,A,I2)') "@SERVER: Accepted socket with fd ", client_fd
3703 : END IF
3704 :
3705 0 : CALL timestop(handle)
3706 :
3707 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
3728 :
3729 0 : CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
3730 :
3731 0 : IF (ionode) THEN
3732 : ! signal the client to quit
3733 0 : CALL writebuffer(client_fd, "QUIT ", header_len)
3734 : ! close the connection
3735 0 : CALL close_socket(client_fd)
3736 0 : CALL close_socket(socket_fd)
3737 :
3738 : ! delete the socket file
3739 0 : IF (file_exists(TRIM(hostname))) THEN
3740 0 : CALL remove_socket_file(TRIM(hostname)//C_NULL_CHAR)
3741 : END IF
3742 : END IF
3743 :
3744 0 : CALL timestop(handle)
3745 :
3746 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb, s_ab
3766 : TYPE(cp_logger_type), POINTER :: logger
3767 :
3768 0 : CALL timeset(routineN, handle)
3769 :
3770 0 : logger => cp_get_default_logger()
3771 0 : iw = cp_logger_get_default_io_unit(logger)
3772 0 : ionode = para_env%is_source()
3773 0 : restricted_orbitals = active_space_env%restricted_orbitals
3774 :
3775 0 : ALLOCATE (eri_aa(active_space_env%nmo_active**4))
3776 0 : CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
3777 0 : IF (active_space_env%nspins == 2) THEN
3778 0 : ALLOCATE (eri_ab(active_space_env%nmo_active**4))
3779 0 : ALLOCATE (eri_bb(active_space_env%nmo_active**4))
3780 0 : IF (restricted_orbitals) THEN
3781 0 : eri_ab(:) = eri_aa
3782 0 : eri_bb(:) = eri_aa
3783 : ELSE
3784 0 : CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
3785 0 : 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 0 : 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 0 : 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 0 : IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len)
3797 : DO
3798 0 : header = ""
3799 0 : CALL para_env%sync()
3800 0 : IF (ionode) THEN
3801 : ! IF (iw > 0) WRITE(iw, *) "@SERVER: Waiting for messages..."
3802 0 : CALL readbuffer(client_fd, header, header_len)
3803 : END IF
3804 0 : CALL para_env%bcast(header, para_env%source)
3805 :
3806 : ! IF (iw > 0) WRITE(iw, *) "@SERVER: Message from client: ", TRIM(header)
3807 :
3808 0 : IF (TRIM(header) == "READY") THEN
3809 : ! if the client is ready, send the data
3810 0 : CALL para_env%sync()
3811 0 : IF (ionode) THEN
3812 0 : CALL writebuffer(client_fd, "TWOBODY ", header_len)
3813 0 : CALL writebuffer(client_fd, active_space_env%nspins)
3814 0 : CALL writebuffer(client_fd, active_space_env%nmo_active)
3815 0 : CALL writebuffer(client_fd, active_space_env%nelec_active)
3816 0 : CALL writebuffer(client_fd, active_space_env%multiplicity)
3817 : ! send the alpha component
3818 0 : CALL writebuffer(client_fd, eri_aa, SIZE(eri_aa))
3819 : ! send the beta part for unrestricted calculations
3820 0 : IF (active_space_env%nspins == 2) THEN
3821 0 : CALL writebuffer(client_fd, eri_ab, SIZE(eri_ab))
3822 0 : CALL writebuffer(client_fd, eri_bb, SIZE(eri_bb))
3823 0 : CALL writebuffer(client_fd, s_ab, SIZE(s_ab))
3824 : END IF
3825 : END IF
3826 0 : ELSE IF (TRIM(header) == "RECEIVED") THEN
3827 : EXIT
3828 : END IF
3829 : END DO
3830 :
3831 0 : DEALLOCATE (eri_aa)
3832 0 : IF (active_space_env%nspins == 2) THEN
3833 0 : DEALLOCATE (eri_ab)
3834 0 : DEALLOCATE (eri_bb)
3835 0 : DEALLOCATE (s_ab)
3836 : END IF
3837 :
3838 0 : CALL para_env%sync()
3839 :
3840 0 : CALL timestop(handle)
3841 :
3842 0 : 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 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
3865 :
3866 : ! Set to .TRUE. to activate debug output
3867 0 : debug = .FALSE.
3868 :
3869 0 : logger => cp_get_default_logger()
3870 0 : iw = cp_logger_get_default_io_unit(logger)
3871 0 : ionode = para_env%is_source()
3872 :
3873 0 : ALLOCATE (p_act_mo_a(active_space_env%nmo_active**2))
3874 0 : ALLOCATE (fock_a(active_space_env%nmo_active**2))
3875 0 : IF (active_space_env%nspins == 2) THEN
3876 0 : ALLOCATE (p_act_mo_b(active_space_env%nmo_active**2))
3877 0 : 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 0 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices)
3883 : END ASSOCIATE
3884 :
3885 0 : IF (active_space_env%nspins == 2) THEN
3886 : ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 2))
3887 0 : 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 0 : IF (ionode) CALL writebuffer(client_fd, "STATUS ", header_len)
3893 : DO
3894 0 : header = ""
3895 :
3896 0 : CALL para_env%sync()
3897 0 : IF (ionode) THEN
3898 : IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Waiting for messages..."
3899 0 : CALL readbuffer(client_fd, header, header_len)
3900 : END IF
3901 0 : 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 0 : IF (TRIM(header) == "READY") THEN
3906 : ! if the client is ready, send the data
3907 0 : CALL para_env%sync()
3908 0 : IF (ionode) THEN
3909 0 : CALL writebuffer(client_fd, "ONEBODY ", header_len)
3910 0 : CALL writebuffer(client_fd, active_space_env%energy_inactive)
3911 : ! send the alpha component
3912 0 : CALL writebuffer(client_fd, fock_a, SIZE(fock_a))
3913 : ! send the beta part for unrestricted calculations
3914 0 : IF (active_space_env%nspins == 2) THEN
3915 0 : CALL writebuffer(client_fd, fock_b, SIZE(fock_b))
3916 : END IF
3917 : END IF
3918 :
3919 0 : ELSE IF (TRIM(header) == "HAVEDATA") THEN
3920 : ! qiskit has data to transfer, let them know we want it and wait for it
3921 0 : CALL para_env%sync()
3922 0 : IF (ionode) THEN
3923 : IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Qiskit has data to transfer"
3924 0 : CALL writebuffer(client_fd, "GETDENSITY ", header_len)
3925 :
3926 : ! read the active energy and density
3927 0 : CALL readbuffer(client_fd, active_space_env%energy_active)
3928 0 : CALL readbuffer(client_fd, p_act_mo_a, SIZE(p_act_mo_a))
3929 0 : IF (active_space_env%nspins == 2) THEN
3930 0 : 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 0 : CALL para_env%bcast(active_space_env%energy_active, para_env%source)
3936 0 : CALL para_env%bcast(p_act_mo_a, para_env%source)
3937 0 : IF (active_space_env%nspins == 2) THEN
3938 0 : 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 0 : 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 0 : IF (active_space_env%nspins == 2) THEN
3946 0 : CALL update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
3947 : ELSE
3948 0 : 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 0 : DEALLOCATE (p_act_mo_a)
3958 0 : DEALLOCATE (fock_a)
3959 0 : IF (active_space_env%nspins == 2) THEN
3960 0 : DEALLOCATE (p_act_mo_b)
3961 0 : DEALLOCATE (fock_b)
3962 : END IF
3963 :
3964 0 : CALL para_env%sync()
3965 :
3966 0 : CALL timestop(handle)
3967 :
3968 0 : END SUBROUTINE send_fock_to_client
3969 : #endif
3970 :
3971 0 : END MODULE qs_active_space_methods
|
