Line data Source code
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 : ! **************************************************************************************************
14 : MODULE mixed_cdft_utils
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE cell_types, ONLY: cell_type
17 : USE cp_array_utils, ONLY: cp_1d_i_p_type,&
18 : cp_1d_r_p_type,&
19 : cp_2d_r_p_type
20 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
21 : cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_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 :
127 : CONTAINS
128 :
129 : ! **************************************************************************************************
130 : !> \brief Parse settings for mixed cdft calculation and check their consistency
131 : !> \param force_env the force_env that holds the CDFT mixed_env
132 : !> \param mixed_env the mixed_env that holds the CDFT states
133 : !> \param mixed_cdft control section for mixed CDFT
134 : !> \param settings container for settings related to the mixed CDFT calculation
135 : !> \param natom the total number of atoms
136 : !> \par History
137 : !> 01.2017 created [Nico Holmberg]
138 : ! **************************************************************************************************
139 78 : 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 78 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: constraint_type
150 78 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: array_sizes
151 : LOGICAL :: is_match
152 78 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
153 : TYPE(cdft_control_type), POINTER :: cdft_control
154 78 : TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:, :) :: atoms
155 78 : 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 78 : NULLIFY (dft_control, qs_env, pw_env, auxbas_pw_pool, force_env_qs, &
163 78 : cdft_control)
164 : ! Allocate storage for temporaries used for checking settings consistency
165 78 : settings%max_nkinds = 30
166 78 : nforce_eval = SIZE(force_env%sub_force_env)
167 234 : ALLOCATE (settings%grid_span(nforce_eval))
168 234 : ALLOCATE (settings%npts(3, nforce_eval))
169 234 : ALLOCATE (settings%cutoff(nforce_eval))
170 156 : ALLOCATE (settings%rel_cutoff(nforce_eval))
171 156 : ALLOCATE (settings%spherical(nforce_eval))
172 234 : ALLOCATE (settings%rs_dims(2, nforce_eval))
173 156 : ALLOCATE (settings%odd(nforce_eval))
174 312 : ALLOCATE (settings%atoms(natom, nforce_eval))
175 78 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
176 96 : ALLOCATE (settings%coeffs(natom, nforce_eval))
177 168 : 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 234 : ALLOCATE (settings%si(6, nforce_eval))
183 234 : ALLOCATE (settings%sb(8, nforce_eval))
184 234 : ALLOCATE (settings%sr(5, nforce_eval))
185 234 : ALLOCATE (settings%cutoffs(settings%max_nkinds, nforce_eval))
186 156 : ALLOCATE (settings%radii(settings%max_nkinds, nforce_eval))
187 258 : settings%grid_span = 0
188 798 : settings%npts = 0
189 258 : settings%cutoff = 0.0_dp
190 258 : settings%rel_cutoff = 0.0_dp
191 258 : settings%spherical = 0
192 78 : settings%is_spherical = .FALSE.
193 618 : settings%rs_dims = 0
194 258 : settings%odd = 0
195 78 : settings%is_odd = .FALSE.
196 656 : settings%atoms = 0
197 1338 : settings%si = 0
198 1158 : settings%sr = 0.0_dp
199 1698 : settings%sb = .FALSE.
200 5658 : settings%cutoffs = 0.0_dp
201 5658 : settings%radii = 0.0_dp
202 : ! Get information from the sub_force_envs
203 258 : DO iforce_eval = 1, nforce_eval
204 180 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
205 154 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
206 154 : IF (mixed_env%do_mixed_qmmm_cdft) THEN
207 12 : qs_env => force_env_qs%qmmm_env%qs_env
208 : ELSE
209 142 : CALL force_env_get(force_env_qs, qs_env=qs_env)
210 : END IF
211 154 : CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
212 154 : 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 0 : "but CDFT constraints were not active in the QS section of all force_evals!")
216 : END IF
217 154 : cdft_control => dft_control%qs_control%cdft_control
218 154 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
219 1540 : settings%bo = auxbas_pw_pool%pw_grid%bounds_local
220 : ! Only the rank 0 process collects info about pw_grid and CDFT
221 232 : IF (force_env_qs%para_env%is_source()) THEN
222 : ! Grid settings
223 90 : settings%grid_span(iforce_eval) = auxbas_pw_pool%pw_grid%grid_span
224 360 : settings%npts(:, iforce_eval) = auxbas_pw_pool%pw_grid%npts
225 90 : settings%cutoff(iforce_eval) = auxbas_pw_pool%pw_grid%cutoff
226 90 : settings%rel_cutoff(iforce_eval) = dft_control%qs_control%relative_cutoff
227 90 : IF (auxbas_pw_pool%pw_grid%spherical) settings%spherical(iforce_eval) = 1
228 270 : settings%rs_dims(:, iforce_eval) = auxbas_pw_pool%pw_grid%para%group%num_pe_cart
229 90 : IF (auxbas_pw_pool%pw_grid%grid_span == HALFSPACE) settings%odd(iforce_eval) = 1
230 : ! Becke constraint atoms/coeffs
231 90 : 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 0 : "Use default values for MIXED\MAPPING.")
235 : END IF
236 251 : settings%atoms(1:cdft_control%natoms, iforce_eval) = cdft_control%atoms
237 90 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
238 72 : settings%coeffs(1:cdft_control%natoms, iforce_eval) = cdft_control%group(1)%coeff
239 : END IF
240 : ! Integer type settings
241 90 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
242 82 : settings%si(1, iforce_eval) = cdft_control%becke_control%cutoff_type
243 82 : settings%si(2, iforce_eval) = cdft_control%becke_control%cavity_shape
244 : END IF
245 90 : settings%si(3, iforce_eval) = dft_control%multiplicity
246 90 : settings%si(4, iforce_eval) = SIZE(cdft_control%group)
247 90 : settings%si(5, iforce_eval) = cdft_control%type
248 90 : IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
249 8 : settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%shape_function
250 8 : settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%gaussian_shape
251 : END IF
252 : ! Logicals
253 90 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
254 82 : settings%sb(1, iforce_eval) = cdft_control%becke_control%cavity_confine
255 82 : settings%sb(2, iforce_eval) = cdft_control%becke_control%should_skip
256 82 : settings%sb(3, iforce_eval) = cdft_control%becke_control%print_cavity
257 82 : settings%sb(4, iforce_eval) = cdft_control%becke_control%in_memory
258 82 : settings%sb(5, iforce_eval) = cdft_control%becke_control%adjust
259 82 : settings%sb(8, iforce_eval) = cdft_control%becke_control%use_bohr
260 : END IF
261 90 : IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
262 8 : settings%sb(8, iforce_eval) = cdft_control%hirshfeld_control%use_bohr
263 : END IF
264 90 : settings%sb(6, iforce_eval) = cdft_control%atomic_charges
265 90 : settings%sb(7, iforce_eval) = qs_env%has_unit_metric
266 : ! Reals
267 90 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
268 82 : settings%sr(1, iforce_eval) = cdft_control%becke_control%rcavity
269 82 : settings%sr(2, iforce_eval) = cdft_control%becke_control%rglobal
270 82 : settings%sr(3, iforce_eval) = cdft_control%becke_control%eps_cavity
271 : END IF
272 90 : IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
273 8 : settings%sr(2, iforce_eval) = cdft_control%hirshfeld_control%radius
274 : END IF
275 90 : settings%sr(4, iforce_eval) = dft_control%qs_control%eps_rho_rspace
276 90 : settings%sr(5, iforce_eval) = pw_env%cube_info(pw_env%auxbas_grid)%max_rad_ga
277 90 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
278 82 : IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
279 50 : nkinds = SIZE(cdft_control%becke_control%cutoffs_tmp)
280 50 : 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 0 : " your input is correct? If yes, please increase max_nkinds and recompile.")
285 : END IF
286 150 : settings%cutoffs(1:nkinds, iforce_eval) = cdft_control%becke_control%cutoffs_tmp(:)
287 : END IF
288 82 : IF (cdft_control%becke_control%adjust) THEN
289 52 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
290 52 : 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 0 : "match number of atomic kinds in the input coordinate file.")
294 : END IF
295 52 : nkinds = SIZE(cdft_control%becke_control%radii_tmp)
296 52 : 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 0 : " your input is correct? If yes, please increase max_nkinds and recompile.")
301 : END IF
302 156 : settings%radii(1:nkinds, iforce_eval) = cdft_control%becke_control%radii_tmp(:)
303 : END IF
304 : END IF
305 116 : IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
306 8 : IF (ASSOCIATED(cdft_control%hirshfeld_control%radii)) THEN
307 0 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
308 0 : 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 0 : "match number of atomic kinds in the input coordinate file.")
312 : END IF
313 0 : nkinds = SIZE(cdft_control%hirshfeld_control%radii)
314 0 : 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 0 : " your input is correct? If yes, please increase max_nkinds and recompile.")
319 : END IF
320 0 : 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 438 : CALL force_env%para_env%sum(settings%grid_span)
327 1518 : CALL force_env%para_env%sum(settings%npts)
328 438 : CALL force_env%para_env%sum(settings%cutoff)
329 438 : CALL force_env%para_env%sum(settings%rel_cutoff)
330 438 : CALL force_env%para_env%sum(settings%spherical)
331 1158 : CALL force_env%para_env%sum(settings%rs_dims)
332 438 : CALL force_env%para_env%sum(settings%odd)
333 78 : is_match = .TRUE.
334 180 : DO iforce_eval = 2, nforce_eval
335 102 : is_match = is_match .AND. (settings%grid_span(1) == settings%grid_span(iforce_eval))
336 102 : is_match = is_match .AND. (settings%npts(1, 1) == settings%npts(1, iforce_eval))
337 102 : is_match = is_match .AND. (settings%cutoff(1) == settings%cutoff(iforce_eval))
338 102 : is_match = is_match .AND. (settings%rel_cutoff(1) == settings%rel_cutoff(iforce_eval))
339 102 : is_match = is_match .AND. (settings%spherical(1) == settings%spherical(iforce_eval))
340 102 : is_match = is_match .AND. (settings%rs_dims(1, 1) == settings%rs_dims(1, iforce_eval))
341 102 : is_match = is_match .AND. (settings%rs_dims(2, 1) == settings%rs_dims(2, iforce_eval))
342 180 : is_match = is_match .AND. (settings%odd(1) == settings%odd(iforce_eval))
343 : END DO
344 78 : IF (.NOT. is_match) THEN
345 : CALL cp_abort(__LOCATION__, &
346 0 : "Mismatch detected in the &MGRID settings of the CDFT force_evals.")
347 : END IF
348 78 : IF (settings%spherical(1) == 1) settings%is_spherical = .TRUE.
349 78 : IF (settings%odd(1) == 1) settings%is_odd = .TRUE.
350 : ! Make sure CDFT settings are consistent
351 1234 : CALL force_env%para_env%sum(settings%atoms)
352 78 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
353 312 : CALL force_env%para_env%sum(settings%coeffs)
354 : END IF
355 78 : settings%ncdft = 0
356 242 : DO i = 1, SIZE(settings%atoms, 1)
357 398 : DO iforce_eval = 2, nforce_eval
358 398 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
359 48 : IF (settings%atoms(i, 1) /= settings%atoms(i, iforce_eval)) is_match = .FALSE.
360 48 : IF (settings%coeffs(i, 1) /= settings%coeffs(i, iforce_eval)) is_match = .FALSE.
361 : END IF
362 : END DO
363 242 : IF (settings%atoms(i, 1) /= 0) settings%ncdft = settings%ncdft + 1
364 : END DO
365 78 : 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 0 : "want to use nonidentical constraint definitions.")
371 : END IF
372 2598 : CALL force_env%para_env%sum(settings%si)
373 2238 : CALL force_env%para_env%sum(settings%sr)
374 702 : DO i = 1, SIZE(settings%sb, 1)
375 624 : CALL force_env%para_env%sum(settings%sb(i, 1))
376 1518 : DO iforce_eval = 2, nforce_eval
377 816 : CALL force_env%para_env%sum(settings%sb(i, iforce_eval))
378 1440 : IF (settings%sb(i, 1) .NEQV. settings%sb(i, iforce_eval)) is_match = .FALSE.
379 : END DO
380 : END DO
381 546 : DO i = 1, SIZE(settings%si, 1)
382 1158 : DO iforce_eval = 2, nforce_eval
383 1080 : IF (settings%si(i, 1) /= settings%si(i, iforce_eval)) is_match = .FALSE.
384 : END DO
385 : END DO
386 468 : DO i = 1, SIZE(settings%sr, 1)
387 978 : DO iforce_eval = 2, nforce_eval
388 900 : IF (settings%sr(i, 1) /= settings%sr(i, iforce_eval)) is_match = .FALSE.
389 : END DO
390 : END DO
391 78 : IF (.NOT. is_match) THEN
392 : CALL cp_abort(__LOCATION__, &
393 0 : "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 78 : IF (mixed_cdft%dlb .AND. .NOT. settings%sb(1, 1)) THEN
397 : CALL cp_abort(__LOCATION__, &
398 0 : "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 78 : IF (mixed_cdft%run_type /= mixed_cdft_parallel) THEN
402 : ! Get array sizes
403 270 : ALLOCATE (array_sizes(nforce_eval, settings%si(4, 1), 2))
404 54 : array_sizes(:, :, :) = 0
405 186 : DO iforce_eval = 1, nforce_eval
406 132 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
407 130 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
408 130 : IF (mixed_env%do_mixed_qmmm_cdft) THEN
409 8 : qs_env => force_env_qs%qmmm_env%qs_env
410 : ELSE
411 122 : CALL force_env_get(force_env_qs, qs_env=qs_env)
412 : END IF
413 130 : CALL get_qs_env(qs_env, dft_control=dft_control)
414 130 : cdft_control => dft_control%qs_control%cdft_control
415 184 : IF (force_env_qs%para_env%is_source()) THEN
416 136 : DO igroup = 1, SIZE(cdft_control%group)
417 68 : array_sizes(iforce_eval, igroup, 1) = SIZE(cdft_control%group(igroup)%atoms)
418 134 : 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 54 : CALL force_env%para_env%sum(array_sizes)
424 438 : IF (ANY(array_sizes(:, :, 1) /= array_sizes(1, 1, 1)) .OR. &
425 : ANY(array_sizes(:, :, 2) /= array_sizes(1, 1, 2))) THEN
426 0 : mixed_cdft%identical_constraints = .FALSE.
427 : END IF
428 : ! Check constraint definitions
429 54 : IF (mixed_cdft%identical_constraints) THEN
430 : ! Prepare temporary storage
431 408 : ALLOCATE (atoms(nforce_eval, settings%si(4, 1)))
432 354 : ALLOCATE (coeff(nforce_eval, settings%si(4, 1)))
433 216 : ALLOCATE (constraint_type(nforce_eval, settings%si(4, 1)))
434 54 : constraint_type(:, :) = 0
435 186 : DO iforce_eval = 1, nforce_eval
436 268 : DO i = 1, settings%si(4, 1)
437 136 : NULLIFY (atoms(iforce_eval, i)%array)
438 : NULLIFY (coeff(iforce_eval, i)%array)
439 408 : ALLOCATE (atoms(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
440 408 : ALLOCATE (coeff(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
441 370 : atoms(iforce_eval, i)%array(:) = 0
442 502 : coeff(iforce_eval, i)%array(:) = 0
443 : END DO
444 : ! Get constraint definitions
445 132 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
446 130 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
447 130 : IF (mixed_env%do_mixed_qmmm_cdft) THEN
448 8 : qs_env => force_env_qs%qmmm_env%qs_env
449 : ELSE
450 122 : CALL force_env_get(force_env_qs, qs_env=qs_env)
451 : END IF
452 130 : CALL get_qs_env(qs_env, dft_control=dft_control)
453 130 : cdft_control => dft_control%qs_control%cdft_control
454 184 : IF (force_env_qs%para_env%is_source()) THEN
455 136 : DO i = 1, settings%si(4, 1)
456 185 : atoms(iforce_eval, i)%array(:) = cdft_control%group(i)%atoms
457 185 : coeff(iforce_eval, i)%array(:) = cdft_control%group(i)%coeff
458 134 : 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 108 : DO i = 1, settings%si(4, 1)
464 192 : DO iforce_eval = 1, nforce_eval
465 604 : CALL force_env%para_env%sum(atoms(iforce_eval, i)%array)
466 604 : CALL force_env%para_env%sum(coeff(iforce_eval, i)%array)
467 192 : CALL force_env%para_env%sum(constraint_type(iforce_eval, i))
468 : END DO
469 134 : DO iforce_eval = 2, nforce_eval
470 206 : DO iatom = 1, SIZE(atoms(1, i)%array)
471 128 : IF (atoms(1, i)%array(iatom) /= atoms(iforce_eval, i)%array(iatom)) THEN
472 0 : mixed_cdft%identical_constraints = .FALSE.
473 : END IF
474 128 : IF (coeff(1, i)%array(iatom) /= coeff(iforce_eval, i)%array(iatom)) THEN
475 2 : mixed_cdft%identical_constraints = .FALSE.
476 : END IF
477 206 : IF (.NOT. mixed_cdft%identical_constraints) EXIT
478 : END DO
479 80 : IF (constraint_type(1, i) /= constraint_type(iforce_eval, i)) THEN
480 0 : mixed_cdft%identical_constraints = .FALSE.
481 : END IF
482 134 : IF (.NOT. mixed_cdft%identical_constraints) EXIT
483 : END DO
484 108 : IF (.NOT. mixed_cdft%identical_constraints) EXIT
485 : END DO
486 : ! Deallocate temporary storage
487 186 : DO iforce_eval = 1, nforce_eval
488 322 : DO i = 1, settings%si(4, 1)
489 136 : DEALLOCATE (atoms(iforce_eval, i)%array)
490 268 : DEALLOCATE (coeff(iforce_eval, i)%array)
491 : END DO
492 : END DO
493 54 : DEALLOCATE (atoms)
494 54 : DEALLOCATE (coeff)
495 54 : DEALLOCATE (constraint_type)
496 : END IF
497 54 : DEALLOCATE (array_sizes)
498 : END IF
499 : ! Deallocate some arrays that are no longer needed
500 78 : IF (mixed_cdft%identical_constraints .AND. mixed_cdft%run_type /= mixed_cdft_parallel_nobuild) THEN
501 246 : DO iforce_eval = 1, nforce_eval
502 172 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
503 148 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
504 148 : IF (mixed_env%do_mixed_qmmm_cdft) THEN
505 12 : qs_env => force_env_qs%qmmm_env%qs_env
506 : ELSE
507 136 : CALL force_env_get(force_env_qs, qs_env=qs_env)
508 : END IF
509 148 : CALL get_qs_env(qs_env, dft_control=dft_control)
510 148 : cdft_control => dft_control%qs_control%cdft_control
511 222 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
512 24 : IF (.NOT. dft_control%qs_control%gapw) THEN
513 48 : DO i = 1, SIZE(cdft_control%group)
514 24 : DEALLOCATE (cdft_control%group(i)%coeff)
515 48 : DEALLOCATE (cdft_control%group(i)%atoms)
516 : END DO
517 24 : IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
518 : END IF
519 124 : ELSE IF (mixed_cdft%run_type == mixed_cdft_serial) THEN
520 124 : IF (iforce_eval == 1) CYCLE
521 150 : DO igroup = 1, SIZE(cdft_control%group)
522 150 : IF (.NOT. dft_control%qs_control%gapw) THEN
523 76 : DEALLOCATE (cdft_control%group(igroup)%coeff)
524 76 : DEALLOCATE (cdft_control%group(igroup)%atoms)
525 : END IF
526 : END DO
527 74 : IF (cdft_control%type == outer_scf_becke_constraint) THEN
528 68 : IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
529 68 : IF (cdft_control%becke_control%cavity_confine) THEN
530 62 : CALL release_hirshfeld_type(cdft_control%becke_control%cavity_env)
531 : END IF
532 68 : IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
533 42 : DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
534 : END IF
535 68 : IF (cdft_control%becke_control%adjust) THEN
536 44 : DEALLOCATE (cdft_control%becke_control%radii_tmp)
537 : END IF
538 : END IF
539 : END IF
540 : END DO
541 : END IF
542 :
543 156 : 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 78 : 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 78 : NULLIFY (cdft_control)
563 78 : is_match = .TRUE.
564 : ! Transfer global settings
565 78 : mixed_cdft%multiplicity = settings%si(3, 1)
566 78 : mixed_cdft%has_unit_metric = settings%sb(7, 1)
567 78 : mixed_cdft%eps_rho_rspace = settings%sr(4, 1)
568 78 : mixed_cdft%nconstraint = settings%si(4, 1)
569 78 : settings%radius = settings%sr(5, 1)
570 : ! Transfer settings only needed if the constraint should be built in parallel
571 78 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
572 24 : IF (settings%sb(6, 1)) THEN
573 : CALL cp_abort(__LOCATION__, &
574 0 : "Calculation of atomic Becke charges not supported with parallel mode mixed CDFT")
575 : END IF
576 24 : IF (mixed_cdft%nconstraint /= 1) THEN
577 : CALL cp_abort(__LOCATION__, &
578 0 : "Parallel mode mixed CDFT does not yet support multiple constraints.")
579 : END IF
580 :
581 24 : IF (settings%si(5, 1) /= outer_scf_becke_constraint) THEN
582 : CALL cp_abort(__LOCATION__, &
583 0 : "Parallel mode mixed CDFT does not support Hirshfeld constraints.")
584 : END IF
585 :
586 72 : ALLOCATE (mixed_cdft%cdft_control)
587 24 : CALL cdft_control_create(mixed_cdft%cdft_control)
588 24 : cdft_control => mixed_cdft%cdft_control
589 72 : ALLOCATE (cdft_control%atoms(settings%ncdft))
590 72 : cdft_control%atoms = settings%atoms(1:settings%ncdft, 1)
591 48 : ALLOCATE (cdft_control%group(1))
592 72 : ALLOCATE (cdft_control%group(1)%atoms(settings%ncdft))
593 72 : ALLOCATE (cdft_control%group(1)%coeff(settings%ncdft))
594 24 : NULLIFY (cdft_control%group(1)%weight)
595 24 : NULLIFY (cdft_control%group(1)%gradients)
596 24 : NULLIFY (cdft_control%group(1)%integrated)
597 72 : cdft_control%group(1)%atoms = cdft_control%atoms
598 72 : cdft_control%group(1)%coeff = settings%coeffs(1:settings%ncdft, 1)
599 24 : cdft_control%natoms = settings%ncdft
600 24 : cdft_control%atomic_charges = settings%sb(6, 1)
601 24 : cdft_control%becke_control%cutoff_type = settings%si(1, 1)
602 24 : cdft_control%becke_control%cavity_confine = settings%sb(1, 1)
603 24 : cdft_control%becke_control%should_skip = settings%sb(2, 1)
604 24 : cdft_control%becke_control%print_cavity = settings%sb(3, 1)
605 24 : cdft_control%becke_control%in_memory = settings%sb(4, 1)
606 24 : cdft_control%becke_control%adjust = settings%sb(5, 1)
607 24 : cdft_control%becke_control%cavity_shape = settings%si(2, 1)
608 24 : cdft_control%becke_control%use_bohr = settings%sb(8, 1)
609 24 : cdft_control%becke_control%rcavity = settings%sr(1, 1)
610 24 : cdft_control%becke_control%rglobal = settings%sr(2, 1)
611 24 : cdft_control%becke_control%eps_cavity = settings%sr(3, 1)
612 24 : nkinds = 0
613 24 : IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
614 2250 : CALL force_env%para_env%sum(settings%cutoffs)
615 558 : DO i = 1, SIZE(settings%cutoffs, 1)
616 540 : IF (settings%cutoffs(i, 1) /= settings%cutoffs(i, 2)) is_match = .FALSE.
617 558 : IF (settings%cutoffs(i, 1) /= 0.0_dp) nkinds = nkinds + 1
618 : END DO
619 18 : IF (.NOT. is_match) THEN
620 : CALL cp_abort(__LOCATION__, &
621 : "Mismatch detected in the &BECKE_CONSTRAINT "// &
622 0 : "&ELEMENT_CUTOFF settings of the two force_evals.")
623 : END IF
624 54 : ALLOCATE (cdft_control%becke_control%cutoffs_tmp(nkinds))
625 54 : cdft_control%becke_control%cutoffs_tmp = settings%cutoffs(1:nkinds, 1)
626 : END IF
627 24 : nkinds = 0
628 24 : IF (cdft_control%becke_control%adjust) THEN
629 2250 : CALL force_env%para_env%sum(settings%radii)
630 558 : DO i = 1, SIZE(settings%radii, 1)
631 540 : IF (settings%radii(i, 1) /= settings%radii(i, 2)) is_match = .FALSE.
632 558 : IF (settings%radii(i, 1) /= 0.0_dp) nkinds = nkinds + 1
633 : END DO
634 18 : IF (.NOT. is_match) THEN
635 : CALL cp_abort(__LOCATION__, &
636 : "Mismatch detected in the &BECKE_CONSTRAINT "// &
637 0 : "&ATOMIC_RADII settings of the two force_evals.")
638 : END IF
639 54 : ALLOCATE (cdft_control%becke_control%radii(nkinds))
640 54 : cdft_control%becke_control%radii = settings%radii(1:nkinds, 1)
641 : END IF
642 : END IF
643 :
644 78 : 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 78 : 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 78 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds
667 : INTEGER, DIMENSION(2, 3) :: bo, bo_mixed
668 : INTEGER, DIMENSION(3) :: higher_grid_layout
669 156 : INTEGER, DIMENSION(:), POINTER :: i_force_eval, mixed_rs_dims, recvbuffer, &
670 156 : recvbuffer2, sendbuffer
671 : LOGICAL :: is_match
672 78 : 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 312 : TYPE(mp_request_type), DIMENSION(3) :: req
678 : TYPE(pw_env_type), POINTER :: pw_env
679 : TYPE(pw_grid_type), POINTER :: pw_grid
680 78 : 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 78 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
684 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
685 78 : POINTER :: rs_descs
686 : TYPE(realspace_grid_input_type) :: input_settings
687 78 : 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 78 : NULLIFY (cell_mix, subsys_mix, force_env_section, subsys_section, &
692 78 : print_section, root_section, kind_section, force_env_sections, &
693 78 : rs_grid_section, auxbas_pw_pool, pw_env, pw_pools, pw_grid, &
694 78 : sendbuffer, qs_env, mixed_rs_dims, i_force_eval, recvbuffer, &
695 78 : recvbuffer2, globenv, atomic_kind_set, qs_kind_set, rs_descs, &
696 78 : rs_grids)
697 :
698 156 : logger => cp_get_default_logger()
699 78 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
700 78 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
701 78 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
702 78 : is_match = .TRUE.
703 78 : nforce_eval = SIZE(force_env%sub_force_env)
704 78 : ncpu = force_env%para_env%num_pe
705 : ! Get infos about the mixed subsys
706 78 : IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
707 : CALL force_env_get(force_env=force_env, &
708 70 : subsys=subsys_mix)
709 : ELSE
710 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
711 8 : cp_subsys=subsys_mix)
712 : END IF
713 : ! Init structures only needed when the CDFT states are treated in parallel
714 78 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
715 : ! Start building the mixed auxbas_pw_pool
716 24 : 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 24 : mixed_cdft%is_pencil = .FALSE. ! Flag to control the first two mappings
731 24 : 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 24 : IF (ncpu/2 > settings%npts(1, 1)) THEN
735 0 : CPABORT("ncpu/2 => nx: decrease ncpu or disable xc_smoothing")
736 : END IF
737 : !
738 24 : ALLOCATE (mixed_rs_dims(2))
739 24 : IF (settings%rs_dims(2, 1) /= 1) mixed_cdft%is_pencil = .TRUE.
740 24 : IF (.NOT. mixed_cdft%is_pencil .AND. ncpu > settings%npts(1, 1)) mixed_cdft%is_special = .TRUE.
741 24 : IF (mixed_cdft%is_special) THEN
742 0 : mixed_rs_dims = [-1, -1]
743 24 : ELSE IF (mixed_cdft%is_pencil) THEN
744 0 : mixed_rs_dims = [settings%rs_dims(1, 1), 2*settings%rs_dims(2, 1)]
745 : ELSE
746 72 : mixed_rs_dims = [2*settings%rs_dims(1, 1), 1]
747 : END IF
748 24 : IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
749 : CALL force_env_get(force_env=force_env, &
750 20 : cell=cell_mix)
751 : ELSE
752 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
753 4 : 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 24 : iounit=iounit)
761 : ! Check if the layout was successfully created
762 24 : IF (mixed_cdft%is_special) THEN
763 0 : IF (.NOT. pw_grid%para%group%num_pe_cart(2) /= 1) is_match = .FALSE.
764 24 : ELSE IF (mixed_cdft%is_pencil) THEN
765 0 : IF (.NOT. pw_grid%para%group%num_pe_cart(1) == mixed_rs_dims(1)) is_match = .FALSE.
766 : ELSE
767 24 : 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 0 : "of processors or disabling xc_smoothing.")
774 : END IF
775 24 : DEALLOCATE (mixed_rs_dims)
776 : ! Create the pool
777 240 : bo_mixed = pw_grid%bounds_local
778 48 : ALLOCATE (pw_pools(1))
779 24 : NULLIFY (pw_pools(1)%pool)
780 24 : CALL pw_pool_create(pw_pools(1)%pool, pw_grid=pw_grid)
781 : ! Initialize Gaussian cavity confinement
782 24 : IF (mixed_cdft%cdft_control%becke_control%cavity_confine) THEN
783 22 : 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 22 : 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 24 : 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 22 : "PRINT%GRID_INFORMATION")
795 22 : 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 22 : print_section=print_section)
800 44 : ALLOCATE (rs_descs(1))
801 374 : ALLOCATE (rs_grids(1))
802 638 : ALLOCATE (mixed_cdft%pw_env%cube_info(1))
803 22 : higher_grid_layout = [-1, -1, -1]
804 22 : 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 22 : pw_grid%orthorhombic, settings%radius)
809 22 : NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
810 22 : CALL force_env_get(force_env, root_section=root_section)
811 22 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
812 22 : 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 22 : i_force_eval(2), i_force_eval(2))
815 22 : 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 22 : higher_grid_layout=higher_grid_layout)
820 22 : NULLIFY (rs_descs(1)%rs_desc)
821 22 : CALL rs_grid_create_descriptor(rs_descs(1)%rs_desc, pw_grid, input_settings)
822 22 : IF (rs_descs(1)%rs_desc%distributed) higher_grid_layout = rs_descs(1)%rs_desc%group_dim
823 22 : CALL rs_grid_create(rs_grids(1), rs_descs(1)%rs_desc)
824 22 : CALL rs_grid_print(rs_grids(1), iounit)
825 22 : mixed_cdft%pw_env%rs_descs => rs_descs
826 22 : 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 22 : i_rep_section=i_force_eval(1))
830 22 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
831 22 : NULLIFY (qs_kind_set)
832 22 : 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 22 : force_env%para_env, force_env_section, silent=.FALSE.)
835 22 : mixed_cdft%qs_kind_set => qs_kind_set
836 22 : DEALLOCATE (i_force_eval)
837 22 : CALL section_vals_release(force_env_section)
838 : END IF
839 : CALL force_env_get(force_env=force_env, &
840 24 : force_env_section=force_env_section)
841 24 : CALL pw_grid_release(pw_grid)
842 24 : mixed_cdft%pw_env%auxbas_grid = 1
843 24 : NULLIFY (mixed_cdft%pw_env%pw_pools)
844 24 : mixed_cdft%pw_env%pw_pools => pw_pools
845 240 : bo = settings%bo
846 : ! Determine which processors need to exchange data when redistributing the weight/gradient
847 24 : IF (.NOT. mixed_cdft%is_special) THEN
848 24 : ALLOCATE (mixed_cdft%dest_list(2))
849 24 : ALLOCATE (mixed_cdft%source_list(2))
850 24 : imap = force_env%para_env%mepos/2
851 72 : 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 24 : MODULO(force_env%para_env%mepos, force_env%para_env%num_pe/2)
854 72 : mixed_cdft%source_list = [imap, imap + 1]
855 : ! Determine bounds of the data that is replicated
856 24 : ALLOCATE (mixed_cdft%recv_bo(4))
857 24 : ALLOCATE (sendbuffer(2), recvbuffer(2), recvbuffer2(2))
858 24 : IF (mixed_cdft%is_pencil) THEN
859 0 : sendbuffer = [bo_mixed(1, 2), bo_mixed(2, 2)]
860 : ELSE
861 72 : 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 24 : request=req(1))
866 : CALL force_env%para_env%irecv(msgout=recvbuffer, source=mixed_cdft%source_list(1), &
867 24 : request=req(2))
868 : CALL force_env%para_env%irecv(msgout=recvbuffer2, source=mixed_cdft%source_list(2), &
869 24 : request=req(3))
870 24 : CALL req(1)%wait()
871 : CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dest_list(2), &
872 24 : request=req(1))
873 24 : CALL mp_waitall(req)
874 120 : mixed_cdft%recv_bo(1:2) = recvbuffer
875 120 : mixed_cdft%recv_bo(3:4) = recvbuffer2
876 24 : DEALLOCATE (sendbuffer, recvbuffer, recvbuffer2)
877 : ELSE
878 0 : IF (mixed_env%do_mixed_qmmm_cdft) THEN
879 0 : qs_env => force_env_qs%qmmm_env%qs_env
880 : ELSE
881 0 : CALL force_env_get(force_env_qs, qs_env=qs_env)
882 : END IF
883 0 : CALL get_qs_env(qs_env, pw_env=pw_env)
884 0 : 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 0 : ALLOCATE (bounds(0:auxbas_pw_pool%pw_grid%para%group%num_pe - 1, 1:2))
888 0 : DO i = 0, auxbas_pw_pool%pw_grid%para%group%num_pe - 1
889 0 : bounds(i, 1:2) = auxbas_pw_pool%pw_grid%para%bo(1:2, 1, i, 1)
890 0 : 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 0 : ntargets = 0
895 0 : offset = -1
896 0 : DO i = 0, auxbas_pw_pool%pw_grid%para%group%num_pe - 1
897 0 : IF ((bounds(i, 1) >= bo_mixed(1, 1) .AND. bounds(i, 1) <= bo_mixed(2, 1)) .OR. &
898 0 : (bounds(i, 2) >= bo_mixed(1, 1) .AND. bounds(i, 2) <= bo_mixed(2, 1))) THEN
899 0 : ntargets = ntargets + 1
900 0 : IF (offset == -1) offset = i
901 0 : ELSE IF (bounds(i, 2) > bo_mixed(2, 1)) THEN
902 : EXIT
903 : ELSE
904 0 : CYCLE
905 : END IF
906 : END DO
907 0 : ALLOCATE (mixed_cdft%dest_list(ntargets))
908 0 : ALLOCATE (mixed_cdft%dest_list_bo(2, ntargets))
909 : ! now determine the actual grid points to send
910 0 : j = 1
911 0 : DO i = offset, offset + ntargets - 1
912 0 : 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 0 : bo_mixed(2, 1) + (bounds(i, 2) - bo_mixed(2, 1))]
915 0 : j = j + 1
916 : END DO
917 0 : 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 0 : mixed_cdft%dest_list_save = mixed_cdft%dest_list
920 0 : 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 0 : DEALLOCATE (bounds)
924 0 : ALLOCATE (bounds(0:pw_pools(1)%pool%pw_grid%para%group%num_pe - 1, 1:4))
925 0 : DO i = 0, pw_pools(1)%pool%pw_grid%para%group%num_pe - 1
926 0 : bounds(i, 1:2) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 1, i, 1)
927 0 : bounds(i, 3:4) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 2, i, 1)
928 0 : bounds(i, 1:2) = bounds(i, 1:2) - pw_pools(1)%pool%pw_grid%npts(1)/2 - 1
929 0 : bounds(i, 3:4) = bounds(i, 3:4) - pw_pools(1)%pool%pw_grid%npts(2)/2 - 1
930 : END DO
931 0 : ntargets = 0
932 0 : offset = -1
933 0 : DO i = 0, pw_pools(1)%pool%pw_grid%para%group%num_pe - 1
934 0 : IF ((bo(1, 1) >= bounds(i, 1) .AND. bo(1, 1) <= bounds(i, 2)) .OR. &
935 0 : (bo(2, 1) >= bounds(i, 1) .AND. bo(2, 1) <= bounds(i, 2))) THEN
936 0 : ntargets = ntargets + 1
937 0 : IF (offset == -1) offset = i
938 0 : ELSE IF (bo(2, 1) < bounds(i, 1)) THEN
939 : EXIT
940 : ELSE
941 0 : CYCLE
942 : END IF
943 : END DO
944 0 : ALLOCATE (mixed_cdft%source_list(ntargets))
945 0 : ALLOCATE (mixed_cdft%source_list_bo(4, ntargets))
946 0 : j = 1
947 0 : DO i = offset, offset + ntargets - 1
948 0 : mixed_cdft%source_list(j) = i
949 0 : 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 0 : bounds(i, 3), bounds(i, 4)]
952 0 : 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 0 : bounds(i, 3), bounds(i, 4)]
955 : ELSE
956 : mixed_cdft%source_list_bo(:, j) = [bounds(i, 1), bo(2, 1), &
957 0 : bounds(i, 3), bounds(i, 4)]
958 : END IF
959 0 : j = j + 1
960 : END DO
961 0 : 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 0 : mixed_cdft%source_list_save = mixed_cdft%source_list
964 0 : mixed_cdft%source_bo_save = mixed_cdft%source_list_bo
965 0 : 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 54 : CALL force_env_get(force_env, root_section=root_section)
972 240 : ALLOCATE (mixed_cdft%sub_logger(nforce_eval - 1))
973 132 : DO i = 1, nforce_eval - 1
974 78 : IF (force_env%para_env%is_source()) THEN
975 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", &
976 39 : c_val=input_file_path)
977 39 : lp = LEN_TRIM(input_file_path)
978 39 : input_file_path(lp + 1:LEN(input_file_path)) = "-r-"//ADJUSTL(cp_to_string(i + 1))
979 39 : lp = LEN_TRIM(input_file_path)
980 39 : 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 39 : unit_number=unit_nr)
984 : ELSE
985 39 : 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 78 : 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 78 : c_val=c_val)
994 78 : IF (c_val /= "") THEN
995 : CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
996 0 : local_filename=TRIM(c_val)//"_localLog")
997 : END IF
998 78 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
999 78 : IF (c_val /= "") THEN
1000 : CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
1001 78 : local_filename=TRIM(c_val)//"_localLog")
1002 : END IF
1003 78 : IF (LEN_TRIM(c_val) > default_string_length) THEN
1004 0 : CPWARN("The mixed CDFT project name will be truncated.")
1005 : END IF
1006 78 : 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 132 : i_val=mixed_cdft%sub_logger(i)%p%iter_info%print_level)
1009 : END DO
1010 54 : IF (mixed_cdft%wfn_overlap_method) THEN
1011 : ! qs_kind_set
1012 6 : NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
1013 6 : CALL force_env_get(force_env, root_section=root_section)
1014 6 : force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
1015 6 : 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 6 : i_force_eval(2), i_force_eval(2))
1018 : subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
1019 6 : i_rep_section=i_force_eval(1))
1020 6 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
1021 6 : NULLIFY (qs_kind_set)
1022 6 : 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 6 : force_env%para_env, force_env_section, silent=.FALSE.)
1025 6 : mixed_cdft%qs_kind_set => qs_kind_set
1026 6 : DEALLOCATE (i_force_eval)
1027 6 : CALL section_vals_release(force_env_section)
1028 6 : mixed_cdft%qs_kind_set => qs_kind_set
1029 : END IF
1030 : CALL force_env_get(force_env=force_env, &
1031 54 : force_env_section=force_env_section)
1032 : END IF
1033 : ! Deallocate settings temporaries
1034 78 : DEALLOCATE (settings%grid_span)
1035 78 : DEALLOCATE (settings%npts)
1036 78 : DEALLOCATE (settings%spherical)
1037 78 : DEALLOCATE (settings%rs_dims)
1038 78 : DEALLOCATE (settings%odd)
1039 78 : DEALLOCATE (settings%atoms)
1040 78 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
1041 24 : DEALLOCATE (settings%coeffs)
1042 : END IF
1043 78 : DEALLOCATE (settings%cutoffs)
1044 78 : DEALLOCATE (settings%radii)
1045 78 : DEALLOCATE (settings%si)
1046 78 : DEALLOCATE (settings%sr)
1047 78 : DEALLOCATE (settings%sb)
1048 78 : DEALLOCATE (settings%cutoff)
1049 78 : DEALLOCATE (settings%rel_cutoff)
1050 : ! Setup mixed blacs_env for redistributing arrays during ET coupling calculation
1051 78 : IF (mixed_env%do_mixed_et) THEN
1052 78 : NULLIFY (root_section)
1053 78 : 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 78 : globenv%blacs_repeatable)
1056 : END IF
1057 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1058 78 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1059 :
1060 390 : 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 100 : 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 100 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ncol_mo, nrow_mo
1077 : LOGICAL :: uniform_occupation
1078 100 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_occupation_numbers
1079 100 : 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 100 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: matrix_p_tmp, mixed_matrix_p_tmp, &
1085 100 : mixed_wmat_tmp, mo_coeff_tmp, wmat_tmp
1086 100 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
1087 100 : 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 100 : NULLIFY (mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
1097 100 : fm_struct_wmat, fm_struct_overlap, fm_struct_tmp, &
1098 100 : mixed_mo_coeff, mixed_matrix_s, density_matrix, blacs_env, w_matrix, force_env_qs)
1099 0 : CPASSERT(ASSOCIATED(force_env))
1100 100 : mixed_env => force_env%mixed_env
1101 100 : nforce_eval = SIZE(force_env%sub_force_env)
1102 100 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1103 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1104 100 : CALL mixed_cdft_work_type_init(mixed_cdft%matrix)
1105 : ! Get nspins and query for non-uniform occupation numbers
1106 300 : ALLOCATE (has_occupation_numbers(nforce_eval))
1107 100 : has_occupation_numbers = .FALSE.
1108 324 : DO iforce_eval = 1, nforce_eval
1109 224 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
1110 186 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
1111 186 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1112 24 : qs_env => force_env_qs%qmmm_env%qs_env
1113 : ELSE
1114 162 : CALL force_env_get(force_env_qs, qs_env=qs_env)
1115 : END IF
1116 186 : CALL get_qs_env(qs_env, dft_control=dft_control)
1117 186 : CPASSERT(ASSOCIATED(dft_control))
1118 186 : nspins = dft_control%nspins
1119 286 : IF (force_env_qs%para_env%is_source()) THEN
1120 150 : has_occupation_numbers(iforce_eval) = ALLOCATED(dft_control%qs_control%cdft_control%occupations)
1121 : END IF
1122 : END DO
1123 100 : CALL force_env%para_env%sum(has_occupation_numbers(1))
1124 224 : DO iforce_eval = 2, nforce_eval
1125 124 : CALL force_env%para_env%sum(has_occupation_numbers(iforce_eval))
1126 224 : IF (has_occupation_numbers(1) .NEQV. has_occupation_numbers(iforce_eval)) THEN
1127 : CALL cp_abort(__LOCATION__, &
1128 0 : "Mixing of uniform and non-uniform occupations is not allowed.")
1129 : END IF
1130 : END DO
1131 100 : uniform_occupation = .NOT. has_occupation_numbers(1)
1132 100 : DEALLOCATE (has_occupation_numbers)
1133 : ! Get number of weight functions per state as well as the type of each constraint
1134 100 : nvar = SIZE(dft_control%qs_control%cdft_control%target)
1135 100 : IF (.NOT. ALLOCATED(mixed_cdft%constraint_type)) THEN
1136 312 : ALLOCATE (mixed_cdft%constraint_type(nvar, nforce_eval))
1137 442 : mixed_cdft%constraint_type(:, :) = 0
1138 78 : IF (mixed_cdft%identical_constraints) THEN
1139 154 : DO ivar = 1, nvar
1140 : mixed_cdft%constraint_type(ivar, :) = &
1141 334 : dft_control%qs_control%cdft_control%group(ivar)%constraint_type
1142 : END DO
1143 : ELSE
1144 : ! Possibly couple spin and charge constraints
1145 6 : DO iforce_eval = 1, nforce_eval
1146 4 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
1147 4 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1148 0 : qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1149 : ELSE
1150 4 : CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1151 : END IF
1152 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
1153 6 : IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1154 4 : DO ivar = 1, nvar
1155 : mixed_cdft%constraint_type(ivar, iforce_eval) = &
1156 4 : dft_control%qs_control%cdft_control%group(ivar)%constraint_type
1157 : END DO
1158 : END IF
1159 : END DO
1160 2 : 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 1030 : ALLOCATE (mixed_cdft%matrix%mixed_mo_coeff(nforce_eval, nspins))
1165 100 : mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
1166 730 : ALLOCATE (mixed_cdft%matrix%w_matrix(nforce_eval, nvar))
1167 100 : w_matrix => mixed_cdft%matrix%w_matrix
1168 100 : CALL dbcsr_init_p(mixed_cdft%matrix%mixed_matrix_s)
1169 100 : mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
1170 100 : IF (mixed_cdft%calculate_metric) THEN
1171 144 : ALLOCATE (mixed_cdft%matrix%density_matrix(nforce_eval, nspins))
1172 14 : density_matrix => mixed_cdft%matrix%density_matrix
1173 : END IF
1174 1560 : ALLOCATE (mo_coeff_tmp(nforce_eval, nspins), wmat_tmp(nforce_eval, nvar))
1175 400 : ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
1176 216 : IF (mixed_cdft%calculate_metric) ALLOCATE (matrix_p_tmp(nforce_eval, nspins))
1177 100 : IF (.NOT. uniform_occupation) THEN
1178 140 : ALLOCATE (mixed_cdft%occupations(nforce_eval, nspins))
1179 126 : ALLOCATE (occno_tmp(nforce_eval, nspins))
1180 : END IF
1181 324 : DO iforce_eval = 1, nforce_eval
1182 : ! Temporary arrays need to be nulled on every process
1183 660 : 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 660 : IF (.NOT. uniform_occupation) THEN
1188 56 : NULLIFY (occno_tmp(iforce_eval, ispin)%array)
1189 : END IF
1190 : END DO
1191 224 : 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 186 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
1195 186 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1196 24 : qs_env => force_env_qs%qmmm_env%qs_env
1197 : ELSE
1198 162 : CALL force_env_get(force_env_qs, qs_env=qs_env)
1199 : END IF
1200 186 : 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 186 : nfullrows_total=nrow_overlap, nfullcols_total=ncol_overlap)
1204 : CALL dbcsr_get_info(dft_control%qs_control%cdft_control%wmat(1)%matrix, &
1205 186 : nfullrows_total=nrow_wmat, nfullcols_total=ncol_wmat)
1206 : ! MO Coefficients
1207 548 : DO ispin = 1, nspins
1208 : CALL cp_fm_get_info(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
1209 362 : 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 362 : //TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1214 : CALL cp_fm_to_fm(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
1215 548 : mo_coeff_tmp(iforce_eval, ispin))
1216 : END DO
1217 186 : 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 186 : square_blocks=.TRUE.)
1222 376 : DO ivar = 1, nvar
1223 190 : CALL cp_fm_create(wmat_tmp(iforce_eval, ivar), fm_struct_tmp, name="w_matrix")
1224 190 : CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%wmat(ivar)%matrix, desymm_tmp)
1225 190 : CALL copy_dbcsr_to_fm(desymm_tmp, wmat_tmp(iforce_eval, ivar))
1226 190 : CALL dbcsr_release(desymm_tmp)
1227 376 : CALL dbcsr_release_p(dft_control%qs_control%cdft_control%wmat(ivar)%matrix)
1228 : END DO
1229 186 : DEALLOCATE (dft_control%qs_control%cdft_control%wmat)
1230 186 : 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 186 : 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 81 : para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
1236 81 : CALL cp_fm_create(matrix_s_tmp, fm_struct_tmp, name="s_matrix")
1237 81 : CALL cp_fm_struct_release(fm_struct_tmp)
1238 81 : CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_s%matrix, desymm_tmp)
1239 81 : CALL copy_dbcsr_to_fm(desymm_tmp, matrix_s_tmp)
1240 81 : CALL dbcsr_release(desymm_tmp)
1241 : END IF
1242 186 : CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_s%matrix)
1243 : ! Density_matrix (dbcsr -> fm)
1244 186 : IF (mixed_cdft%calculate_metric) THEN
1245 72 : 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 48 : para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
1250 48 : CALL cp_fm_create(matrix_p_tmp(iforce_eval, ispin), fm_struct_tmp, name="dm_matrix")
1251 48 : CALL cp_fm_struct_release(fm_struct_tmp)
1252 48 : CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix, desymm_tmp)
1253 48 : CALL copy_dbcsr_to_fm(desymm_tmp, matrix_p_tmp(iforce_eval, ispin))
1254 48 : CALL dbcsr_release(desymm_tmp)
1255 72 : CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix)
1256 : END DO
1257 24 : DEALLOCATE (dft_control%qs_control%cdft_control%matrix_p)
1258 : END IF
1259 : ! Occupation numbers
1260 472 : IF (.NOT. uniform_occupation) THEN
1261 84 : DO ispin = 1, nspins
1262 56 : IF (ncol_mo(ispin) /= SIZE(dft_control%qs_control%cdft_control%occupations(ispin)%array)) THEN
1263 0 : CPABORT("Array dimensions dont match.")
1264 : END IF
1265 56 : IF (force_env_qs%para_env%is_source()) THEN
1266 84 : ALLOCATE (occno_tmp(iforce_eval, ispin)%array(ncol_mo(ispin)))
1267 126 : occno_tmp(iforce_eval, ispin)%array = dft_control%qs_control%cdft_control%occupations(ispin)%array
1268 : END IF
1269 84 : DEALLOCATE (dft_control%qs_control%cdft_control%occupations(ispin)%array)
1270 : END DO
1271 28 : 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 100 : 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 100 : 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 630 : ALLOCATE (mixed_wmat_tmp(nforce_eval, nvar))
1283 100 : IF (mixed_cdft%calculate_metric) THEN
1284 130 : ALLOCATE (mixed_matrix_p_tmp(nforce_eval, nspins))
1285 : END IF
1286 324 : DO iforce_eval = 1, nforce_eval
1287 : ! MO coefficients
1288 660 : DO ispin = 1, nspins
1289 436 : NULLIFY (fm_struct_mo)
1290 : CALL cp_fm_struct_create(fm_struct_mo, nrow_global=nrow_mo(ispin), ncol_global=ncol_mo(ispin), &
1291 436 : 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 436 : //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 436 : mixed_cdft%blacs_env%para_env)
1299 436 : CALL cp_fm_release(mo_coeff_tmp(iforce_eval, ispin))
1300 660 : CALL cp_fm_struct_release(fm_struct_mo)
1301 : END DO
1302 : ! Weight
1303 452 : DO ivar = 1, nvar
1304 228 : NULLIFY (w_matrix(iforce_eval, ivar)%matrix)
1305 228 : 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 228 : 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 228 : mixed_cdft%blacs_env%para_env)
1312 228 : 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 228 : w_matrix(iforce_eval, ivar)%matrix)
1316 452 : CALL cp_fm_release(mixed_wmat_tmp(iforce_eval, ivar))
1317 : END DO
1318 : ! Density matrix (fm -> dbcsr)
1319 324 : IF (mixed_cdft%calculate_metric) THEN
1320 90 : DO ispin = 1, nspins
1321 60 : NULLIFY (density_matrix(iforce_eval, ispin)%matrix)
1322 60 : 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 60 : 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 60 : mixed_cdft%blacs_env%para_env)
1330 60 : 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 60 : density_matrix(iforce_eval, ispin)%matrix)
1333 90 : CALL cp_fm_release(mixed_matrix_p_tmp(iforce_eval, ispin))
1334 : END DO
1335 : END IF
1336 : END DO
1337 100 : CALL cp_fm_struct_release(fm_struct_wmat)
1338 100 : DEALLOCATE (mo_coeff_tmp, wmat_tmp, mixed_wmat_tmp)
1339 100 : IF (mixed_cdft%calculate_metric) THEN
1340 14 : DEALLOCATE (matrix_p_tmp)
1341 14 : 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 100 : name="OVERLAP_MATRIX")
1347 100 : CALL cp_fm_struct_release(fm_struct_overlap)
1348 : CALL cp_fm_copy_general(matrix_s_tmp, &
1349 : mixed_matrix_s_tmp, &
1350 100 : mixed_cdft%blacs_env%para_env)
1351 100 : CALL cp_fm_release(matrix_s_tmp)
1352 100 : CALL copy_fm_to_dbcsr_bc(mixed_matrix_s_tmp, mixed_matrix_s)
1353 100 : CALL cp_fm_release(mixed_matrix_s_tmp)
1354 : ! Occupation numbers
1355 100 : IF (.NOT. uniform_occupation) THEN
1356 42 : DO iforce_eval = 1, nforce_eval
1357 98 : DO ispin = 1, nspins
1358 168 : ALLOCATE (mixed_cdft%occupations(iforce_eval, ispin)%array(ncol_mo(ispin)))
1359 252 : mixed_cdft%occupations(iforce_eval, ispin)%array = 0.0_dp
1360 56 : IF (ASSOCIATED(occno_tmp(iforce_eval, ispin)%array)) THEN
1361 126 : mixed_cdft%occupations(iforce_eval, ispin)%array = occno_tmp(iforce_eval, ispin)%array
1362 28 : DEALLOCATE (occno_tmp(iforce_eval, ispin)%array)
1363 : END IF
1364 476 : CALL force_env%para_env%sum(mixed_cdft%occupations(iforce_eval, ispin)%array)
1365 : END DO
1366 : END DO
1367 14 : DEALLOCATE (occno_tmp)
1368 : END IF
1369 100 : DEALLOCATE (ncol_mo, nrow_mo)
1370 :
1371 300 : 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 100 : 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 100 : NULLIFY (print_section, mixed_cdft)
1389 :
1390 100 : logger => cp_get_default_logger()
1391 100 : CPASSERT(ASSOCIATED(force_env))
1392 : CALL force_env_get(force_env=force_env, &
1393 100 : force_env_section=force_env_section)
1394 100 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1395 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1396 100 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1397 100 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1398 : !
1399 100 : CPASSERT(ALLOCATED(mixed_cdft%results%strength))
1400 100 : CPASSERT(ALLOCATED(mixed_cdft%results%W_diagonal))
1401 100 : CPASSERT(ALLOCATED(mixed_cdft%results%S))
1402 100 : CPASSERT(ALLOCATED(mixed_cdft%results%energy))
1403 100 : nforce_eval = SIZE(force_env%sub_force_env)
1404 100 : nvar = SIZE(mixed_cdft%results%strength, 1)
1405 100 : npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
1406 100 : IF (iounit > 0) THEN
1407 : WRITE (iounit, '(/,T3,A,T66)') &
1408 50 : '------------------------- CDFT coupling information --------------------------'
1409 : WRITE (iounit, '(T3,A,T66,(3X,F12.2))') &
1410 50 : 'Information at step (fs):', mixed_cdft%sim_step*mixed_cdft%sim_dt
1411 141 : DO ipermutation = 1, npermutations
1412 91 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1413 91 : WRITE (iounit, '(/,T3,A)') REPEAT('#', 44)
1414 91 : WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### CDFT states I =', istate, ' and J = ', jstate, ' ######'
1415 91 : WRITE (iounit, '(T3,A)') REPEAT('#', 44)
1416 183 : DO ivar = 1, nvar
1417 92 : IF (ivar > 1) THEN
1418 1 : WRITE (iounit, '(A)') ''
1419 : END IF
1420 92 : WRITE (iounit, '(T3,A,T60,(3X,I18))') 'Atomic group:', ivar
1421 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1422 92 : 'Strength of constraint I:', mixed_cdft%results%strength(ivar, istate)
1423 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1424 92 : 'Strength of constraint J:', mixed_cdft%results%strength(ivar, jstate)
1425 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1426 92 : 'Final value of constraint I:', mixed_cdft%results%W_diagonal(ivar, istate)
1427 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1428 183 : '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 91 : 'Overlap between states I and J:', mixed_cdft%results%S(istate, jstate)
1432 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1433 91 : 'Charge transfer energy (J-I) (Hartree):', (mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate))
1434 91 : WRITE (iounit, *)
1435 91 : IF (ALLOCATED(mixed_cdft%results%rotation)) THEN
1436 89 : IF (ABS(mixed_cdft%results%rotation(ipermutation))*1.0E3_dp >= 0.1_dp) THEN
1437 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1438 87 : 'Diabatic electronic coupling (rotation, mHartree):', &
1439 174 : ABS(mixed_cdft%results%rotation(ipermutation)*1.0E3_dp)
1440 : ELSE
1441 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1442 2 : 'Diabatic electronic coupling (rotation, microHartree):', &
1443 4 : ABS(mixed_cdft%results%rotation(ipermutation)*1.0E6_dp)
1444 : END IF
1445 : END IF
1446 91 : IF (ALLOCATED(mixed_cdft%results%lowdin)) THEN
1447 10 : IF (ABS(mixed_cdft%results%lowdin(ipermutation))*1.0E3_dp >= 0.1_dp) THEN
1448 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1449 9 : 'Diabatic electronic coupling (Lowdin, mHartree):', &
1450 18 : ABS(mixed_cdft%results%lowdin(ipermutation)*1.0E3_dp)
1451 : ELSE
1452 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1453 1 : 'Diabatic electronic coupling (Lowdin, microHartree):', &
1454 2 : ABS(mixed_cdft%results%lowdin(ipermutation)*1.0E6_dp)
1455 : END IF
1456 : END IF
1457 91 : IF (ALLOCATED(mixed_cdft%results%wfn)) THEN
1458 6 : IF (mixed_cdft%results%wfn(ipermutation)*1.0E3_dp >= 0.1_dp) THEN
1459 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1460 5 : 'Diabatic electronic coupling (wfn overlap, mHartree):', &
1461 10 : ABS(mixed_cdft%results%wfn(ipermutation)*1.0E3_dp)
1462 : ELSE
1463 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1464 1 : 'Diabatic electronic coupling (wfn overlap, microHartree):', &
1465 2 : ABS(mixed_cdft%results%wfn(ipermutation)*1.0E6_dp)
1466 : END IF
1467 : END IF
1468 91 : IF (ALLOCATED(mixed_cdft%results%nonortho)) THEN
1469 : WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
1470 50 : 'Diabatic electronic coupling (nonorthogonal, Hartree):', mixed_cdft%results%nonortho(ipermutation)
1471 : END IF
1472 232 : IF (ALLOCATED(mixed_cdft%results%metric)) THEN
1473 9 : WRITE (iounit, *)
1474 9 : IF (SIZE(mixed_cdft%results%metric, 2) == 1) THEN
1475 : WRITE (iounit, '(T3,A,T66,(3X,F12.6))') &
1476 0 : 'Coupling reliability metric (0 is ideal):', mixed_cdft%results%metric(ipermutation, 1)
1477 : ELSE
1478 : WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
1479 9 : 'Coupling reliability metric (0 is ideal):', &
1480 18 : 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 50 : '------------------------------------------------------------------------------'
1486 : END IF
1487 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1488 100 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1489 :
1490 100 : 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 100 : 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 100 : NULLIFY (mixed_cdft)
1504 :
1505 100 : CPASSERT(ASSOCIATED(force_env))
1506 100 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1507 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1508 100 : CALL mixed_cdft_result_type_release(mixed_cdft%results)
1509 :
1510 100 : 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 1281 : 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 1281 : npermutations = n*(n - 1)/2 ! Size of upper triangular part
1528 1281 : IF (ipermutation > npermutations) THEN
1529 0 : CPABORT("Permutation index out of bounds")
1530 : END IF
1531 1281 : kpermutation = 0
1532 2157 : DO krow = 1, n
1533 7599 : DO kcol = krow + 1, n
1534 6723 : kpermutation = kpermutation + 1
1535 7599 : IF (kpermutation == ipermutation) THEN
1536 1281 : i = krow
1537 1281 : j = kcol
1538 1281 : 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 34 : 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 34 : n1 = SIZE(fun, 1)
1566 34 : n2 = SIZE(fun, 2)
1567 34 : n3 = SIZE(fun, 3)
1568 34 : nzeroed_total = 0
1569 34 : IF (.NOT. just_zero) THEN
1570 34 : CPASSERT(PRESENT(bounds))
1571 34 : CPASSERT(PRESENT(work))
1572 : lb = 1
1573 : lb_final = .FALSE.
1574 : ub_final = .FALSE.
1575 : END IF
1576 1586 : DO i3 = 1, n3
1577 1552 : IF (.NOT. just_zero) nzeroed = 0
1578 75920 : DO i2 = 1, n2
1579 1956496 : DO i1 = 1, n1
1580 1954944 : IF (fun(i1, i2, i3) < th) THEN
1581 983300 : IF (.NOT. just_zero) THEN
1582 983300 : nzeroed = nzeroed + 1
1583 983300 : nzeroed_total = nzeroed_total + 1
1584 : ELSE
1585 0 : fun(i1, i2, i3) = 0.0_dp
1586 : END IF
1587 : END IF
1588 : END DO
1589 : END DO
1590 1586 : IF (.NOT. just_zero) THEN
1591 1552 : IF (nzeroed == (n2*n1)) THEN
1592 80 : IF (.NOT. lb_final) THEN
1593 : lb = i3
1594 56 : ELSE IF (.NOT. ub_final) THEN
1595 8 : ub = i3
1596 8 : 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 34 : IF (.NOT. just_zero) THEN
1605 34 : IF (.NOT. ub_final) ub = n3
1606 34 : bounds(1) = lb
1607 34 : bounds(2) = ub
1608 102 : bounds = bounds - (n3/2) - 1
1609 34 : work = n3*n2*n1 - nzeroed_total
1610 : END IF
1611 :
1612 34 : 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 8 : 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 8 : 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 8 : NULLIFY (force_env_section, block_section)
1640 0 : CPASSERT(ASSOCIATED(force_env))
1641 8 : nforce_eval = SIZE(force_env%sub_force_env)
1642 :
1643 : CALL force_env_get(force_env=force_env, &
1644 8 : force_env_section=force_env_section)
1645 8 : block_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%BLOCK_DIAGONALIZE")
1646 :
1647 8 : CALL section_vals_get(block_section, explicit=explicit)
1648 8 : IF (.NOT. explicit) THEN
1649 : CALL cp_abort(__LOCATION__, &
1650 : "Block diagonalization of CDFT Hamiltonian was requested, but the "// &
1651 0 : "corresponding input section is missing!")
1652 : END IF
1653 :
1654 8 : CALL section_vals_val_get(block_section, "BLOCK", n_rep_val=nblk)
1655 44 : ALLOCATE (blocks(nblk))
1656 28 : DO i = 1, nblk
1657 20 : NULLIFY (blocks(i)%array)
1658 20 : CALL section_vals_val_get(block_section, "BLOCK", i_rep_val=i, i_vals=tmplist)
1659 20 : IF (SIZE(tmplist) < 1) THEN
1660 0 : CPABORT("Each BLOCK must contain at least 1 state.")
1661 : END IF
1662 60 : ALLOCATE (blocks(i)%array(SIZE(tmplist)))
1663 66 : blocks(i)%array(:) = tmplist(:)
1664 : END DO
1665 8 : CALL section_vals_val_get(block_section, "IGNORE_EXCITED", l_val=ignore_excited)
1666 8 : CALL section_vals_val_get(block_section, "RECURSIVE_DIAGONALIZATION", l_val=do_recursive)
1667 : ! Check that the requested states exist
1668 28 : DO i = 1, nblk
1669 66 : DO j = 1, SIZE(blocks(i)%array)
1670 58 : IF (blocks(i)%array(j) < 1 .OR. blocks(i)%array(j) > nforce_eval) THEN
1671 0 : CPABORT("Requested state does not exist.")
1672 : END IF
1673 : END DO
1674 : END DO
1675 : ! Check for duplicates
1676 8 : has_duplicates = .FALSE.
1677 28 : DO i = 1, nblk
1678 : ! Within same block
1679 58 : DO j = 1, SIZE(blocks(i)%array)
1680 76 : DO k = j + 1, SIZE(blocks(i)%array)
1681 56 : IF (blocks(i)%array(j) == blocks(i)%array(k)) has_duplicates = .TRUE.
1682 : END DO
1683 : END DO
1684 : ! Within different blocks
1685 46 : DO j = i + 1, nblk
1686 74 : DO k = 1, SIZE(blocks(i)%array)
1687 122 : DO l = 1, SIZE(blocks(j)%array)
1688 104 : 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 8 : IF (has_duplicates) CPABORT("Duplicate states are not allowed.")
1694 8 : nrecursion = 1
1695 8 : IF (do_recursive) THEN
1696 2 : IF (MODULO(nblk, 2) /= 0) THEN
1697 : CALL cp_warn(__LOCATION__, &
1698 : "Number of blocks not divisible with 2. Recursive diagonalization not possible. "// &
1699 0 : "Calculation proceeds without.")
1700 0 : nrecursion = 1
1701 : ELSE
1702 2 : nrecursion = nblk/2
1703 : END IF
1704 2 : IF (nrecursion /= 1 .AND. .NOT. ignore_excited) THEN
1705 : CALL cp_abort(__LOCATION__, &
1706 0 : "Keyword IGNORE_EXCITED must be active for recursive diagonalization.")
1707 : END IF
1708 : END IF
1709 :
1710 32 : 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 10 : 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 10 : CPASSERT(ASSOCIATED(mixed_cdft))
1732 :
1733 10 : nblk = SIZE(blocks)
1734 88 : ALLOCATE (H_block(nblk), S_block(nblk))
1735 34 : DO i = 1, nblk
1736 24 : NULLIFY (H_block(i)%array)
1737 24 : NULLIFY (S_block(i)%array)
1738 96 : ALLOCATE (H_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1739 96 : ALLOCATE (S_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1740 24 : icol = 0
1741 70 : DO j = 1, SIZE(blocks(i)%array)
1742 46 : irow = 0
1743 46 : icol = icol + 1
1744 160 : DO k = 1, SIZE(blocks(i)%array)
1745 90 : irow = irow + 1
1746 90 : H_block(i)%array(irow, icol) = mixed_cdft%results%H(blocks(i)%array(k), blocks(i)%array(j))
1747 136 : 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 170 : 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 0 : " is repulsive.")
1755 : END IF
1756 : END DO
1757 :
1758 10 : 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 10 : 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 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
1777 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat_copy, S_mat_copy
1778 :
1779 : EXTERNAL :: dsygv
1780 :
1781 10 : nblk = SIZE(blocks)
1782 54 : ALLOCATE (eigenvalues(nblk))
1783 34 : DO i = 1, nblk
1784 24 : NULLIFY (eigenvalues(i)%array)
1785 72 : ALLOCATE (eigenvalues(i)%array(SIZE(blocks(i)%array)))
1786 70 : eigenvalues(i)%array = 0.0_dp
1787 : ! Workspace query
1788 24 : ALLOCATE (work(1))
1789 24 : info = 0
1790 96 : ALLOCATE (H_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1791 72 : ALLOCATE (S_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
1792 160 : H_mat_copy(:, :) = H_block(i)%array(:, :) ! Need explicit copies because dsygv destroys original values
1793 160 : 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 24 : S_mat_copy, SIZE(blocks(i)%array), eigenvalues(i)%array, work, -1, info)
1796 24 : work_array_size = NINT(work(1))
1797 24 : DEALLOCATE (H_mat_copy, S_mat_copy)
1798 : ! Allocate work array
1799 24 : DEALLOCATE (work)
1800 72 : ALLOCATE (work(work_array_size))
1801 24 : work = 0.0_dp
1802 : ! Solve Hc = eSc
1803 24 : info = 0
1804 : CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_block(i)%array, SIZE(blocks(i)%array), &
1805 24 : S_block(i)%array, SIZE(blocks(i)%array), eigenvalues(i)%array, work, work_array_size, info)
1806 24 : IF (info /= 0) THEN
1807 0 : IF (info > SIZE(blocks(i)%array)) THEN
1808 0 : CPABORT("Matrix S is not positive definite")
1809 : ELSE
1810 0 : CPABORT("Diagonalization of H matrix failed.")
1811 : END IF
1812 : END IF
1813 34 : DEALLOCATE (work)
1814 : END DO
1815 :
1816 10 : 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 10 : 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 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, H_offdiag, S_mat, S_offdiag
1843 :
1844 : EXTERNAL :: dsygv
1845 :
1846 70 : ALLOCATE (H_mat(n, n), S_mat(n, n))
1847 10 : nblk = SIZE(blocks)
1848 10 : ignore_excited = (nblk == n)
1849 : ! The diagonal contains the eigenvalues of each block
1850 10 : IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Eigenvalues of the block diagonalized states"
1851 10 : H_mat(:, :) = 0.0_dp
1852 10 : S_mat(:, :) = 0.0_dp
1853 10 : k = 1
1854 34 : DO i = 1, nblk
1855 24 : IF (iounit > 0) WRITE (iounit, '(T6,A,I3)') "Block", i
1856 42 : DO j = 1, SIZE(eigenvalues(i)%array)
1857 28 : H_mat(k, k) = eigenvalues(i)%array(j)
1858 28 : S_mat(k, k) = 1.0_dp
1859 28 : k = k + 1
1860 28 : IF (iounit > 0) THEN
1861 14 : IF (j == 1) THEN
1862 12 : 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 2 : 'Excited state (', j - 1, ' ) energy:', eigenvalues(i)%array(j)
1866 : END IF
1867 : END IF
1868 32 : 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 10 : npermutations = nblk*(nblk - 1)/2
1873 10 : IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Interactions between block diagonalized states"
1874 30 : DO ipermutation = 1, npermutations
1875 20 : CALL map_permutation_to_states(nblk, ipermutation, i, j)
1876 : ! Get the untransformed off-diagonal block
1877 80 : ALLOCATE (H_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
1878 60 : ALLOCATE (S_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
1879 58 : icol = 0
1880 58 : DO k = 1, SIZE(blocks(j)%array)
1881 38 : irow = 0
1882 38 : icol = icol + 1
1883 134 : DO l = 1, SIZE(blocks(i)%array)
1884 76 : irow = irow + 1
1885 76 : H_offdiag(irow, icol) = mixed_cdft%results%H(blocks(i)%array(l), blocks(j)%array(k))
1886 114 : 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 134 : 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 0 : " and "//TRIM(ADJUSTL(cp_to_string(j)))//" is repulsive.")
1894 : END IF
1895 : ! Now transform: C_i^T * H * C_j
1896 838 : H_offdiag(:, :) = MATMUL(H_offdiag, H_block(j)%array)
1897 858 : H_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), H_offdiag)
1898 838 : S_offdiag(:, :) = MATMUL(S_offdiag, H_block(j)%array)
1899 858 : 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 30 : IF (ANY(S_offdiag < 0.0_dp)) THEN
1905 58 : DO l = 1, SIZE(S_offdiag, 2)
1906 134 : DO k = 1, SIZE(S_offdiag, 1)
1907 114 : IF (S_offdiag(k, l) < 0.0_dp) THEN
1908 36 : S_offdiag(k, l) = -1.0_dp*S_offdiag(k, l)
1909 36 : H_offdiag(k, l) = -1.0_dp*H_offdiag(k, l)
1910 : END IF
1911 : END DO
1912 : END DO
1913 : END IF
1914 20 : IF (ignore_excited) THEN
1915 18 : H_mat(i, j) = H_offdiag(1, 1)
1916 18 : H_mat(j, i) = H_mat(i, j)
1917 18 : S_mat(i, j) = S_offdiag(1, 1)
1918 18 : S_mat(j, i) = S_mat(i, j)
1919 : ELSE
1920 2 : irow = 1
1921 2 : icol = 1
1922 2 : DO k = 1, i - 1
1923 2 : irow = irow + SIZE(blocks(k)%array)
1924 : END DO
1925 4 : DO k = 1, j - 1
1926 4 : icol = icol + SIZE(blocks(k)%array)
1927 : END DO
1928 14 : H_mat(irow:irow + SIZE(H_offdiag, 1) - 1, icol:icol + SIZE(H_offdiag, 2) - 1) = H_offdiag(:, :)
1929 14 : H_mat(icol:icol + SIZE(H_offdiag, 2) - 1, irow:irow + SIZE(H_offdiag, 1) - 1) = TRANSPOSE(H_offdiag)
1930 14 : S_mat(irow:irow + SIZE(H_offdiag, 1) - 1, icol:icol + SIZE(H_offdiag, 2) - 1) = S_offdiag(:, :)
1931 14 : S_mat(icol:icol + SIZE(H_offdiag, 2) - 1, irow:irow + SIZE(H_offdiag, 1) - 1) = TRANSPOSE(S_offdiag)
1932 : END IF
1933 20 : IF (iounit > 0) THEN
1934 10 : WRITE (iounit, '(/,T3,A)') REPEAT('#', 39)
1935 10 : WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### Blocks I =', i, ' and J = ', j, ' ######'
1936 10 : WRITE (iounit, '(T3,A)') REPEAT('#', 39)
1937 10 : WRITE (iounit, '(T3,A)') 'Interaction energies'
1938 21 : DO irow = 1, SIZE(H_offdiag, 1)
1939 20 : ilabel = "(ground state)"
1940 20 : IF (irow > 1) THEN
1941 10 : IF (ignore_excited) EXIT
1942 1 : WRITE (tmp, '(I3)') irow - 1
1943 1 : ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
1944 : END IF
1945 34 : DO icol = 1, SIZE(H_offdiag, 2)
1946 21 : jlabel = "(ground state)"
1947 21 : IF (icol > 1) THEN
1948 10 : IF (ignore_excited) EXIT
1949 2 : WRITE (tmp, '(I3)') icol - 1
1950 2 : jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
1951 : END IF
1952 24 : WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', H_offdiag(irow, icol)
1953 : END DO
1954 : END DO
1955 10 : WRITE (iounit, '(T3,A)') 'Overlaps'
1956 21 : DO irow = 1, SIZE(H_offdiag, 1)
1957 20 : ilabel = "(ground state)"
1958 20 : IF (irow > 1) THEN
1959 10 : IF (ignore_excited) EXIT
1960 1 : ilabel = "(excited state)"
1961 1 : WRITE (tmp, '(I3)') irow - 1
1962 1 : ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
1963 : END IF
1964 34 : DO icol = 1, SIZE(H_offdiag, 2)
1965 21 : jlabel = "(ground state)"
1966 21 : IF (icol > 1) THEN
1967 10 : IF (ignore_excited) EXIT
1968 2 : WRITE (tmp, '(I3)') icol - 1
1969 2 : jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
1970 : END IF
1971 24 : 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 50 : DEALLOCATE (H_offdiag, S_offdiag)
1976 : END DO
1977 10 : CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat, S=S_mat)
1978 : ! Deallocate work
1979 10 : DEALLOCATE (H_mat, S_mat)
1980 :
1981 10 : END SUBROUTINE mixed_cdft_assemble_block_diag
1982 :
1983 : END MODULE mixed_cdft_utils
|
