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