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