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