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