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