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