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