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