(git:ccc2433)
mixed_cdft_utils.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 Utility subroutines for mixed CDFT calculations
10 !> \par History
11 !> separated from mixed_cdft_methods [01.2017]
12 !> \author Nico Holmberg [01.2017]
13 ! **************************************************************************************************
14 MODULE mixed_cdft_utils
15  USE atomic_kind_types, ONLY: atomic_kind_type
16  USE cell_types, ONLY: cell_type
17  USE cp_array_utils, ONLY: cp_1d_i_p_type,&
18  cp_1d_r_p_type,&
19  cp_2d_r_p_type
20  USE cp_blacs_env, ONLY: cp_blacs_env_create,&
21  cp_blacs_env_type
22  USE cp_control_types, ONLY: dft_control_type
23  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
24  copy_fm_to_dbcsr_bc
25  USE cp_files, ONLY: open_file
26  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
27  cp_fm_struct_release,&
28  cp_fm_struct_type
29  USE cp_fm_types, ONLY: cp_fm_copy_general,&
30  cp_fm_create,&
31  cp_fm_get_info,&
32  cp_fm_release,&
33  cp_fm_to_fm,&
34  cp_fm_type
35  USE cp_log_handling, ONLY: cp_get_default_logger,&
36  cp_logger_create,&
37  cp_logger_set,&
38  cp_logger_type,&
39  cp_to_string
40  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
41  cp_print_key_unit_nr
42  USE cp_realspace_grid_init, ONLY: init_input_type
43  USE cp_subsys_types, ONLY: cp_subsys_get,&
44  cp_subsys_type
45  USE cube_utils, ONLY: init_cube_info,&
46  return_cube_max_iradius
47  USE d3_poly, ONLY: init_d3_poly_module
48  USE dbcsr_api, ONLY: dbcsr_desymmetrize,&
49  dbcsr_get_info,&
50  dbcsr_init_p,&
51  dbcsr_p_type,&
52  dbcsr_release,&
53  dbcsr_release_p,&
54  dbcsr_type
55  USE force_env_types, ONLY: force_env_get,&
56  force_env_type,&
57  multiple_fe_list
58  USE gaussian_gridlevels, ONLY: init_gaussian_gridlevel
59  USE global_types, ONLY: global_environment_type
60  USE hirshfeld_types, ONLY: create_hirshfeld_type,&
61  release_hirshfeld_type,&
62  set_hirshfeld_info
63  USE input_constants, ONLY: becke_cutoff_element,&
64  mixed_cdft_parallel,&
65  mixed_cdft_parallel_nobuild,&
66  mixed_cdft_serial,&
67  outer_scf_becke_constraint,&
68  outer_scf_hirshfeld_constraint,&
69  shape_function_gaussian
70  USE input_section_types, ONLY: section_vals_duplicate,&
71  section_vals_get,&
72  section_vals_get_subs_vals,&
73  section_vals_release,&
74  section_vals_type,&
75  section_vals_val_get
76  USE kinds, ONLY: default_path_length,&
77  dp
78  USE message_passing, ONLY: mp_request_type,&
79  mp_waitall
80  USE mixed_cdft_types, ONLY: mixed_cdft_result_type_release,&
81  mixed_cdft_result_type_set,&
82  mixed_cdft_settings_type,&
83  mixed_cdft_type,&
84  mixed_cdft_work_type_init
85  USE mixed_environment_types, ONLY: get_mixed_env,&
86  mixed_environment_type
87  USE pw_env_methods, ONLY: pw_env_create
88  USE pw_env_types, ONLY: pw_env_get,&
89  pw_env_type
90  USE pw_grid_types, ONLY: halfspace,&
91  pw_grid_type
92  USE pw_grids, ONLY: do_pw_grid_blocked_false,&
93  pw_grid_create,&
94  pw_grid_release,&
95  pw_grid_setup
96  USE pw_pool_types, ONLY: pw_pool_create,&
97  pw_pool_p_type,&
98  pw_pool_type
99  USE qs_cdft_types, ONLY: cdft_control_create,&
100  cdft_control_type
101  USE qs_environment_types, ONLY: get_qs_env,&
102  qs_environment_type
103  USE qs_kind_types, ONLY: create_qs_kind_set,&
104  qs_kind_type
105  USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
106  realspace_grid_input_type,&
107  realspace_grid_type,&
108  rs_grid_create,&
109  rs_grid_create_descriptor,&
110  rs_grid_print
111 #include "./base/base_uses.f90"
112 
113  IMPLICIT NONE
114  PRIVATE
115 
116 ! *** Public subroutines ***
117 
118  PUBLIC :: mixed_cdft_parse_settings, mixed_cdft_transfer_settings, &
119  mixed_cdft_init_structures, mixed_cdft_redistribute_arrays, &
120  mixed_cdft_print_couplings, map_permutation_to_states, hfun_zero, &
121  mixed_cdft_release_work, mixed_cdft_read_block_diag, &
122  mixed_cdft_get_blocks, mixed_cdft_diagonalize_blocks, &
123  mixed_cdft_assemble_block_diag
124 
125  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_utils'
126 
127 CONTAINS
128 
129 ! **************************************************************************************************
130 !> \brief Parse settings for mixed cdft calculation and check their consistency
131 !> \param force_env the force_env that holds the CDFT mixed_env
132 !> \param mixed_env the mixed_env that holds the CDFT states
133 !> \param mixed_cdft control section for mixed CDFT
134 !> \param settings container for settings related to the mixed CDFT calculation
135 !> \param natom the total number of atoms
136 !> \par History
137 !> 01.2017 created [Nico Holmberg]
138 ! **************************************************************************************************
139  SUBROUTINE mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, &
140  settings, natom)
141  TYPE(force_env_type), POINTER :: force_env
142  TYPE(mixed_environment_type), POINTER :: mixed_env
143  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
144  TYPE(mixed_cdft_settings_type) :: settings
145  INTEGER :: natom
146 
147  INTEGER :: i, iatom, iforce_eval, igroup, &
148  nforce_eval, nkinds
149  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: constraint_type
150  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: array_sizes
151  LOGICAL :: is_match
152  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
153  TYPE(cdft_control_type), POINTER :: cdft_control
154  TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:, :) :: atoms
155  TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: coeff
156  TYPE(dft_control_type), POINTER :: dft_control
157  TYPE(force_env_type), POINTER :: force_env_qs
158  TYPE(pw_env_type), POINTER :: pw_env
159  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
160  TYPE(qs_environment_type), POINTER :: qs_env
161 
162  NULLIFY (dft_control, qs_env, pw_env, auxbas_pw_pool, force_env_qs, &
163  cdft_control)
164  ! Allocate storage for temporaries used for checking settings consistency
165  settings%max_nkinds = 30
166  nforce_eval = SIZE(force_env%sub_force_env)
167  ALLOCATE (settings%grid_span(nforce_eval))
168  ALLOCATE (settings%npts(3, nforce_eval))
169  ALLOCATE (settings%cutoff(nforce_eval))
170  ALLOCATE (settings%rel_cutoff(nforce_eval))
171  ALLOCATE (settings%spherical(nforce_eval))
172  ALLOCATE (settings%rs_dims(2, nforce_eval))
173  ALLOCATE (settings%odd(nforce_eval))
174  ALLOCATE (settings%atoms(natom, nforce_eval))
175  IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
176  ALLOCATE (settings%coeffs(natom, nforce_eval))
177  settings%coeffs = 0.0_dp
178  END IF
179  ! Some of the checked settings are only defined for certain types of constraints
180  ! We nonetheless use arrays that are large enough to contain settings for all constraints
181  ! This is not completely optimal...
182  ALLOCATE (settings%si(6, nforce_eval))
183  ALLOCATE (settings%sb(8, nforce_eval))
184  ALLOCATE (settings%sr(5, nforce_eval))
185  ALLOCATE (settings%cutoffs(settings%max_nkinds, nforce_eval))
186  ALLOCATE (settings%radii(settings%max_nkinds, nforce_eval))
187  settings%grid_span = 0
188  settings%npts = 0
189  settings%cutoff = 0.0_dp
190  settings%rel_cutoff = 0.0_dp
191  settings%spherical = 0
192  settings%is_spherical = .false.
193  settings%rs_dims = 0
194  settings%odd = 0
195  settings%is_odd = .false.
196  settings%atoms = 0
197  settings%si = 0
198  settings%sr = 0.0_dp
199  settings%sb = .false.
200  settings%cutoffs = 0.0_dp
201  settings%radii = 0.0_dp
202  ! Get information from the sub_force_envs
203  DO iforce_eval = 1, nforce_eval
204  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
205  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
206  IF (mixed_env%do_mixed_qmmm_cdft) THEN
207  qs_env => force_env_qs%qmmm_env%qs_env
208  ELSE
209  CALL force_env_get(force_env_qs, qs_env=qs_env)
210  END IF
211  CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
212  IF (.NOT. dft_control%qs_control%cdft) &
213  CALL cp_abort(__location__, &
214  "A mixed CDFT simulation with multiple force_evals was requested, "// &
215  "but CDFT constraints were not active in the QS section of all force_evals!")
216  cdft_control => dft_control%qs_control%cdft_control
217  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
218  settings%bo = auxbas_pw_pool%pw_grid%bounds_local
219  ! Only the rank 0 process collects info about pw_grid and CDFT
220  IF (force_env_qs%para_env%is_source()) THEN
221  ! Grid settings
222  settings%grid_span(iforce_eval) = auxbas_pw_pool%pw_grid%grid_span
223  settings%npts(:, iforce_eval) = auxbas_pw_pool%pw_grid%npts
224  settings%cutoff(iforce_eval) = auxbas_pw_pool%pw_grid%cutoff
225  settings%rel_cutoff(iforce_eval) = dft_control%qs_control%relative_cutoff
226  IF (auxbas_pw_pool%pw_grid%spherical) settings%spherical(iforce_eval) = 1
227  settings%rs_dims(:, iforce_eval) = auxbas_pw_pool%pw_grid%para%rs_dims
228  IF (auxbas_pw_pool%pw_grid%grid_span == halfspace) settings%odd(iforce_eval) = 1
229  ! Becke constraint atoms/coeffs
230  IF (cdft_control%natoms .GT. SIZE(settings%atoms, 1)) &
231  CALL cp_abort(__location__, &
232  "More CDFT constraint atoms than defined in mixed section. "// &
233  "Use default values for MIXED\MAPPING.")
234  settings%atoms(1:cdft_control%natoms, iforce_eval) = cdft_control%atoms
235  IF (mixed_cdft%run_type == mixed_cdft_parallel) &
236  settings%coeffs(1:cdft_control%natoms, iforce_eval) = cdft_control%group(1)%coeff
237  ! Integer type settings
238  IF (cdft_control%type == outer_scf_becke_constraint) THEN
239  settings%si(1, iforce_eval) = cdft_control%becke_control%cutoff_type
240  settings%si(2, iforce_eval) = cdft_control%becke_control%cavity_shape
241  END IF
242  settings%si(3, iforce_eval) = dft_control%multiplicity
243  settings%si(4, iforce_eval) = SIZE(cdft_control%group)
244  settings%si(5, iforce_eval) = cdft_control%type
245  IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
246  settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%shape_function
247  settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%gaussian_shape
248  END IF
249  ! Logicals
250  IF (cdft_control%type == outer_scf_becke_constraint) THEN
251  settings%sb(1, iforce_eval) = cdft_control%becke_control%cavity_confine
252  settings%sb(2, iforce_eval) = cdft_control%becke_control%should_skip
253  settings%sb(3, iforce_eval) = cdft_control%becke_control%print_cavity
254  settings%sb(4, iforce_eval) = cdft_control%becke_control%in_memory
255  settings%sb(5, iforce_eval) = cdft_control%becke_control%adjust
256  settings%sb(8, iforce_eval) = cdft_control%becke_control%use_bohr
257  END IF
258  IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
259  settings%sb(8, iforce_eval) = cdft_control%hirshfeld_control%use_bohr
260  END IF
261  settings%sb(6, iforce_eval) = cdft_control%atomic_charges
262  settings%sb(7, iforce_eval) = qs_env%has_unit_metric
263  ! Reals
264  IF (cdft_control%type == outer_scf_becke_constraint) THEN
265  settings%sr(1, iforce_eval) = cdft_control%becke_control%rcavity
266  settings%sr(2, iforce_eval) = cdft_control%becke_control%rglobal
267  settings%sr(3, iforce_eval) = cdft_control%becke_control%eps_cavity
268  END IF
269  IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
270  settings%sr(2, iforce_eval) = cdft_control%hirshfeld_control%radius
271  END IF
272  settings%sr(4, iforce_eval) = dft_control%qs_control%eps_rho_rspace
273  settings%sr(5, iforce_eval) = pw_env%cube_info(pw_env%auxbas_grid)%max_rad_ga
274  IF (cdft_control%type == outer_scf_becke_constraint) THEN
275  IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
276  nkinds = SIZE(cdft_control%becke_control%cutoffs_tmp)
277  IF (nkinds .GT. settings%max_nkinds) &
278  CALL cp_abort(__location__, &
279  "More than "//trim(cp_to_string(settings%max_nkinds))// &
280  " unique elements were defined in BECKE_CONSTRAINT\ELEMENT_CUTOFF. Are you sure"// &
281  " your input is correct? If yes, please increase max_nkinds and recompile.")
282  settings%cutoffs(1:nkinds, iforce_eval) = cdft_control%becke_control%cutoffs_tmp(:)
283  END IF
284  IF (cdft_control%becke_control%adjust) THEN
285  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
286  IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%becke_control%radii_tmp)) &
287  CALL cp_abort(__location__, &
288  "Length of keyword BECKE_CONSTRAINT\ATOMIC_RADII does not "// &
289  "match number of atomic kinds in the input coordinate file.")
290  nkinds = SIZE(cdft_control%becke_control%radii_tmp)
291  IF (nkinds .GT. settings%max_nkinds) &
292  CALL cp_abort(__location__, &
293  "More than "//trim(cp_to_string(settings%max_nkinds))// &
294  " unique elements were defined in BECKE_CONSTRAINT\ATOMIC_RADII. Are you sure"// &
295  " your input is correct? If yes, please increase max_nkinds and recompile.")
296  settings%radii(1:nkinds, iforce_eval) = cdft_control%becke_control%radii_tmp(:)
297  END IF
298  END IF
299  IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
300  IF (ASSOCIATED(cdft_control%hirshfeld_control%radii)) THEN
301  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
302  IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%hirshfeld_control%radii)) &
303  CALL cp_abort(__location__, &
304  "Length of keyword HIRSHFELD_CONSTRAINT&RADII does not "// &
305  "match number of atomic kinds in the input coordinate file.")
306  nkinds = SIZE(cdft_control%hirshfeld_control%radii)
307  IF (nkinds .GT. settings%max_nkinds) &
308  CALL cp_abort(__location__, &
309  "More than "//trim(cp_to_string(settings%max_nkinds))// &
310  " unique elements were defined in HIRSHFELD_CONSTRAINT&RADII. Are you sure"// &
311  " your input is correct? If yes, please increase max_nkinds and recompile.")
312  settings%radii(1:nkinds, iforce_eval) = cdft_control%hirshfeld_control%radii(:)
313  END IF
314  END IF
315  END IF
316  END DO
317  ! Make sure the grids are consistent
318  CALL force_env%para_env%sum(settings%grid_span)
319  CALL force_env%para_env%sum(settings%npts)
320  CALL force_env%para_env%sum(settings%cutoff)
321  CALL force_env%para_env%sum(settings%rel_cutoff)
322  CALL force_env%para_env%sum(settings%spherical)
323  CALL force_env%para_env%sum(settings%rs_dims)
324  CALL force_env%para_env%sum(settings%odd)
325  is_match = .true.
326  DO iforce_eval = 2, nforce_eval
327  is_match = is_match .AND. (settings%grid_span(1) == settings%grid_span(iforce_eval))
328  is_match = is_match .AND. (settings%npts(1, 1) == settings%npts(1, iforce_eval))
329  is_match = is_match .AND. (settings%cutoff(1) == settings%cutoff(iforce_eval))
330  is_match = is_match .AND. (settings%rel_cutoff(1) == settings%rel_cutoff(iforce_eval))
331  is_match = is_match .AND. (settings%spherical(1) == settings%spherical(iforce_eval))
332  is_match = is_match .AND. (settings%rs_dims(1, 1) == settings%rs_dims(1, iforce_eval))
333  is_match = is_match .AND. (settings%rs_dims(2, 1) == settings%rs_dims(2, iforce_eval))
334  is_match = is_match .AND. (settings%odd(1) == settings%odd(iforce_eval))
335  END DO
336  IF (.NOT. is_match) &
337  CALL cp_abort(__location__, &
338  "Mismatch detected in the &MGRID settings of the CDFT force_evals.")
339  IF (settings%spherical(1) == 1) settings%is_spherical = .true.
340  IF (settings%odd(1) == 1) settings%is_odd = .true.
341  ! Make sure CDFT settings are consistent
342  CALL force_env%para_env%sum(settings%atoms)
343  IF (mixed_cdft%run_type == mixed_cdft_parallel) &
344  CALL force_env%para_env%sum(settings%coeffs)
345  settings%ncdft = 0
346  DO i = 1, SIZE(settings%atoms, 1)
347  DO iforce_eval = 2, nforce_eval
348  IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
349  IF (settings%atoms(i, 1) /= settings%atoms(i, iforce_eval)) is_match = .false.
350  IF (settings%coeffs(i, 1) /= settings%coeffs(i, iforce_eval)) is_match = .false.
351  END IF
352  END DO
353  IF (settings%atoms(i, 1) /= 0) settings%ncdft = settings%ncdft + 1
354  END DO
355  IF (.NOT. is_match .AND. mixed_cdft%run_type == mixed_cdft_parallel) &
356  CALL cp_abort(__location__, &
357  "Mismatch detected in the &CDFT section of the CDFT force_evals. "// &
358  "Parallel mode mixed CDFT requires identical constraint definitions in both CDFT states. "// &
359  "Switch to serial mode or disable keyword PARALLEL_BUILD if you "// &
360  "want to use nonidentical constraint definitions.")
361  CALL force_env%para_env%sum(settings%si)
362  CALL force_env%para_env%sum(settings%sr)
363  DO i = 1, SIZE(settings%sb, 1)
364  CALL force_env%para_env%sum(settings%sb(i, 1))
365  DO iforce_eval = 2, nforce_eval
366  CALL force_env%para_env%sum(settings%sb(i, iforce_eval))
367  IF (settings%sb(i, 1) .NEQV. settings%sb(i, iforce_eval)) is_match = .false.
368  END DO
369  END DO
370  DO i = 1, SIZE(settings%si, 1)
371  DO iforce_eval = 2, nforce_eval
372  IF (settings%si(i, 1) /= settings%si(i, iforce_eval)) is_match = .false.
373  END DO
374  END DO
375  DO i = 1, SIZE(settings%sr, 1)
376  DO iforce_eval = 2, nforce_eval
377  IF (settings%sr(i, 1) /= settings%sr(i, iforce_eval)) is_match = .false.
378  END DO
379  END DO
380  IF (.NOT. is_match) &
381  CALL cp_abort(__location__, &
382  "Mismatch detected in the &CDFT settings of the CDFT force_evals.")
383  ! Some CDFT features are currently disabled for mixed calculations: check that these features were not requested
384  IF (mixed_cdft%dlb .AND. .NOT. settings%sb(1, 1)) &
385  CALL cp_abort(__location__, &
386  "Parallel mode mixed CDFT load balancing requires Gaussian cavity confinement.")
387  ! Check for identical constraints in case of run type serial/parallel_nobuild
388  IF (mixed_cdft%run_type /= mixed_cdft_parallel) THEN
389  ! Get array sizes
390  ALLOCATE (array_sizes(nforce_eval, settings%si(4, 1), 2))
391  array_sizes(:, :, :) = 0
392  DO iforce_eval = 1, nforce_eval
393  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
394  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
395  IF (mixed_env%do_mixed_qmmm_cdft) THEN
396  qs_env => force_env_qs%qmmm_env%qs_env
397  ELSE
398  CALL force_env_get(force_env_qs, qs_env=qs_env)
399  END IF
400  CALL get_qs_env(qs_env, dft_control=dft_control)
401  cdft_control => dft_control%qs_control%cdft_control
402  IF (force_env_qs%para_env%is_source()) THEN
403  DO igroup = 1, SIZE(cdft_control%group)
404  array_sizes(iforce_eval, igroup, 1) = SIZE(cdft_control%group(igroup)%atoms)
405  array_sizes(iforce_eval, igroup, 2) = SIZE(cdft_control%group(igroup)%coeff)
406  END DO
407  END IF
408  END DO
409  ! Sum up array sizes and check consistency
410  CALL force_env%para_env%sum(array_sizes)
411  IF (any(array_sizes(:, :, 1) /= array_sizes(1, 1, 1)) .OR. &
412  any(array_sizes(:, :, 2) /= array_sizes(1, 1, 2))) &
413  mixed_cdft%identical_constraints = .false.
414  ! Check constraint definitions
415  IF (mixed_cdft%identical_constraints) THEN
416  ! Prepare temporary storage
417  ALLOCATE (atoms(nforce_eval, settings%si(4, 1)))
418  ALLOCATE (coeff(nforce_eval, settings%si(4, 1)))
419  ALLOCATE (constraint_type(nforce_eval, settings%si(4, 1)))
420  constraint_type(:, :) = 0
421  DO iforce_eval = 1, nforce_eval
422  DO i = 1, settings%si(4, 1)
423  NULLIFY (atoms(iforce_eval, i)%array)
424  NULLIFY (coeff(iforce_eval, i)%array)
425  ALLOCATE (atoms(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
426  ALLOCATE (coeff(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
427  atoms(iforce_eval, i)%array(:) = 0
428  coeff(iforce_eval, i)%array(:) = 0
429  END DO
430  ! Get constraint definitions
431  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
432  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
433  IF (mixed_env%do_mixed_qmmm_cdft) THEN
434  qs_env => force_env_qs%qmmm_env%qs_env
435  ELSE
436  CALL force_env_get(force_env_qs, qs_env=qs_env)
437  END IF
438  CALL get_qs_env(qs_env, dft_control=dft_control)
439  cdft_control => dft_control%qs_control%cdft_control
440  IF (force_env_qs%para_env%is_source()) THEN
441  DO i = 1, settings%si(4, 1)
442  atoms(iforce_eval, i)%array(:) = cdft_control%group(i)%atoms
443  coeff(iforce_eval, i)%array(:) = cdft_control%group(i)%coeff
444  constraint_type(iforce_eval, i) = cdft_control%group(i)%constraint_type
445  END DO
446  END IF
447  END DO
448  ! Sum up constraint definitions and check consistency
449  DO i = 1, settings%si(4, 1)
450  DO iforce_eval = 1, nforce_eval
451  CALL force_env%para_env%sum(atoms(iforce_eval, i)%array)
452  CALL force_env%para_env%sum(coeff(iforce_eval, i)%array)
453  CALL force_env%para_env%sum(constraint_type(iforce_eval, i))
454  END DO
455  DO iforce_eval = 2, nforce_eval
456  DO iatom = 1, SIZE(atoms(1, i)%array)
457  IF (atoms(1, i)%array(iatom) /= atoms(iforce_eval, i)%array(iatom)) &
458  mixed_cdft%identical_constraints = .false.
459  IF (coeff(1, i)%array(iatom) /= coeff(iforce_eval, i)%array(iatom)) &
460  mixed_cdft%identical_constraints = .false.
461  IF (.NOT. mixed_cdft%identical_constraints) EXIT
462  END DO
463  IF (constraint_type(1, i) /= constraint_type(iforce_eval, i)) &
464  mixed_cdft%identical_constraints = .false.
465  IF (.NOT. mixed_cdft%identical_constraints) EXIT
466  END DO
467  IF (.NOT. mixed_cdft%identical_constraints) EXIT
468  END DO
469  ! Deallocate temporary storage
470  DO iforce_eval = 1, nforce_eval
471  DO i = 1, settings%si(4, 1)
472  DEALLOCATE (atoms(iforce_eval, i)%array)
473  DEALLOCATE (coeff(iforce_eval, i)%array)
474  END DO
475  END DO
476  DEALLOCATE (atoms)
477  DEALLOCATE (coeff)
478  DEALLOCATE (constraint_type)
479  END IF
480  DEALLOCATE (array_sizes)
481  END IF
482  ! Deallocate some arrays that are no longer needed
483  IF (mixed_cdft%identical_constraints .AND. mixed_cdft%run_type /= mixed_cdft_parallel_nobuild) THEN
484  DO iforce_eval = 1, nforce_eval
485  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
486  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
487  IF (mixed_env%do_mixed_qmmm_cdft) THEN
488  qs_env => force_env_qs%qmmm_env%qs_env
489  ELSE
490  CALL force_env_get(force_env_qs, qs_env=qs_env)
491  END IF
492  CALL get_qs_env(qs_env, dft_control=dft_control)
493  cdft_control => dft_control%qs_control%cdft_control
494  IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
495  IF (.NOT. dft_control%qs_control%gapw) THEN
496  DO i = 1, SIZE(cdft_control%group)
497  DEALLOCATE (cdft_control%group(i)%coeff)
498  DEALLOCATE (cdft_control%group(i)%atoms)
499  END DO
500  IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
501  END IF
502  ELSE IF (mixed_cdft%run_type == mixed_cdft_serial) THEN
503  IF (iforce_eval == 1) cycle
504  DO igroup = 1, SIZE(cdft_control%group)
505  IF (.NOT. dft_control%qs_control%gapw) THEN
506  DEALLOCATE (cdft_control%group(igroup)%coeff)
507  DEALLOCATE (cdft_control%group(igroup)%atoms)
508  END IF
509  END DO
510  IF (cdft_control%type == outer_scf_becke_constraint) THEN
511  IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
512  IF (cdft_control%becke_control%cavity_confine) &
513  CALL release_hirshfeld_type(cdft_control%becke_control%cavity_env)
514  IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) &
515  DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
516  IF (cdft_control%becke_control%adjust) &
517  DEALLOCATE (cdft_control%becke_control%radii_tmp)
518  END IF
519  END IF
520  END DO
521  END IF
522 
523  END SUBROUTINE mixed_cdft_parse_settings
524 
525 ! **************************************************************************************************
526 !> \brief Transfer settings to mixed_cdft
527 !> \param force_env the force_env that holds the CDFT states
528 !> \param mixed_cdft the control section for mixed CDFT calculations
529 !> \param settings container for settings related to the mixed CDFT calculation
530 !> \par History
531 !> 01.2017 created [Nico Holmberg]
532 ! **************************************************************************************************
533  SUBROUTINE mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
534  TYPE(force_env_type), POINTER :: force_env
535  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
536  TYPE(mixed_cdft_settings_type) :: settings
537 
538  INTEGER :: i, nkinds
539  LOGICAL :: is_match
540  TYPE(cdft_control_type), POINTER :: cdft_control
541 
542  NULLIFY (cdft_control)
543  is_match = .true.
544  ! Transfer global settings
545  mixed_cdft%multiplicity = settings%si(3, 1)
546  mixed_cdft%has_unit_metric = settings%sb(7, 1)
547  mixed_cdft%eps_rho_rspace = settings%sr(4, 1)
548  mixed_cdft%nconstraint = settings%si(4, 1)
549  settings%radius = settings%sr(5, 1)
550  ! Transfer settings only needed if the constraint should be built in parallel
551  IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
552  IF (settings%sb(6, 1)) &
553  CALL cp_abort(__location__, &
554  "Calculation of atomic Becke charges not supported with parallel mode mixed CDFT")
555  IF (mixed_cdft%nconstraint /= 1) &
556  CALL cp_abort(__location__, &
557  "Parallel mode mixed CDFT does not yet support multiple constraints.")
558 
559  IF (settings%si(5, 1) /= outer_scf_becke_constraint) &
560  CALL cp_abort(__location__, &
561  "Parallel mode mixed CDFT does not support Hirshfeld constraints.")
562 
563  ALLOCATE (mixed_cdft%cdft_control)
564  CALL cdft_control_create(mixed_cdft%cdft_control)
565  cdft_control => mixed_cdft%cdft_control
566  ALLOCATE (cdft_control%atoms(settings%ncdft))
567  cdft_control%atoms = settings%atoms(1:settings%ncdft, 1)
568  ALLOCATE (cdft_control%group(1))
569  ALLOCATE (cdft_control%group(1)%atoms(settings%ncdft))
570  ALLOCATE (cdft_control%group(1)%coeff(settings%ncdft))
571  NULLIFY (cdft_control%group(1)%weight)
572  NULLIFY (cdft_control%group(1)%gradients)
573  NULLIFY (cdft_control%group(1)%integrated)
574  cdft_control%group(1)%atoms = cdft_control%atoms
575  cdft_control%group(1)%coeff = settings%coeffs(1:settings%ncdft, 1)
576  cdft_control%natoms = settings%ncdft
577  cdft_control%atomic_charges = settings%sb(6, 1)
578  cdft_control%becke_control%cutoff_type = settings%si(1, 1)
579  cdft_control%becke_control%cavity_confine = settings%sb(1, 1)
580  cdft_control%becke_control%should_skip = settings%sb(2, 1)
581  cdft_control%becke_control%print_cavity = settings%sb(3, 1)
582  cdft_control%becke_control%in_memory = settings%sb(4, 1)
583  cdft_control%becke_control%adjust = settings%sb(5, 1)
584  cdft_control%becke_control%cavity_shape = settings%si(2, 1)
585  cdft_control%becke_control%use_bohr = settings%sb(8, 1)
586  cdft_control%becke_control%rcavity = settings%sr(1, 1)
587  cdft_control%becke_control%rglobal = settings%sr(2, 1)
588  cdft_control%becke_control%eps_cavity = settings%sr(3, 1)
589  nkinds = 0
590  IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
591  CALL force_env%para_env%sum(settings%cutoffs)
592  DO i = 1, SIZE(settings%cutoffs, 1)
593  IF (settings%cutoffs(i, 1) /= settings%cutoffs(i, 2)) is_match = .false.
594  IF (settings%cutoffs(i, 1) /= 0.0_dp) nkinds = nkinds + 1
595  END DO
596  IF (.NOT. is_match) &
597  CALL cp_abort(__location__, &
598  "Mismatch detected in the &BECKE_CONSTRAINT "// &
599  "&ELEMENT_CUTOFF settings of the two force_evals.")
600  ALLOCATE (cdft_control%becke_control%cutoffs_tmp(nkinds))
601  cdft_control%becke_control%cutoffs_tmp = settings%cutoffs(1:nkinds, 1)
602  END IF
603  nkinds = 0
604  IF (cdft_control%becke_control%adjust) THEN
605  CALL force_env%para_env%sum(settings%radii)
606  DO i = 1, SIZE(settings%radii, 1)
607  IF (settings%radii(i, 1) /= settings%radii(i, 2)) is_match = .false.
608  IF (settings%radii(i, 1) /= 0.0_dp) nkinds = nkinds + 1
609  END DO
610  IF (.NOT. is_match) &
611  CALL cp_abort(__location__, &
612  "Mismatch detected in the &BECKE_CONSTRAINT "// &
613  "&ATOMIC_RADII settings of the two force_evals.")
614  ALLOCATE (cdft_control%becke_control%radii(nkinds))
615  cdft_control%becke_control%radii = settings%radii(1:nkinds, 1)
616  END IF
617  END IF
618 
619  END SUBROUTINE mixed_cdft_transfer_settings
620 
621 ! **************************************************************************************************
622 !> \brief Initialize all the structures needed for a mixed CDFT calculation
623 !> \param force_env the force_env that holds the CDFT mixed_env
624 !> \param force_env_qs the force_env that holds the qs_env, which is CDFT state specific
625 !> \param mixed_env the mixed_env that holds the CDFT states
626 !> \param mixed_cdft the control section for mixed CDFT calculations
627 !> \param settings container for settings related to the mixed CDFT calculation
628 !> \par History
629 !> 01.2017 created [Nico Holmberg]
630 ! **************************************************************************************************
631  SUBROUTINE mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
632  TYPE(force_env_type), POINTER :: force_env, force_env_qs
633  TYPE(mixed_environment_type), POINTER :: mixed_env
634  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
635  TYPE(mixed_cdft_settings_type) :: settings
636 
637  CHARACTER(len=default_path_length) :: c_val, input_file_path, output_file_path
638  INTEGER :: i, imap, iounit, j, lp, n_force_eval, &
639  ncpu, nforce_eval, ntargets, offset, &
640  unit_nr
641  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds
642  INTEGER, DIMENSION(2, 3) :: bo, bo_mixed
643  INTEGER, DIMENSION(3) :: higher_grid_layout
644  INTEGER, DIMENSION(:), POINTER :: i_force_eval, mixed_rs_dims, recvbuffer, &
645  recvbuffer2, sendbuffer
646  LOGICAL :: is_match
647  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
648  TYPE(cell_type), POINTER :: cell_mix
649  TYPE(cp_logger_type), POINTER :: logger
650  TYPE(cp_subsys_type), POINTER :: subsys_mix
651  TYPE(global_environment_type), POINTER :: globenv
652  TYPE(mp_request_type), DIMENSION(3) :: req
653  TYPE(pw_env_type), POINTER :: pw_env
654  TYPE(pw_grid_type), POINTER :: pw_grid
655  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
656  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
657  TYPE(qs_environment_type), POINTER :: qs_env
658  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
659  TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
660  POINTER :: rs_descs
661  TYPE(realspace_grid_input_type) :: input_settings
662  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
663  TYPE(section_vals_type), POINTER :: force_env_section, force_env_sections, kind_section, &
664  print_section, root_section, rs_grid_section, subsys_section
665 
666  NULLIFY (cell_mix, subsys_mix, force_env_section, subsys_section, &
667  print_section, root_section, kind_section, force_env_sections, &
668  rs_grid_section, auxbas_pw_pool, pw_env, pw_pools, pw_grid, &
669  sendbuffer, qs_env, mixed_rs_dims, i_force_eval, recvbuffer, &
670  recvbuffer2, globenv, atomic_kind_set, qs_kind_set, rs_descs, &
671  rs_grids)
672 
673  logger => cp_get_default_logger()
674  CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
675  print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
676  iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
677  is_match = .true.
678  nforce_eval = SIZE(force_env%sub_force_env)
679  ncpu = force_env%para_env%num_pe
680  ! Get infos about the mixed subsys
681  IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
682  CALL force_env_get(force_env=force_env, &
683  subsys=subsys_mix)
684  ELSE
685  CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
686  cp_subsys=subsys_mix)
687  END IF
688  ! Init structures only needed when the CDFT states are treated in parallel
689  IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
690  ! Start building the mixed auxbas_pw_pool
691  CALL pw_grid_create(pw_grid, force_env%para_env)
692  CALL pw_env_create(mixed_cdft%pw_env)
693  ! Decide what kind of layout to use and setup the grid
694  ! Processor mappings currently supported:
695  ! (2np,1) --> (np,1)
696  ! (nx,2ny) --> (nx,ny)
697  ! (nx,ny) --> (nx*ny/2,1) (required when xc_smooth is in use and with intermediate proc counts)
698  !
699  ! For cases 2 and 3, dlb redistributes YZ slices from overloaded processors to underloaded processors
700  ! For case 1, XZ slices are redistributed
701  ! TODO: Unify mappings. Now we essentially have separate code for cases 1-2 and 3.
702  ! This leads to very messy code especially with dlb turned on...
703  ! In terms of memory usage, it would be beneficial to replace case 1 with 3
704  ! and implement a similar arbitrary mapping to replace case 2
705 
706  mixed_cdft%is_pencil = .false. ! Flag to control the first two mappings
707  mixed_cdft%is_special = .false. ! Flag to control the last mapping
708  ! With xc smoothing, the grid is always (ncpu/2,1) distributed
709  ! and correct behavior cannot be guaranteed for ncpu/2 > nx, so we abort...
710  IF (ncpu/2 .GT. settings%npts(1, 1)) &
711  cpabort("ncpu/2 => nx: decrease ncpu or disable xc_smoothing")
712  !
713  ALLOCATE (mixed_rs_dims(2))
714  IF (settings%rs_dims(2, 1) /= 1) mixed_cdft%is_pencil = .true.
715  IF (.NOT. mixed_cdft%is_pencil .AND. ncpu .GT. settings%npts(1, 1)) mixed_cdft%is_special = .true.
716  IF (mixed_cdft%is_special) THEN
717  mixed_rs_dims = (/-1, -1/)
718  ELSE IF (mixed_cdft%is_pencil) THEN
719  mixed_rs_dims = (/settings%rs_dims(1, 1), 2*settings%rs_dims(2, 1)/)
720  ELSE
721  mixed_rs_dims = (/2*settings%rs_dims(1, 1), 1/)
722  END IF
723  IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
724  CALL force_env_get(force_env=force_env, &
725  cell=cell_mix)
726  ELSE
727  CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
728  cell=cell_mix)
729  END IF
730  CALL pw_grid_setup(cell_mix%hmat, pw_grid, grid_span=settings%grid_span(1), &
731  npts=settings%npts(:, 1), cutoff=settings%cutoff(1), &
732  spherical=settings%is_spherical, odd=settings%is_odd, &
733  fft_usage=.true., ncommensurate=0, icommensurate=1, &
734  blocked=do_pw_grid_blocked_false, rs_dims=mixed_rs_dims, &
735  iounit=iounit)
736  ! Check if the layout was successfully created
737  IF (mixed_cdft%is_special) THEN
738  IF (.NOT. pw_grid%para%rs_dims(2) /= 1) is_match = .false.
739  ELSE IF (mixed_cdft%is_pencil) THEN
740  IF (.NOT. pw_grid%para%rs_dims(1) == mixed_rs_dims(1)) is_match = .false.
741  ELSE
742  IF (.NOT. pw_grid%para%rs_dims(2) == 1) is_match = .false.
743  END IF
744  IF (.NOT. is_match) &
745  CALL cp_abort(__location__, &
746  "Unable to create a suitable grid distribution "// &
747  "for mixed CDFT calculations. Try decreasing the total number "// &
748  "of processors or disabling xc_smoothing.")
749  DEALLOCATE (mixed_rs_dims)
750  ! Create the pool
751  bo_mixed = pw_grid%bounds_local
752  ALLOCATE (pw_pools(1))
753  NULLIFY (pw_pools(1)%pool)
754  CALL pw_pool_create(pw_pools(1)%pool, pw_grid=pw_grid)
755  ! Initialize Gaussian cavity confinement
756  IF (mixed_cdft%cdft_control%becke_control%cavity_confine) THEN
757  CALL create_hirshfeld_type(mixed_cdft%cdft_control%becke_control%cavity_env)
758  CALL set_hirshfeld_info(mixed_cdft%cdft_control%becke_control%cavity_env, &
759  shape_function_type=shape_function_gaussian, iterative=.false., &
760  radius_type=mixed_cdft%cdft_control%becke_control%cavity_shape, &
761  use_bohr=mixed_cdft%cdft_control%becke_control%use_bohr)
762  END IF
763  ! Gaussian confinement/wavefunction overlap method needs qs_kind_set
764  ! Gaussian cavity confinement also needs the auxbas_rs_grid
765  IF (mixed_cdft%cdft_control%becke_control%cavity_confine .OR. &
766  mixed_cdft%wfn_overlap_method) THEN
767  print_section => section_vals_get_subs_vals(force_env_section, &
768  "PRINT%GRID_INFORMATION")
769  ALLOCATE (mixed_cdft%pw_env%gridlevel_info)
770  CALL init_gaussian_gridlevel(mixed_cdft%pw_env%gridlevel_info, &
771  ngrid_levels=1, cutoff=settings%cutoff, &
772  rel_cutoff=settings%rel_cutoff(1), &
773  print_section=print_section)
774  ALLOCATE (rs_descs(1))
775  ALLOCATE (rs_grids(1))
776  ALLOCATE (mixed_cdft%pw_env%cube_info(1))
777  higher_grid_layout = (/-1, -1, -1/)
778  CALL init_d3_poly_module()
779  CALL init_cube_info(mixed_cdft%pw_env%cube_info(1), &
780  pw_grid%dr(:), pw_grid%dh(:, :), &
781  pw_grid%dh_inv(:, :), &
782  pw_grid%orthorhombic, settings%radius)
783  NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
784  CALL force_env_get(force_env, root_section=root_section)
785  force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
786  CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, n_force_eval)
787  CALL section_vals_duplicate(force_env_sections, force_env_section, &
788  i_force_eval(2), i_force_eval(2))
789  rs_grid_section => section_vals_get_subs_vals(force_env_section, "DFT%MGRID%RS_GRID")
790  CALL init_input_type(input_settings, &
791  nsmax=2*max(1, return_cube_max_iradius(mixed_cdft%pw_env%cube_info(1))) + 1, &
792  rs_grid_section=rs_grid_section, ilevel=1, &
793  higher_grid_layout=higher_grid_layout)
794  NULLIFY (rs_descs(1)%rs_desc)
795  CALL rs_grid_create_descriptor(rs_descs(1)%rs_desc, pw_grid, input_settings)
796  IF (rs_descs(1)%rs_desc%distributed) higher_grid_layout = rs_descs(1)%rs_desc%group_dim
797  CALL rs_grid_create(rs_grids(1), rs_descs(1)%rs_desc)
798  CALL rs_grid_print(rs_grids(1), iounit)
799  mixed_cdft%pw_env%rs_descs => rs_descs
800  mixed_cdft%pw_env%rs_grids => rs_grids
801  ! qs_kind_set
802  subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
803  i_rep_section=i_force_eval(1))
804  kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
805  NULLIFY (qs_kind_set)
806  CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
807  CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
808  force_env%para_env, force_env_section)
809  mixed_cdft%qs_kind_set => qs_kind_set
810  DEALLOCATE (i_force_eval)
811  CALL section_vals_release(force_env_section)
812  END IF
813  CALL force_env_get(force_env=force_env, &
814  force_env_section=force_env_section)
815  CALL pw_grid_release(pw_grid)
816  mixed_cdft%pw_env%auxbas_grid = 1
817  NULLIFY (mixed_cdft%pw_env%pw_pools)
818  mixed_cdft%pw_env%pw_pools => pw_pools
819  bo = settings%bo
820  ! Determine which processors need to exchange data when redistributing the weight/gradient
821  IF (.NOT. mixed_cdft%is_special) THEN
822  ALLOCATE (mixed_cdft%dest_list(2))
823  ALLOCATE (mixed_cdft%source_list(2))
824  imap = force_env%para_env%mepos/2
825  mixed_cdft%dest_list = (/imap, imap + force_env%para_env%num_pe/2/)
826  imap = mod(force_env%para_env%mepos, force_env%para_env%num_pe/2) + &
827  modulo(force_env%para_env%mepos, force_env%para_env%num_pe/2)
828  mixed_cdft%source_list = (/imap, imap + 1/)
829  ! Determine bounds of the data that is replicated
830  ALLOCATE (mixed_cdft%recv_bo(4))
831  ALLOCATE (sendbuffer(2), recvbuffer(2), recvbuffer2(2))
832  IF (mixed_cdft%is_pencil) THEN
833  sendbuffer = (/bo_mixed(1, 2), bo_mixed(2, 2)/)
834  ELSE
835  sendbuffer = (/bo_mixed(1, 1), bo_mixed(2, 1)/)
836  END IF
837  ! Communicate bounds in steps
838  CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dest_list(1), &
839  request=req(1))
840  CALL force_env%para_env%irecv(msgout=recvbuffer, source=mixed_cdft%source_list(1), &
841  request=req(2))
842  CALL force_env%para_env%irecv(msgout=recvbuffer2, source=mixed_cdft%source_list(2), &
843  request=req(3))
844  CALL req(1)%wait()
845  CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dest_list(2), &
846  request=req(1))
847  CALL mp_waitall(req)
848  mixed_cdft%recv_bo(1:2) = recvbuffer
849  mixed_cdft%recv_bo(3:4) = recvbuffer2
850  DEALLOCATE (sendbuffer, recvbuffer, recvbuffer2)
851  ELSE
852  IF (mixed_env%do_mixed_qmmm_cdft) THEN
853  qs_env => force_env_qs%qmmm_env%qs_env
854  ELSE
855  CALL force_env_get(force_env_qs, qs_env=qs_env)
856  END IF
857  CALL get_qs_env(qs_env, pw_env=pw_env)
858  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
859  ! work out the pw grid points each proc holds in the two (identical) parallel proc groups
860  ! note we only care about the x dir since we assume the y dir is not subdivided
861  ALLOCATE (bounds(0:auxbas_pw_pool%pw_grid%para%group_size - 1, 1:2))
862  DO i = 0, auxbas_pw_pool%pw_grid%para%group_size - 1
863  bounds(i, 1:2) = auxbas_pw_pool%pw_grid%para%bo(1:2, 1, i, 1)
864  bounds(i, 1:2) = bounds(i, 1:2) - auxbas_pw_pool%pw_grid%npts(1)/2 - 1
865  END DO
866  ! work out which procs to send my grid points
867  ! first get the number of target procs per group
868  ntargets = 0
869  offset = -1
870  DO i = 0, auxbas_pw_pool%pw_grid%para%group_size - 1
871  IF ((bounds(i, 1) .GE. bo_mixed(1, 1) .AND. bounds(i, 1) .LE. bo_mixed(2, 1)) .OR. &
872  (bounds(i, 2) .GE. bo_mixed(1, 1) .AND. bounds(i, 2) .LE. bo_mixed(2, 1))) THEN
873  ntargets = ntargets + 1
874  IF (offset == -1) offset = i
875  ELSE IF (bounds(i, 2) .GT. bo_mixed(2, 1)) THEN
876  EXIT
877  ELSE
878  cycle
879  END IF
880  END DO
881  ALLOCATE (mixed_cdft%dest_list(ntargets))
882  ALLOCATE (mixed_cdft%dest_list_bo(2, ntargets))
883  ! now determine the actual grid points to send
884  j = 1
885  DO i = offset, offset + ntargets - 1
886  mixed_cdft%dest_list(j) = i
887  mixed_cdft%dest_list_bo(:, j) = (/bo_mixed(1, 1) + (bounds(i, 1) - bo_mixed(1, 1)), &
888  bo_mixed(2, 1) + (bounds(i, 2) - bo_mixed(2, 1))/)
889  j = j + 1
890  END DO
891  ALLOCATE (mixed_cdft%dest_list_save(ntargets), mixed_cdft%dest_bo_save(2, ntargets))
892  ! We need to store backups of these arrays since they might get reallocated during dlb
893  mixed_cdft%dest_list_save = mixed_cdft%dest_list
894  mixed_cdft%dest_bo_save = mixed_cdft%dest_list_bo
895  ! finally determine which procs will send me grid points
896  ! now we need info about y dir also
897  DEALLOCATE (bounds)
898  ALLOCATE (bounds(0:pw_pools(1)%pool%pw_grid%para%group_size - 1, 1:4))
899  DO i = 0, pw_pools(1)%pool%pw_grid%para%group_size - 1
900  bounds(i, 1:2) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 1, i, 1)
901  bounds(i, 3:4) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 2, i, 1)
902  bounds(i, 1:2) = bounds(i, 1:2) - pw_pools(1)%pool%pw_grid%npts(1)/2 - 1
903  bounds(i, 3:4) = bounds(i, 3:4) - pw_pools(1)%pool%pw_grid%npts(2)/2 - 1
904  END DO
905  ntargets = 0
906  offset = -1
907  DO i = 0, pw_pools(1)%pool%pw_grid%para%group_size - 1
908  IF ((bo(1, 1) .GE. bounds(i, 1) .AND. bo(1, 1) .LE. bounds(i, 2)) .OR. &
909  (bo(2, 1) .GE. bounds(i, 1) .AND. bo(2, 1) .LE. bounds(i, 2))) THEN
910  ntargets = ntargets + 1
911  IF (offset == -1) offset = i
912  ELSE IF (bo(2, 1) .LT. bounds(i, 1)) THEN
913  EXIT
914  ELSE
915  cycle
916  END IF
917  END DO
918  ALLOCATE (mixed_cdft%source_list(ntargets))
919  ALLOCATE (mixed_cdft%source_list_bo(4, ntargets))
920  j = 1
921  DO i = offset, offset + ntargets - 1
922  mixed_cdft%source_list(j) = i
923  IF (bo(1, 1) .GE. bounds(i, 1) .AND. bo(2, 1) .LE. bounds(i, 2)) THEN
924  mixed_cdft%source_list_bo(:, j) = (/bo(1, 1), bo(2, 1), &
925  bounds(i, 3), bounds(i, 4)/)
926  ELSE IF (bo(1, 1) .GE. bounds(i, 1) .AND. bo(1, 1) .LE. bounds(i, 2)) THEN
927  mixed_cdft%source_list_bo(:, j) = (/bo(1, 1), bounds(i, 2), &
928  bounds(i, 3), bounds(i, 4)/)
929  ELSE
930  mixed_cdft%source_list_bo(:, j) = (/bounds(i, 1), bo(2, 1), &
931  bounds(i, 3), bounds(i, 4)/)
932  END IF
933  j = j + 1
934  END DO
935  ALLOCATE (mixed_cdft%source_list_save(ntargets), mixed_cdft%source_bo_save(4, ntargets))
936  ! We need to store backups of these arrays since they might get reallocated during dlb
937  mixed_cdft%source_list_save = mixed_cdft%source_list
938  mixed_cdft%source_bo_save = mixed_cdft%source_list_bo
939  DEALLOCATE (bounds)
940  END IF
941  ELSE
942  ! Create loggers to redirect the output of all CDFT states to different files
943  ! even when the states are treated in serial (the initial print of QS data [basis set etc] for
944  ! all states unfortunately goes to the first log file)
945  CALL force_env_get(force_env, root_section=root_section)
946  ALLOCATE (mixed_cdft%sub_logger(nforce_eval - 1))
947  DO i = 1, nforce_eval - 1
948  IF (force_env%para_env%is_source()) THEN
949  CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", &
950  c_val=input_file_path)
951  lp = len_trim(input_file_path)
952  input_file_path(lp + 1:len(input_file_path)) = "-r-"//adjustl(cp_to_string(i + 1))
953  lp = len_trim(input_file_path)
954  output_file_path = input_file_path(1:lp)//".out"
955  CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
956  file_action="WRITE", file_position="APPEND", &
957  unit_number=unit_nr)
958  ELSE
959  unit_nr = -1
960  END IF
961  CALL cp_logger_create(mixed_cdft%sub_logger(i)%p, &
962  para_env=force_env%para_env, &
963  default_global_unit_nr=unit_nr, &
964  close_global_unit_on_dealloc=.false.)
965  ! Try to use better names for the local log if it is not too late
966  CALL section_vals_val_get(root_section, "GLOBAL%OUTPUT_FILE_NAME", &
967  c_val=c_val)
968  IF (c_val /= "") THEN
969  CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
970  local_filename=trim(c_val)//"_localLog")
971  END IF
972  CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
973  IF (c_val /= "") THEN
974  CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
975  local_filename=trim(c_val)//"_localLog")
976  END IF
977  mixed_cdft%sub_logger(i)%p%iter_info%project_name = c_val
978  CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
979  i_val=mixed_cdft%sub_logger(i)%p%iter_info%print_level)
980  END DO
981  IF (mixed_cdft%wfn_overlap_method) THEN
982  ! qs_kind_set
983  NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
984  CALL force_env_get(force_env, root_section=root_section)
985  force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
986  CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, n_force_eval)
987  CALL section_vals_duplicate(force_env_sections, force_env_section, &
988  i_force_eval(2), i_force_eval(2))
989  subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
990  i_rep_section=i_force_eval(1))
991  kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
992  NULLIFY (qs_kind_set)
993  CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
994  CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
995  force_env%para_env, force_env_section)
996  mixed_cdft%qs_kind_set => qs_kind_set
997  DEALLOCATE (i_force_eval)
998  CALL section_vals_release(force_env_section)
999  mixed_cdft%qs_kind_set => qs_kind_set
1000  END IF
1001  CALL force_env_get(force_env=force_env, &
1002  force_env_section=force_env_section)
1003  END IF
1004  ! Deallocate settings temporaries
1005  DEALLOCATE (settings%grid_span)
1006  DEALLOCATE (settings%npts)
1007  DEALLOCATE (settings%spherical)
1008  DEALLOCATE (settings%rs_dims)
1009  DEALLOCATE (settings%odd)
1010  DEALLOCATE (settings%atoms)
1011  IF (mixed_cdft%run_type == mixed_cdft_parallel) &
1012  DEALLOCATE (settings%coeffs)
1013  DEALLOCATE (settings%cutoffs)
1014  DEALLOCATE (settings%radii)
1015  DEALLOCATE (settings%si)
1016  DEALLOCATE (settings%sr)
1017  DEALLOCATE (settings%sb)
1018  DEALLOCATE (settings%cutoff)
1019  DEALLOCATE (settings%rel_cutoff)
1020  ! Setup mixed blacs_env for redistributing arrays during ET coupling calculation
1021  IF (mixed_env%do_mixed_et) THEN
1022  NULLIFY (root_section)
1023  CALL force_env_get(force_env, globenv=globenv, root_section=root_section)
1024  CALL cp_blacs_env_create(mixed_cdft%blacs_env, force_env%para_env, globenv%blacs_grid_layout, &
1025  globenv%blacs_repeatable)
1026  END IF
1027  CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1028  "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1029 
1030  END SUBROUTINE mixed_cdft_init_structures
1031 
1032 ! **************************************************************************************************
1033 !> \brief Redistribute arrays needed for an ET coupling calculation from individual CDFT states to
1034 !> the mixed CDFT env, that is, move the arrays to the correct blacs context. For parallel
1035 !> simulations, the array processor distributions also change from N to 2N processors.
1036 !> \param force_env the force_env that holds the CDFT states
1037 !> \par History
1038 !> 01.2017 created [Nico Holmberg]
1039 ! **************************************************************************************************
1040  SUBROUTINE mixed_cdft_redistribute_arrays(force_env)
1041  TYPE(force_env_type), POINTER :: force_env
1042 
1043  INTEGER :: iforce_eval, ispin, ivar, ncol_overlap, &
1044  ncol_wmat, nforce_eval, nrow_overlap, &
1045  nrow_wmat, nspins, nvar
1046  INTEGER, ALLOCATABLE, DIMENSION(:) :: ncol_mo, nrow_mo
1047  LOGICAL :: uniform_occupation
1048  LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_occupation_numbers
1049  TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: occno_tmp
1050  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1051  TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, fm_struct_overlap, &
1052  fm_struct_tmp, fm_struct_wmat
1053  TYPE(cp_fm_type) :: matrix_s_tmp, mixed_matrix_s_tmp
1054  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: matrix_p_tmp, mixed_matrix_p_tmp, &
1055  mixed_wmat_tmp, mo_coeff_tmp, wmat_tmp
1056  TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
1057  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix, w_matrix
1058  TYPE(dbcsr_type) :: desymm_tmp
1059  TYPE(dbcsr_type), POINTER :: mixed_matrix_s
1060  TYPE(dft_control_type), POINTER :: dft_control
1061  TYPE(force_env_type), POINTER :: force_env_qs
1062  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1063  TYPE(mixed_environment_type), POINTER :: mixed_env
1064  TYPE(qs_environment_type), POINTER :: qs_env
1065 
1066  NULLIFY (mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
1067  fm_struct_wmat, fm_struct_overlap, fm_struct_tmp, &
1068  mixed_mo_coeff, mixed_matrix_s, density_matrix, blacs_env, w_matrix, force_env_qs)
1069  cpassert(ASSOCIATED(force_env))
1070  mixed_env => force_env%mixed_env
1071  nforce_eval = SIZE(force_env%sub_force_env)
1072  CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1073  cpassert(ASSOCIATED(mixed_cdft))
1074  CALL mixed_cdft_work_type_init(mixed_cdft%matrix)
1075  ! Get nspins and query for non-uniform occupation numbers
1076  ALLOCATE (has_occupation_numbers(nforce_eval))
1077  has_occupation_numbers = .false.
1078  DO iforce_eval = 1, nforce_eval
1079  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1080  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
1081  IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1082  qs_env => force_env_qs%qmmm_env%qs_env
1083  ELSE
1084  CALL force_env_get(force_env_qs, qs_env=qs_env)
1085  END IF
1086  CALL get_qs_env(qs_env, dft_control=dft_control)
1087  cpassert(ASSOCIATED(dft_control))
1088  nspins = dft_control%nspins
1089  IF (force_env_qs%para_env%is_source()) &
1090  has_occupation_numbers(iforce_eval) = ALLOCATED(dft_control%qs_control%cdft_control%occupations)
1091  END DO
1092  CALL force_env%para_env%sum(has_occupation_numbers(1))
1093  DO iforce_eval = 2, nforce_eval
1094  CALL force_env%para_env%sum(has_occupation_numbers(iforce_eval))
1095  IF (has_occupation_numbers(1) .NEQV. has_occupation_numbers(iforce_eval)) &
1096  CALL cp_abort(__location__, &
1097  "Mixing of uniform and non-uniform occupations is not allowed.")
1098  END DO
1099  uniform_occupation = .NOT. has_occupation_numbers(1)
1100  DEALLOCATE (has_occupation_numbers)
1101  ! Get number of weight functions per state as well as the type of each constraint
1102  nvar = SIZE(dft_control%qs_control%cdft_control%target)
1103  IF (.NOT. ALLOCATED(mixed_cdft%constraint_type)) THEN
1104  ALLOCATE (mixed_cdft%constraint_type(nvar, nforce_eval))
1105  mixed_cdft%constraint_type(:, :) = 0
1106  IF (mixed_cdft%identical_constraints) THEN
1107  DO ivar = 1, nvar
1108  mixed_cdft%constraint_type(ivar, :) = &
1109  dft_control%qs_control%cdft_control%group(ivar)%constraint_type
1110  END DO
1111  ELSE
1112  ! Possibly couple spin and charge constraints
1113  DO iforce_eval = 1, nforce_eval
1114  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1115  IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1116  qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1117  ELSE
1118  CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1119  END IF
1120  CALL get_qs_env(qs_env, dft_control=dft_control)
1121  IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1122  DO ivar = 1, nvar
1123  mixed_cdft%constraint_type(ivar, iforce_eval) = &
1124  dft_control%qs_control%cdft_control%group(ivar)%constraint_type
1125  END DO
1126  END IF
1127  END DO
1128  CALL force_env%para_env%sum(mixed_cdft%constraint_type)
1129  END IF
1130  END IF
1131  ! Transfer data from sub_force_envs to temporaries
1132  ALLOCATE (mixed_cdft%matrix%mixed_mo_coeff(nforce_eval, nspins))
1133  mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
1134  ALLOCATE (mixed_cdft%matrix%w_matrix(nforce_eval, nvar))
1135  w_matrix => mixed_cdft%matrix%w_matrix
1136  CALL dbcsr_init_p(mixed_cdft%matrix%mixed_matrix_s)
1137  mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
1138  IF (mixed_cdft%calculate_metric) THEN
1139  ALLOCATE (mixed_cdft%matrix%density_matrix(nforce_eval, nspins))
1140  density_matrix => mixed_cdft%matrix%density_matrix
1141  END IF
1142  ALLOCATE (mo_coeff_tmp(nforce_eval, nspins), wmat_tmp(nforce_eval, nvar))
1143  ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
1144  IF (mixed_cdft%calculate_metric) ALLOCATE (matrix_p_tmp(nforce_eval, nspins))
1145  IF (.NOT. uniform_occupation) THEN
1146  ALLOCATE (mixed_cdft%occupations(nforce_eval, nspins))
1147  ALLOCATE (occno_tmp(nforce_eval, nspins))
1148  END IF
1149  DO iforce_eval = 1, nforce_eval
1150  ! Temporary arrays need to be nulled on every process
1151  DO ispin = 1, nspins
1152  ! Valgrind 3.12/gfortran 4.8.4 oddly complains here (unconditional jump)
1153  ! if mixed_cdft%calculate_metric = .FALSE. and the need to null the array
1154  ! is queried with IF (mixed_cdft%calculate_metric) &
1155  IF (.NOT. uniform_occupation) &
1156  NULLIFY (occno_tmp(iforce_eval, ispin)%array)
1157  END DO
1158  IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1159  ! From this point onward, we access data local to the sub_force_envs
1160  ! Get qs_env
1161  force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
1162  IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1163  qs_env => force_env_qs%qmmm_env%qs_env
1164  ELSE
1165  CALL force_env_get(force_env_qs, qs_env=qs_env)
1166  END IF
1167  CALL get_qs_env(qs_env, dft_control=dft_control, blacs_env=blacs_env)
1168  ! Store dimensions of the transferred arrays
1169  CALL dbcsr_get_info(dft_control%qs_control%cdft_control%matrix_s%matrix, &
1170  nfullrows_total=nrow_overlap, nfullcols_total=ncol_overlap)
1171  CALL dbcsr_get_info(dft_control%qs_control%cdft_control%wmat(1)%matrix, &
1172  nfullrows_total=nrow_wmat, nfullcols_total=ncol_wmat)
1173  ! MO Coefficients
1174  DO ispin = 1, nspins
1175  CALL cp_fm_get_info(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
1176  ncol_global=ncol_mo(ispin), nrow_global=nrow_mo(ispin))
1177  CALL cp_fm_create(matrix=mo_coeff_tmp(iforce_eval, ispin), &
1178  matrix_struct=dft_control%qs_control%cdft_control%mo_coeff(ispin)%matrix_struct, &
1179  name="MO_COEFF_"//trim(adjustl(cp_to_string(iforce_eval)))//"_" &
1180  //trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1181  CALL cp_fm_to_fm(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
1182  mo_coeff_tmp(iforce_eval, ispin))
1183  END DO
1184  CALL cp_fm_release(dft_control%qs_control%cdft_control%mo_coeff)
1185  ! Matrix representation(s) of the weight function(s) (dbcsr -> fm)
1186  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_wmat, ncol_global=ncol_wmat, context=blacs_env, &
1187  para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env, &
1188  square_blocks=.true.)
1189  DO ivar = 1, nvar
1190  CALL cp_fm_create(wmat_tmp(iforce_eval, ivar), fm_struct_tmp, name="w_matrix")
1191  CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%wmat(ivar)%matrix, desymm_tmp)
1192  CALL copy_dbcsr_to_fm(desymm_tmp, wmat_tmp(iforce_eval, ivar))
1193  CALL dbcsr_release(desymm_tmp)
1194  CALL dbcsr_release_p(dft_control%qs_control%cdft_control%wmat(ivar)%matrix)
1195  END DO
1196  DEALLOCATE (dft_control%qs_control%cdft_control%wmat)
1197  CALL cp_fm_struct_release(fm_struct_tmp)
1198  ! Overlap matrix is the same for all sub_force_envs, so we just copy the first one (dbcsr -> fm)
1199  IF (iforce_eval == 1) THEN
1200  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_overlap, &
1201  ncol_global=ncol_overlap, context=blacs_env, &
1202  para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
1203  CALL cp_fm_create(matrix_s_tmp, fm_struct_tmp, name="s_matrix")
1204  CALL cp_fm_struct_release(fm_struct_tmp)
1205  CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_s%matrix, desymm_tmp)
1206  CALL copy_dbcsr_to_fm(desymm_tmp, matrix_s_tmp)
1207  CALL dbcsr_release(desymm_tmp)
1208  END IF
1209  CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_s%matrix)
1210  ! Density_matrix (dbcsr -> fm)
1211  IF (mixed_cdft%calculate_metric) THEN
1212  DO ispin = 1, nspins
1213  ! Size AOxAO
1214  CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_overlap, &
1215  ncol_global=ncol_overlap, context=blacs_env, &
1216  para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
1217  CALL cp_fm_create(matrix_p_tmp(iforce_eval, ispin), fm_struct_tmp, name="dm_matrix")
1218  CALL cp_fm_struct_release(fm_struct_tmp)
1219  CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix, desymm_tmp)
1220  CALL copy_dbcsr_to_fm(desymm_tmp, matrix_p_tmp(iforce_eval, ispin))
1221  CALL dbcsr_release(desymm_tmp)
1222  CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix)
1223  END DO
1224  DEALLOCATE (dft_control%qs_control%cdft_control%matrix_p)
1225  END IF
1226  ! Occupation numbers
1227  IF (.NOT. uniform_occupation) THEN
1228  DO ispin = 1, nspins
1229  IF (ncol_mo(ispin) /= SIZE(dft_control%qs_control%cdft_control%occupations(ispin)%array)) &
1230  cpabort("Array dimensions dont match.")
1231  IF (force_env_qs%para_env%is_source()) THEN
1232  ALLOCATE (occno_tmp(iforce_eval, ispin)%array(ncol_mo(ispin)))
1233  occno_tmp(iforce_eval, ispin)%array = dft_control%qs_control%cdft_control%occupations(ispin)%array
1234  END IF
1235  DEALLOCATE (dft_control%qs_control%cdft_control%occupations(ispin)%array)
1236  END DO
1237  DEALLOCATE (dft_control%qs_control%cdft_control%occupations)
1238  END IF
1239  END DO
1240  ! Create needed fm structs
1241  CALL cp_fm_struct_create(fm_struct_wmat, nrow_global=nrow_wmat, ncol_global=ncol_wmat, &
1242  context=mixed_cdft%blacs_env, para_env=force_env%para_env)
1243  CALL cp_fm_struct_create(fm_struct_overlap, nrow_global=nrow_overlap, ncol_global=ncol_overlap, &
1244  context=mixed_cdft%blacs_env, para_env=force_env%para_env)
1245  ! Redistribute arrays with copy_general (this is not optimal for dbcsr matrices but...)
1246  ! We use this method for the serial case (mixed_cdft%run_type == mixed_cdft_serial) as well to move the arrays to the
1247  ! correct blacs_env, which is impossible using a simple copy of the arrays
1248  ALLOCATE (mixed_wmat_tmp(nforce_eval, nvar))
1249  IF (mixed_cdft%calculate_metric) &
1250  ALLOCATE (mixed_matrix_p_tmp(nforce_eval, nspins))
1251  DO iforce_eval = 1, nforce_eval
1252  ! MO coefficients
1253  DO ispin = 1, nspins
1254  NULLIFY (fm_struct_mo)
1255  CALL cp_fm_struct_create(fm_struct_mo, nrow_global=nrow_mo(ispin), ncol_global=ncol_mo(ispin), &
1256  context=mixed_cdft%blacs_env, para_env=force_env%para_env)
1257  CALL cp_fm_create(matrix=mixed_mo_coeff(iforce_eval, ispin), &
1258  matrix_struct=fm_struct_mo, &
1259  name="MO_COEFF_"//trim(adjustl(cp_to_string(iforce_eval)))//"_" &
1260  //trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1261  CALL cp_fm_copy_general(mo_coeff_tmp(iforce_eval, ispin), &
1262  mixed_mo_coeff(iforce_eval, ispin), &
1263  mixed_cdft%blacs_env%para_env)
1264  CALL cp_fm_release(mo_coeff_tmp(iforce_eval, ispin))
1265  CALL cp_fm_struct_release(fm_struct_mo)
1266  END DO
1267  ! Weight
1268  DO ivar = 1, nvar
1269  NULLIFY (w_matrix(iforce_eval, ivar)%matrix)
1270  CALL dbcsr_init_p(w_matrix(iforce_eval, ivar)%matrix)
1271  CALL cp_fm_create(matrix=mixed_wmat_tmp(iforce_eval, ivar), &
1272  matrix_struct=fm_struct_wmat, &
1273  name="WEIGHT_"//trim(adjustl(cp_to_string(iforce_eval)))//"_MATRIX")
1274  CALL cp_fm_copy_general(wmat_tmp(iforce_eval, ivar), &
1275  mixed_wmat_tmp(iforce_eval, ivar), &
1276  mixed_cdft%blacs_env%para_env)
1277  CALL cp_fm_release(wmat_tmp(iforce_eval, ivar))
1278  ! (fm -> dbcsr)
1279  CALL copy_fm_to_dbcsr_bc(mixed_wmat_tmp(iforce_eval, ivar), &
1280  w_matrix(iforce_eval, ivar)%matrix)
1281  CALL cp_fm_release(mixed_wmat_tmp(iforce_eval, ivar))
1282  END DO
1283  ! Density matrix (fm -> dbcsr)
1284  IF (mixed_cdft%calculate_metric) THEN
1285  DO ispin = 1, nspins
1286  NULLIFY (density_matrix(iforce_eval, ispin)%matrix)
1287  CALL dbcsr_init_p(density_matrix(iforce_eval, ispin)%matrix)
1288  CALL cp_fm_create(matrix=mixed_matrix_p_tmp(iforce_eval, ispin), &
1289  matrix_struct=fm_struct_overlap, &
1290  name="DENSITY_"//trim(adjustl(cp_to_string(iforce_eval)))//"_"// &
1291  trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1292  CALL cp_fm_copy_general(matrix_p_tmp(iforce_eval, ispin), &
1293  mixed_matrix_p_tmp(iforce_eval, ispin), &
1294  mixed_cdft%blacs_env%para_env)
1295  CALL cp_fm_release(matrix_p_tmp(iforce_eval, ispin))
1296  CALL copy_fm_to_dbcsr_bc(mixed_matrix_p_tmp(iforce_eval, ispin), &
1297  density_matrix(iforce_eval, ispin)%matrix)
1298  CALL cp_fm_release(mixed_matrix_p_tmp(iforce_eval, ispin))
1299  END DO
1300  END IF
1301  END DO
1302  CALL cp_fm_struct_release(fm_struct_wmat)
1303  DEALLOCATE (mo_coeff_tmp, wmat_tmp, mixed_wmat_tmp)
1304  IF (mixed_cdft%calculate_metric) THEN
1305  DEALLOCATE (matrix_p_tmp)
1306  DEALLOCATE (mixed_matrix_p_tmp)
1307  END IF
1308  ! Overlap (fm -> dbcsr)
1309  CALL cp_fm_create(matrix=mixed_matrix_s_tmp, &
1310  matrix_struct=fm_struct_overlap, &
1311  name="OVERLAP_MATRIX")
1312  CALL cp_fm_struct_release(fm_struct_overlap)
1313  CALL cp_fm_copy_general(matrix_s_tmp, &
1314  mixed_matrix_s_tmp, &
1315  mixed_cdft%blacs_env%para_env)
1316  CALL cp_fm_release(matrix_s_tmp)
1317  CALL copy_fm_to_dbcsr_bc(mixed_matrix_s_tmp, mixed_matrix_s)
1318  CALL cp_fm_release(mixed_matrix_s_tmp)
1319  ! Occupation numbers
1320  IF (.NOT. uniform_occupation) THEN
1321  DO iforce_eval = 1, nforce_eval
1322  DO ispin = 1, nspins
1323  ALLOCATE (mixed_cdft%occupations(iforce_eval, ispin)%array(ncol_mo(ispin)))
1324  mixed_cdft%occupations(iforce_eval, ispin)%array = 0.0_dp
1325  IF (ASSOCIATED(occno_tmp(iforce_eval, ispin)%array)) THEN
1326  mixed_cdft%occupations(iforce_eval, ispin)%array = occno_tmp(iforce_eval, ispin)%array
1327  DEALLOCATE (occno_tmp(iforce_eval, ispin)%array)
1328  END IF
1329  CALL force_env%para_env%sum(mixed_cdft%occupations(iforce_eval, ispin)%array)
1330  END DO
1331  END DO
1332  DEALLOCATE (occno_tmp)
1333  END IF
1334  DEALLOCATE (ncol_mo, nrow_mo)
1335 
1336  END SUBROUTINE mixed_cdft_redistribute_arrays
1337 ! **************************************************************************************************
1338 !> \brief Routine to print out the electronic coupling(s) between CDFT states.
1339 !> \param force_env the force_env that holds the CDFT states
1340 !> \par History
1341 !> 11.17 created [Nico Holmberg]
1342 ! **************************************************************************************************
1343  SUBROUTINE mixed_cdft_print_couplings(force_env)
1344  TYPE(force_env_type), POINTER :: force_env
1345 
1346  INTEGER :: iounit, ipermutation, istate, ivar, &
1347  jstate, nforce_eval, npermutations, &
1348  nvar
1349  TYPE(cp_logger_type), POINTER :: logger
1350  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1351  TYPE(section_vals_type), POINTER :: force_env_section, print_section
1352 
1353  NULLIFY (print_section, mixed_cdft)
1354 
1355  logger => cp_get_default_logger()
1356  cpassert(ASSOCIATED(force_env))
1357  CALL force_env_get(force_env=force_env, &
1358  force_env_section=force_env_section)
1359  CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1360  cpassert(ASSOCIATED(mixed_cdft))
1361  print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1362  iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1363  !
1364  cpassert(ALLOCATED(mixed_cdft%results%strength))
1365  cpassert(ALLOCATED(mixed_cdft%results%W_diagonal))
1366  cpassert(ALLOCATED(mixed_cdft%results%S))
1367  cpassert(ALLOCATED(mixed_cdft%results%energy))
1368  nforce_eval = SIZE(force_env%sub_force_env)
1369  nvar = SIZE(mixed_cdft%results%strength, 1)
1370  npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
1371  IF (iounit > 0) THEN
1372  WRITE (iounit, '(/,T3,A,T66)') &
1373  '------------------------- CDFT coupling information --------------------------'
1374  WRITE (iounit, '(T3,A,T66,(3X,F12.2))') &
1375  'Information at step (fs):', mixed_cdft%sim_step*mixed_cdft%sim_dt
1376  DO ipermutation = 1, npermutations
1377  CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1378  WRITE (iounit, '(/,T3,A)') repeat('#', 44)
1379  WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### CDFT states I =', istate, ' and J = ', jstate, ' ######'
1380  WRITE (iounit, '(T3,A)') repeat('#', 44)
1381  DO ivar = 1, nvar
1382  IF (ivar > 1) &
1383  WRITE (iounit, '(A)') ''
1384  WRITE (iounit, '(T3,A,T60,(3X,I18))') 'Atomic group:', ivar
1385  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1386  'Strength of constraint I:', mixed_cdft%results%strength(ivar, istate)
1387  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1388  'Strength of constraint J:', mixed_cdft%results%strength(ivar, jstate)
1389  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1390  'Final value of constraint I:', mixed_cdft%results%W_diagonal(ivar, istate)
1391  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1392  'Final value of constraint J:', mixed_cdft%results%W_diagonal(ivar, jstate)
1393  END DO
1394  WRITE (iounit, '(/,T3,A,T60,(3X,F18.12))') &
1395  'Overlap between states I and J:', mixed_cdft%results%S(istate, jstate)
1396  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1397  'Charge transfer energy (J-I) (Hartree):', (mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate))
1398  WRITE (iounit, *)
1399  IF (ALLOCATED(mixed_cdft%results%rotation)) THEN
1400  IF (abs(mixed_cdft%results%rotation(ipermutation))*1.0e3_dp .GE. 0.1_dp) THEN
1401  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1402  'Diabatic electronic coupling (rotation, mHartree):', &
1403  abs(mixed_cdft%results%rotation(ipermutation)*1.0e3_dp)
1404  ELSE
1405  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1406  'Diabatic electronic coupling (rotation, microHartree):', &
1407  abs(mixed_cdft%results%rotation(ipermutation)*1.0e6_dp)
1408  END IF
1409  END IF
1410  IF (ALLOCATED(mixed_cdft%results%lowdin)) THEN
1411  IF (abs(mixed_cdft%results%lowdin(ipermutation))*1.0e3_dp .GE. 0.1_dp) THEN
1412  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1413  'Diabatic electronic coupling (Lowdin, mHartree):', &
1414  abs(mixed_cdft%results%lowdin(ipermutation)*1.0e3_dp)
1415  ELSE
1416  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1417  'Diabatic electronic coupling (Lowdin, microHartree):', &
1418  abs(mixed_cdft%results%lowdin(ipermutation)*1.0e6_dp)
1419  END IF
1420  END IF
1421  IF (ALLOCATED(mixed_cdft%results%wfn)) THEN
1422  IF (mixed_cdft%results%wfn(ipermutation)*1.0e3_dp .GE. 0.1_dp) THEN
1423  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1424  'Diabatic electronic coupling (wfn overlap, mHartree):', &
1425  abs(mixed_cdft%results%wfn(ipermutation)*1.0e3_dp)
1426  ELSE
1427  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1428  'Diabatic electronic coupling (wfn overlap, microHartree):', &
1429  abs(mixed_cdft%results%wfn(ipermutation)*1.0e6_dp)
1430  END IF
1431  END IF
1432  IF (ALLOCATED(mixed_cdft%results%nonortho)) THEN
1433  WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1434  'Diabatic electronic coupling (nonorthogonal, Hartree):', mixed_cdft%results%nonortho(ipermutation)
1435  END IF
1436  IF (ALLOCATED(mixed_cdft%results%metric)) THEN
1437  WRITE (iounit, *)
1438  IF (SIZE(mixed_cdft%results%metric, 2) == 1) THEN
1439  WRITE (iounit, '(T3,A,T66,(3X,F12.6))') &
1440  'Coupling reliability metric (0 is ideal):', mixed_cdft%results%metric(ipermutation, 1)
1441  ELSE
1442  WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
1443  'Coupling reliability metric (0 is ideal):', &
1444  mixed_cdft%results%metric(ipermutation, 1), mixed_cdft%results%metric(ipermutation, 2)
1445  END IF
1446  END IF
1447  END DO
1448  WRITE (iounit, '(T3,A)') &
1449  '------------------------------------------------------------------------------'
1450  END IF
1451  CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1452  "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1453 
1454  END SUBROUTINE mixed_cdft_print_couplings
1455 
1456 ! **************************************************************************************************
1457 !> \brief Release storage reserved for mixed CDFT matrices
1458 !> \param force_env the force_env that holds the CDFT states
1459 !> \par History
1460 !> 11.17 created [Nico Holmberg]
1461 ! **************************************************************************************************
1462  SUBROUTINE mixed_cdft_release_work(force_env)
1463  TYPE(force_env_type), POINTER :: force_env
1464 
1465  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1466 
1467  NULLIFY (mixed_cdft)
1468 
1469  cpassert(ASSOCIATED(force_env))
1470  CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1471  cpassert(ASSOCIATED(mixed_cdft))
1472  CALL mixed_cdft_result_type_release(mixed_cdft%results)
1473 
1474  END SUBROUTINE mixed_cdft_release_work
1475 
1476 ! **************************************************************************************************
1477 !> \brief Given the size of a symmetric matrix and a permutation index, returns indices (i, j) of the
1478 !> off-diagonal element that corresponds to the permutation index. Assumes that the permutation
1479 !> index was computed by going through the upper triangular part of the input matrix row-by-row.
1480 !> \param n the size of the symmetric matrix
1481 !> \param ipermutation the permutation index
1482 !> \param i the row index corresponding to ipermutation
1483 !> \param j the column index corresponding to ipermutation
1484 ! **************************************************************************************************
1485  SUBROUTINE map_permutation_to_states(n, ipermutation, i, j)
1486  INTEGER, INTENT(IN) :: n, ipermutation
1487  INTEGER, INTENT(OUT) :: i, j
1488 
1489  INTEGER :: kcol, kpermutation, krow, npermutations
1490 
1491  npermutations = n*(n - 1)/2 ! Size of upper triangular part
1492  IF (ipermutation > npermutations) &
1493  cpabort("Permutation index out of bounds")
1494  kpermutation = 0
1495  DO krow = 1, n
1496  DO kcol = krow + 1, n
1497  kpermutation = kpermutation + 1
1498  IF (kpermutation == ipermutation) THEN
1499  i = krow
1500  j = kcol
1501  RETURN
1502  END IF
1503  END DO
1504  END DO
1505 
1506  END SUBROUTINE map_permutation_to_states
1507 
1508 ! **************************************************************************************************
1509 !> \brief Determine confinement bounds along confinement dir (hardcoded to be z)
1510 !> and determine the number of nonzero entries
1511 !> Optionally zero entries below a given threshold
1512 !> \param fun input 3D potential (real space)
1513 !> \param th threshold for screening values
1514 !> \param just_zero determines if fun should only be zeroed without returning bounds/work
1515 !> \param bounds the confinement bounds: fun is nonzero only between these values along 3rd dimension
1516 !> \param work an estimate of the total number of grid points where fun is nonzero
1517 ! **************************************************************************************************
1518  SUBROUTINE hfun_zero(fun, th, just_zero, bounds, work)
1519  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: fun
1520  REAL(kind=dp), INTENT(IN) :: th
1521  LOGICAL :: just_zero
1522  INTEGER, OPTIONAL :: bounds(2), work
1523 
1524  INTEGER :: i1, i2, i3, lb, n1, n2, n3, nzeroed, &
1525  nzeroed_total, ub
1526  LOGICAL :: lb_final, ub_final
1527 
1528  n1 = SIZE(fun, 1)
1529  n2 = SIZE(fun, 2)
1530  n3 = SIZE(fun, 3)
1531  nzeroed_total = 0
1532  IF (.NOT. just_zero) THEN
1533  cpassert(PRESENT(bounds))
1534  cpassert(PRESENT(work))
1535  lb = 1
1536  lb_final = .false.
1537  ub_final = .false.
1538  END IF
1539  DO i3 = 1, n3
1540  IF (.NOT. just_zero) nzeroed = 0
1541  DO i2 = 1, n2
1542  DO i1 = 1, n1
1543  IF (fun(i1, i2, i3) < th) THEN
1544  IF (.NOT. just_zero) THEN
1545  nzeroed = nzeroed + 1
1546  nzeroed_total = nzeroed_total + 1
1547  ELSE
1548  fun(i1, i2, i3) = 0.0_dp
1549  END IF
1550  END IF
1551  END DO
1552  END DO
1553  IF (.NOT. just_zero) THEN
1554  IF (nzeroed == (n2*n1)) THEN
1555  IF (.NOT. lb_final) THEN
1556  lb = i3
1557  ELSE IF (.NOT. ub_final) THEN
1558  ub = i3
1559  ub_final = .true.
1560  END IF
1561  ELSE
1562  IF (.NOT. lb_final) lb_final = .true.
1563  IF (ub_final) ub_final = .false. ! Safeguard against "holes"
1564  END IF
1565  END IF
1566  END DO
1567  IF (.NOT. just_zero) THEN
1568  IF (.NOT. ub_final) ub = n3
1569  bounds(1) = lb
1570  bounds(2) = ub
1571  bounds = bounds - (n3/2) - 1
1572  work = n3*n2*n1 - nzeroed_total
1573  END IF
1574 
1575  END SUBROUTINE hfun_zero
1576 
1577 ! **************************************************************************************************
1578 !> \brief Read input section related to block diagonalization of the mixed CDFT Hamiltonian matrix.
1579 !> \param force_env the force_env that holds the CDFT states
1580 !> \param blocks list of CDFT states defining the matrix blocks
1581 !> \param ignore_excited flag that determines if excited states resulting from the block
1582 !> diagonalization process should be ignored
1583 !> \param nrecursion integer that determines how many steps of recursive block diagonalization
1584 !> is performed (1 if disabled)
1585 !> \par History
1586 !> 01.18 created [Nico Holmberg]
1587 ! **************************************************************************************************
1588  SUBROUTINE mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
1589  TYPE(force_env_type), POINTER :: force_env
1590  TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:), &
1591  INTENT(OUT) :: blocks
1592  LOGICAL, INTENT(OUT) :: ignore_excited
1593  INTEGER, INTENT(OUT) :: nrecursion
1594 
1595  INTEGER :: i, j, k, l, nblk, nforce_eval
1596  INTEGER, DIMENSION(:), POINTER :: tmplist
1597  LOGICAL :: do_recursive, explicit, has_duplicates
1598  TYPE(section_vals_type), POINTER :: block_section, force_env_section
1599 
1600  EXTERNAL :: dsygv
1601 
1602  NULLIFY (force_env_section, block_section)
1603  cpassert(ASSOCIATED(force_env))
1604  nforce_eval = SIZE(force_env%sub_force_env)
1605 
1606  CALL force_env_get(force_env=force_env, &
1607  force_env_section=force_env_section)
1608  block_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%BLOCK_DIAGONALIZE")
1609 
1610  CALL section_vals_get(block_section, explicit=explicit)
1611  IF (.NOT. explicit) &
1612  CALL cp_abort(__location__, &
1613  "Block diagonalization of CDFT Hamiltonian was requested, but the "// &
1614  "corresponding input section is missing!")
1615 
1616  CALL section_vals_val_get(block_section, "BLOCK", n_rep_val=nblk)
1617  ALLOCATE (blocks(nblk))
1618  DO i = 1, nblk
1619  NULLIFY (blocks(i)%array)
1620  CALL section_vals_val_get(block_section, "BLOCK", i_rep_val=i, i_vals=tmplist)
1621  IF (SIZE(tmplist) < 1) &
1622  cpabort("Each BLOCK must contain at least 1 state.")
1623  ALLOCATE (blocks(i)%array(SIZE(tmplist)))
1624  blocks(i)%array(:) = tmplist(:)
1625  END DO
1626  CALL section_vals_val_get(block_section, "IGNORE_EXCITED", l_val=ignore_excited)
1627  CALL section_vals_val_get(block_section, "RECURSIVE_DIAGONALIZATION", l_val=do_recursive)
1628  ! Check that the requested states exist
1629  DO i = 1, nblk
1630  DO j = 1, SIZE(blocks(i)%array)
1631  IF (blocks(i)%array(j) < 1 .OR. blocks(i)%array(j) > nforce_eval) &
1632  cpabort("Requested state does not exist.")
1633  END DO
1634  END DO
1635  ! Check for duplicates
1636  has_duplicates = .false.
1637  DO i = 1, nblk
1638  ! Within same block
1639  DO j = 1, SIZE(blocks(i)%array)
1640  DO k = j + 1, SIZE(blocks(i)%array)
1641  IF (blocks(i)%array(j) == blocks(i)%array(k)) has_duplicates = .true.
1642  END DO
1643  END DO
1644  ! Within different blocks
1645  DO j = i + 1, nblk
1646  DO k = 1, SIZE(blocks(i)%array)
1647  DO l = 1, SIZE(blocks(j)%array)
1648  IF (blocks(i)%array(k) == blocks(j)%array(l)) has_duplicates = .true.
1649  END DO
1650  END DO
1651  END DO
1652  END DO
1653  IF (has_duplicates) cpabort("Duplicate states are not allowed.")
1654  nrecursion = 1
1655  IF (do_recursive) THEN
1656  IF (modulo(nblk, 2) /= 0) THEN
1657  CALL cp_warn(__location__, &
1658  "Number of blocks not divisible with 2. Recursive diagonalization not possible. "// &
1659  "Calculation proceeds without.")
1660  nrecursion = 1
1661  ELSE
1662  nrecursion = nblk/2
1663  END IF
1664  IF (nrecursion /= 1 .AND. .NOT. ignore_excited) &
1665  CALL cp_abort(__location__, &
1666  "Keyword IGNORE_EXCITED must be active for recursive diagonalization.")
1667  END IF
1668 
1669  END SUBROUTINE mixed_cdft_read_block_diag
1670 
1671 ! **************************************************************************************************
1672 !> \brief Assembles the matrix blocks from the mixed CDFT Hamiltonian.
1673 !> \param mixed_cdft the env that holds the CDFT states
1674 !> \param blocks list of CDFT states defining the matrix blocks
1675 !> \param H_block list of Hamiltonian matrix blocks
1676 !> \param S_block list of overlap matrix blocks
1677 !> \par History
1678 !> 01.18 created [Nico Holmberg]
1679 ! **************************************************************************************************
1680  SUBROUTINE mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block)
1681  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1682  TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1683  TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:), &
1684  INTENT(OUT) :: h_block, s_block
1685 
1686  INTEGER :: i, icol, irow, j, k, nblk
1687 
1688  EXTERNAL :: dsygv
1689 
1690  cpassert(ASSOCIATED(mixed_cdft))
1691 
1692  nblk = SIZE(blocks)
1693  ALLOCATE (h_block(nblk), s_block(nblk))
1694  DO i = 1, nblk
1695  NULLIFY (h_block(i)%array)
1696  NULLIFY (s_block(i)%array)
1697  ALLOCATE (h_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1698  ALLOCATE (s_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1699  icol = 0
1700  DO j = 1, SIZE(blocks(i)%array)
1701  irow = 0
1702  icol = icol + 1
1703  DO k = 1, SIZE(blocks(i)%array)
1704  irow = irow + 1
1705  h_block(i)%array(irow, icol) = mixed_cdft%results%H(blocks(i)%array(k), blocks(i)%array(j))
1706  s_block(i)%array(irow, icol) = mixed_cdft%results%S(blocks(i)%array(k), blocks(i)%array(j))
1707  END DO
1708  END DO
1709  ! Check that none of the interaction energies is repulsive
1710  IF (any(h_block(i)%array .GE. 0.0_dp)) &
1711  CALL cp_abort(__location__, &
1712  "At least one of the interaction energies within block "//trim(adjustl(cp_to_string(i)))// &
1713  " is repulsive.")
1714  END DO
1715 
1716  END SUBROUTINE mixed_cdft_get_blocks
1717 
1718 ! **************************************************************************************************
1719 !> \brief Diagonalizes each of the matrix blocks.
1720 !> \param blocks list of CDFT states defining the matrix blocks
1721 !> \param H_block list of Hamiltonian matrix blocks
1722 !> \param S_block list of overlap matrix blocks
1723 !> \param eigenvalues list of eigenvalues for each block
1724 !> \par History
1725 !> 01.18 created [Nico Holmberg]
1726 ! **************************************************************************************************
1727  SUBROUTINE mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues)
1728  TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1729  TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: h_block, s_block
1730  TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:), &
1731  INTENT(OUT) :: eigenvalues
1732 
1733  INTEGER :: i, info, nblk, work_array_size
1734  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
1735  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: h_mat_copy, s_mat_copy
1736 
1737  EXTERNAL :: dsygv
1738 
1739  nblk = SIZE(blocks)
1740  ALLOCATE (eigenvalues(nblk))
1741  DO i = 1, nblk
1742  NULLIFY (eigenvalues(i)%array)
1743  ALLOCATE (eigenvalues(i)%array(SIZE(blocks(i)%array)))
1744  eigenvalues(i)%array = 0.0_dp
1745  ! Workspace query
1746  ALLOCATE (work(1))
1747  info = 0
1748  ALLOCATE (h_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1749  ALLOCATE (s_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1750  h_mat_copy(:, :) = h_block(i)%array(:, :) ! Need explicit copies because dsygv destroys original values
1751  s_mat_copy(:, :) = s_block(i)%array(:, :)
1752  CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), h_mat_copy, SIZE(blocks(i)%array), &
1753  s_mat_copy, SIZE(blocks(i)%array), eigenvalues(i)%array, work, -1, info)
1754  work_array_size = nint(work(1))
1755  DEALLOCATE (h_mat_copy, s_mat_copy)
1756  ! Allocate work array
1757  DEALLOCATE (work)
1758  ALLOCATE (work(work_array_size))
1759  work = 0.0_dp
1760  ! Solve Hc = eSc
1761  info = 0
1762  CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), h_block(i)%array, SIZE(blocks(i)%array), &
1763  s_block(i)%array, SIZE(blocks(i)%array), eigenvalues(i)%array, work, work_array_size, info)
1764  IF (info /= 0) THEN
1765  IF (info > SIZE(blocks(i)%array)) THEN
1766  cpabort("Matrix S is not positive definite")
1767  ELSE
1768  cpabort("Diagonalization of H matrix failed.")
1769  END IF
1770  END IF
1771  DEALLOCATE (work)
1772  END DO
1773 
1774  END SUBROUTINE mixed_cdft_diagonalize_blocks
1775 
1776 ! **************************************************************************************************
1777 !> \brief Assembles the new block diagonalized mixed CDFT Hamiltonian and overlap matrices.
1778 !> \param mixed_cdft the env that holds the CDFT states
1779 !> \param blocks list of CDFT states defining the matrix blocks
1780 !> \param H_block list of Hamiltonian matrix blocks
1781 !> \param eigenvalues list of eigenvalues for each block
1782 !> \param n size of the new Hamiltonian and overlap matrices
1783 !> \param iounit the output unit
1784 !> \par History
1785 !> 01.18 created [Nico Holmberg]
1786 ! **************************************************************************************************
1787  SUBROUTINE mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, &
1788  n, iounit)
1789  TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1790  TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1791  TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: h_block
1792  TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1793  INTEGER :: n, iounit
1794 
1795  CHARACTER(LEN=20) :: ilabel, jlabel
1796  CHARACTER(LEN=3) :: tmp
1797  INTEGER :: i, icol, ipermutation, irow, j, k, l, &
1798  nblk, npermutations
1799  LOGICAL :: ignore_excited
1800  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: h_mat, h_offdiag, s_mat, s_offdiag
1801 
1802  EXTERNAL :: dsygv
1803 
1804  ALLOCATE (h_mat(n, n), s_mat(n, n))
1805  nblk = SIZE(blocks)
1806  ignore_excited = (nblk == n)
1807  ! The diagonal contains the eigenvalues of each block
1808  IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Eigenvalues of the block diagonalized states"
1809  h_mat(:, :) = 0.0_dp
1810  s_mat(:, :) = 0.0_dp
1811  k = 1
1812  DO i = 1, nblk
1813  IF (iounit > 0) WRITE (iounit, '(T6,A,I3)') "Block", i
1814  DO j = 1, SIZE(eigenvalues(i)%array)
1815  h_mat(k, k) = eigenvalues(i)%array(j)
1816  s_mat(k, k) = 1.0_dp
1817  k = k + 1
1818  IF (iounit > 0) THEN
1819  IF (j == 1) THEN
1820  WRITE (iounit, '(T9,A,T58,(3X,F20.14))') 'Ground state energy:', eigenvalues(i)%array(j)
1821  ELSE
1822  WRITE (iounit, '(T9,A,I2,A,T58,(3X,F20.14))') &
1823  'Excited state (', j - 1, ' ) energy:', eigenvalues(i)%array(j)
1824  END IF
1825  END IF
1826  IF (ignore_excited .AND. j == 1) EXIT
1827  END DO
1828  END DO
1829  ! Transform the off-diagonal blocks using the eigenvectors of each block
1830  npermutations = nblk*(nblk - 1)/2
1831  IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Interactions between block diagonalized states"
1832  DO ipermutation = 1, npermutations
1833  CALL map_permutation_to_states(nblk, ipermutation, i, j)
1834  ! Get the untransformed off-diagonal block
1835  ALLOCATE (h_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
1836  ALLOCATE (s_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
1837  icol = 0
1838  DO k = 1, SIZE(blocks(j)%array)
1839  irow = 0
1840  icol = icol + 1
1841  DO l = 1, SIZE(blocks(i)%array)
1842  irow = irow + 1
1843  h_offdiag(irow, icol) = mixed_cdft%results%H(blocks(i)%array(l), blocks(j)%array(k))
1844  s_offdiag(irow, icol) = mixed_cdft%results%S(blocks(i)%array(l), blocks(j)%array(k))
1845  END DO
1846  END DO
1847  ! Check that none of the interaction energies is repulsive
1848  IF (any(h_offdiag .GE. 0.0_dp)) &
1849  CALL cp_abort(__location__, &
1850  "At least one of the interaction energies between blocks "//trim(adjustl(cp_to_string(i)))// &
1851  " and "//trim(adjustl(cp_to_string(j)))//" is repulsive.")
1852  ! Now transform: C_i^T * H * C_j
1853  h_offdiag(:, :) = matmul(h_offdiag, h_block(j)%array)
1854  h_offdiag(:, :) = matmul(transpose(h_block(i)%array), h_offdiag)
1855  s_offdiag(:, :) = matmul(s_offdiag, h_block(j)%array)
1856  s_offdiag(:, :) = matmul(transpose(h_block(i)%array), s_offdiag)
1857  ! Make sure the transformation preserves the sign of elements in the S and H matrices
1858  ! The S/H matrices contain only positive/negative values so that any sign flipping occurs in the
1859  ! same elements in both matrices
1860  ! Check for sign flipping using the S matrix
1861  IF (any(s_offdiag .LT. 0.0_dp)) THEN
1862  DO l = 1, SIZE(s_offdiag, 2)
1863  DO k = 1, SIZE(s_offdiag, 1)
1864  IF (s_offdiag(k, l) .LT. 0.0_dp) THEN
1865  s_offdiag(k, l) = -1.0_dp*s_offdiag(k, l)
1866  h_offdiag(k, l) = -1.0_dp*h_offdiag(k, l)
1867  END IF
1868  END DO
1869  END DO
1870  END IF
1871  IF (ignore_excited) THEN
1872  h_mat(i, j) = h_offdiag(1, 1)
1873  h_mat(j, i) = h_mat(i, j)
1874  s_mat(i, j) = s_offdiag(1, 1)
1875  s_mat(j, i) = s_mat(i, j)
1876  ELSE
1877  irow = 1
1878  icol = 1
1879  DO k = 1, i - 1
1880  irow = irow + SIZE(blocks(k)%array)
1881  END DO
1882  DO k = 1, j - 1
1883  icol = icol + SIZE(blocks(k)%array)
1884  END DO
1885  h_mat(irow:irow + SIZE(h_offdiag, 1) - 1, icol:icol + SIZE(h_offdiag, 2) - 1) = h_offdiag(:, :)
1886  h_mat(icol:icol + SIZE(h_offdiag, 2) - 1, irow:irow + SIZE(h_offdiag, 1) - 1) = transpose(h_offdiag)
1887  s_mat(irow:irow + SIZE(h_offdiag, 1) - 1, icol:icol + SIZE(h_offdiag, 2) - 1) = s_offdiag(:, :)
1888  s_mat(icol:icol + SIZE(h_offdiag, 2) - 1, irow:irow + SIZE(h_offdiag, 1) - 1) = transpose(s_offdiag)
1889  END IF
1890  IF (iounit > 0) THEN
1891  WRITE (iounit, '(/,T3,A)') repeat('#', 39)
1892  WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### Blocks I =', i, ' and J = ', j, ' ######'
1893  WRITE (iounit, '(T3,A)') repeat('#', 39)
1894  WRITE (iounit, '(T3,A)') 'Interaction energies'
1895  DO irow = 1, SIZE(h_offdiag, 1)
1896  ilabel = "(ground state)"
1897  IF (irow > 1) THEN
1898  IF (ignore_excited) EXIT
1899  WRITE (tmp, '(I3)') irow - 1
1900  ilabel = "(excited state "//trim(adjustl(tmp))//")"
1901  END IF
1902  DO icol = 1, SIZE(h_offdiag, 2)
1903  jlabel = "(ground state)"
1904  IF (icol > 1) THEN
1905  IF (ignore_excited) EXIT
1906  WRITE (tmp, '(I3)') icol - 1
1907  jlabel = "(excited state "//trim(adjustl(tmp))//")"
1908  END IF
1909  WRITE (iounit, '(T6,A,T58,(3X,F20.14))') trim(ilabel)//'-'//trim(jlabel)//':', h_offdiag(irow, icol)
1910  END DO
1911  END DO
1912  WRITE (iounit, '(T3,A)') 'Overlaps'
1913  DO irow = 1, SIZE(h_offdiag, 1)
1914  ilabel = "(ground state)"
1915  IF (irow > 1) THEN
1916  IF (ignore_excited) EXIT
1917  ilabel = "(excited state)"
1918  WRITE (tmp, '(I3)') irow - 1
1919  ilabel = "(excited state "//trim(adjustl(tmp))//")"
1920  END IF
1921  DO icol = 1, SIZE(h_offdiag, 2)
1922  jlabel = "(ground state)"
1923  IF (icol > 1) THEN
1924  IF (ignore_excited) EXIT
1925  WRITE (tmp, '(I3)') icol - 1
1926  jlabel = "(excited state "//trim(adjustl(tmp))//")"
1927  END IF
1928  WRITE (iounit, '(T6,A,T58,(3X,F20.14))') trim(ilabel)//'-'//trim(jlabel)//':', s_offdiag(irow, icol)
1929  END DO
1930  END DO
1931  END IF
1932  DEALLOCATE (h_offdiag, s_offdiag)
1933  END DO
1934  CALL mixed_cdft_result_type_set(mixed_cdft%results, h=h_mat, s=s_mat)
1935  ! Deallocate work
1936  DEALLOCATE (h_mat, s_mat)
1937 
1938  END SUBROUTINE mixed_cdft_assemble_block_diag
1939 
1940 END MODULE mixed_cdft_utils
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
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::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::copy_fm_to_dbcsr_bc
subroutine, public copy_fm_to_dbcsr_bc(fm, bc_mat)
Copy a BLACS matrix to a dbcsr matrix with a special block-cyclic distribution, which requires no com...
Definition: cp_dbcsr_operations.F:143
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_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_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
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_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_set
subroutine, public cp_logger_set(logger, local_filename, global_filename)
sets various attributes of the given logger
Definition: cp_log_handling.F:717
cp_log_handling::cp_logger_create
subroutine, public cp_logger_create(logger, para_env, print_level, default_global_unit_nr, default_local_unit_nr, global_filename, local_filename, close_global_unit_on_dealloc, iter_info, close_local_unit_on_dealloc, suffix, template_logger)
initializes a logger
Definition: cp_log_handling.F:285
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::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::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_realspace_grid_init
Definition: cp_realspace_grid_init.F:15
cp_realspace_grid_init::init_input_type
subroutine, public init_input_type(input_settings, nsmax, rs_grid_section, ilevel, higher_grid_layout)
parses an input section to assign the proper values to the input type
Definition: cp_realspace_grid_init.F:48
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, 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)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition: cube_utils.F:18
cube_utils::return_cube_max_iradius
integer function, public return_cube_max_iradius(info)
...
Definition: cube_utils.F:175
cube_utils::init_cube_info
subroutine, public init_cube_info(info, dr, dh, dh_inv, ortho, max_radius)
...
Definition: cube_utils.F:212
d3_poly
Routines to efficiently handle dense polynomial in 3 variables up to a given degree....
Definition: d3_poly.F:23
d3_poly::init_d3_poly_module
subroutine, public init_d3_poly_module()
initialization of the cache, is called by functions in this module that use cached values
Definition: d3_poly.F:74
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
force_env_types::multiple_fe_list
subroutine, public multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
returns the order of the multiple force_env
Definition: force_env_types.F:685
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Definition: force_env_types.F:329
gaussian_gridlevels
Definition: gaussian_gridlevels.F:13
gaussian_gridlevels::init_gaussian_gridlevel
subroutine, public init_gaussian_gridlevel(gridlevel_info, ngrid_levels, cutoff, rel_cutoff, print_section)
...
Definition: gaussian_gridlevels.F:59
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
hirshfeld_types
The types needed for the calculation of Hirshfeld charges and related functions.
Definition: hirshfeld_types.F:15
hirshfeld_types::create_hirshfeld_type
subroutine, public create_hirshfeld_type(hirshfeld_env)
...
Definition: hirshfeld_types.F:66
hirshfeld_types::set_hirshfeld_info
subroutine, public set_hirshfeld_info(hirshfeld_env, shape_function_type, iterative, ref_charge, fnorm, radius_type, use_bohr)
Set values of a Hirshfeld env.
Definition: hirshfeld_types.F:181
hirshfeld_types::release_hirshfeld_type
subroutine, public release_hirshfeld_type(hirshfeld_env)
...
Definition: hirshfeld_types.F:89
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition: input_constants.F:772
input_constants::becke_cutoff_element
integer, parameter, public becke_cutoff_element
Definition: input_constants.F:758
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition: input_constants.F:718
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition: input_constants.F:718
input_constants::mixed_cdft_parallel
integer, parameter, public mixed_cdft_parallel
Definition: input_constants.F:772
input_constants::shape_function_gaussian
integer, parameter, public shape_function_gaussian
Definition: input_constants.F:626
input_constants::mixed_cdft_parallel_nobuild
integer, parameter, public mixed_cdft_parallel_nobuild
Definition: input_constants.F:772
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_duplicate
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
Definition: input_section_types.F:1732
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
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition: input_section_types.F:646
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mixed_cdft_types
Types for mixed CDFT calculations.
Definition: mixed_cdft_types.F:14
mixed_cdft_types::mixed_cdft_result_type_release
subroutine, public mixed_cdft_result_type_release(results)
Releases all arrays within the mixed CDFT result container.
Definition: mixed_cdft_types.F:563
mixed_cdft_types::mixed_cdft_work_type_init
subroutine, public mixed_cdft_work_type_init(matrix)
Initializes the mixed_cdft_work_type.
Definition: mixed_cdft_types.F:587
mixed_cdft_types::mixed_cdft_result_type_set
subroutine, public mixed_cdft_result_type_set(results, lowdin, wfn, nonortho, metric, rotation, H, S, Wad, Wda, W_diagonal, energy, strength, S_minushalf)
Updates arrays within the mixed CDFT result container.
Definition: mixed_cdft_types.F:481
mixed_cdft_utils
Utility subroutines for mixed CDFT calculations.
Definition: mixed_cdft_utils.F:14
mixed_cdft_utils::map_permutation_to_states
subroutine, public map_permutation_to_states(n, ipermutation, i, j)
Given the size of a symmetric matrix and a permutation index, returns indices (i, j) of the off-diago...
Definition: mixed_cdft_utils.F:1486
mixed_cdft_utils::mixed_cdft_assemble_block_diag
subroutine, public mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, n, iounit)
Assembles the new block diagonalized mixed CDFT Hamiltonian and overlap matrices.
Definition: mixed_cdft_utils.F:1789
mixed_cdft_utils::mixed_cdft_diagonalize_blocks
subroutine, public mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues)
Diagonalizes each of the matrix blocks.
Definition: mixed_cdft_utils.F:1728
mixed_cdft_utils::mixed_cdft_init_structures
subroutine, public mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
Initialize all the structures needed for a mixed CDFT calculation.
Definition: mixed_cdft_utils.F:632
mixed_cdft_utils::mixed_cdft_transfer_settings
subroutine, public mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
Transfer settings to mixed_cdft.
Definition: mixed_cdft_utils.F:534
mixed_cdft_utils::mixed_cdft_print_couplings
subroutine, public mixed_cdft_print_couplings(force_env)
Routine to print out the electronic coupling(s) between CDFT states.
Definition: mixed_cdft_utils.F:1344
mixed_cdft_utils::mixed_cdft_get_blocks
subroutine, public mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block)
Assembles the matrix blocks from the mixed CDFT Hamiltonian.
Definition: mixed_cdft_utils.F:1681
mixed_cdft_utils::mixed_cdft_read_block_diag
subroutine, public mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
Read input section related to block diagonalization of the mixed CDFT Hamiltonian matrix.
Definition: mixed_cdft_utils.F:1589
mixed_cdft_utils::mixed_cdft_redistribute_arrays
subroutine, public mixed_cdft_redistribute_arrays(force_env)
Redistribute arrays needed for an ET coupling calculation from individual CDFT states to the mixed CD...
Definition: mixed_cdft_utils.F:1041
mixed_cdft_utils::hfun_zero
subroutine, public hfun_zero(fun, th, just_zero, bounds, work)
Determine confinement bounds along confinement dir (hardcoded to be z) and determine the number of no...
Definition: mixed_cdft_utils.F:1519
mixed_cdft_utils::mixed_cdft_release_work
subroutine, public mixed_cdft_release_work(force_env)
Release storage reserved for mixed CDFT matrices.
Definition: mixed_cdft_utils.F:1463
mixed_cdft_utils::mixed_cdft_parse_settings
subroutine, public mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, settings, natom)
Parse settings for mixed cdft calculation and check their consistency.
Definition: mixed_cdft_utils.F:141
mixed_environment_types
Definition: mixed_environment_types.F:11
mixed_environment_types::get_mixed_env
subroutine, public get_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, mixed_energy, para_env, sub_para_env, subsys, input, results, cdft_control)
Get the MIXED environment.
Definition: mixed_environment_types.F:130
pw_env_methods
methods of pw_env that have dependence on qs_env
Definition: pw_env_methods.F:17
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_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_grid_types
Definition: pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition: pw_grid_types.F:28
pw_grids
This module defines the grid data type and some basic operations on it.
Definition: pw_grids.F:36
pw_grids::do_pw_grid_blocked_false
integer, parameter, public do_pw_grid_blocked_false
Definition: pw_grids.F:78
pw_grids::pw_grid_release
subroutine, public pw_grid_release(pw_grid)
releases the given pw grid
Definition: pw_grids.F:2133
pw_grids::pw_grid_setup
subroutine, public pw_grid_setup(cell_hmat, pw_grid, grid_span, cutoff, bounds, bounds_local, npts, spherical, odd, fft_usage, ncommensurate, icommensurate, blocked, ref_grid, rs_dims, iounit)
sets up a pw_grid
Definition: pw_grids.F:286
pw_grids::pw_grid_create
subroutine, public pw_grid_create(pw_grid, pe_group, local)
Initialize a PW grid with all defaults.
Definition: pw_grids.F:93
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_pool_types::pw_pool_create
subroutine, public pw_pool_create(pool, pw_grid, max_cache)
creates a pool for pw
Definition: pw_pool_types.F:173
qs_cdft_types
Defines CDFT control structures.
Definition: qs_cdft_types.F:14
qs_cdft_types::cdft_control_create
subroutine, public cdft_control_create(cdft_control)
create the cdft_control_type
Definition: qs_cdft_types.F:353
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_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_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::create_qs_kind_set
subroutine, public create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, force_env_section)
Read an atomic kind set data set from the input file.
Definition: qs_kind_types.F:2600
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::rs_grid_print
subroutine, public rs_grid_print(rs, iounit)
Print information on grids to output.
Definition: realspace_grid_types.F:583
realspace_grid_types::rs_grid_create
subroutine, public rs_grid_create(rs, desc)
...
Definition: realspace_grid_types.F:482
realspace_grid_types::rs_grid_create_descriptor
subroutine, public rs_grid_create_descriptor(desc, pw_grid, input_settings, border_points)
Determine the setup of real space grids - this is divided up into the creation of a descriptor and th...
Definition: realspace_grid_types.F:203