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