Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods for mixed CDFT calculations
10 : !> \par History
11 : !> Separated CDFT routines from mixed_environment_utils
12 : !> \author Nico Holmberg [01.2017]
13 : ! **************************************************************************************************
14 : MODULE mixed_cdft_methods
15 : USE ao_util, ONLY: exp_radius_very_extended
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_array_utils, ONLY: cp_1d_i_p_type,&
21 : cp_1d_r_p_type,&
22 : cp_2d_r_p_type
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_init_p, dbcsr_p_type, dbcsr_release, &
26 : dbcsr_release_p, dbcsr_scale, dbcsr_type
27 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
28 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
29 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
30 : cp_fm_invert,&
31 : cp_fm_transpose
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_to_fm,&
40 : cp_fm_type,&
41 : cp_fm_write_formatted
42 : USE cp_log_handling, ONLY: cp_get_default_logger,&
43 : cp_logger_type,&
44 : cp_to_string
45 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
46 : cp_print_key_unit_nr
47 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
48 : USE cp_subsys_types, ONLY: cp_subsys_get,&
49 : cp_subsys_type
50 : USE cp_units, ONLY: cp_unit_from_cp2k
51 : USE force_env_types, ONLY: force_env_get,&
52 : force_env_type,&
53 : use_qmmm,&
54 : use_qmmmx,&
55 : use_qs_force
56 : USE grid_api, ONLY: GRID_FUNC_AB,&
57 : collocate_pgf_product
58 : USE hirshfeld_methods, ONLY: create_shape_function
59 : USE hirshfeld_types, ONLY: hirshfeld_type
60 : USE input_constants, ONLY: &
61 : becke_cutoff_element, becke_cutoff_global, cdft_alpha_constraint, cdft_beta_constraint, &
62 : cdft_charge_constraint, cdft_magnetization_constraint, mix_cdft, mixed_cdft_parallel, &
63 : mixed_cdft_parallel_nobuild, mixed_cdft_serial, outer_scf_becke_constraint
64 : USE input_section_types, ONLY: section_get_lval,&
65 : section_vals_get,&
66 : section_vals_get_subs_vals,&
67 : section_vals_type,&
68 : section_vals_val_get
69 : USE kinds, ONLY: default_path_length,&
70 : dp
71 : USE machine, ONLY: m_walltime
72 : USE mathlib, ONLY: diamat_all
73 : USE memory_utilities, ONLY: reallocate
74 : USE message_passing, ONLY: mp_request_type,&
75 : mp_testall,&
76 : mp_waitall
77 : USE mixed_cdft_types, ONLY: mixed_cdft_result_type_set,&
78 : mixed_cdft_settings_type,&
79 : mixed_cdft_type,&
80 : mixed_cdft_type_create,&
81 : mixed_cdft_work_type_release
82 : USE mixed_cdft_utils, ONLY: &
83 : hfun_zero, map_permutation_to_states, mixed_cdft_assemble_block_diag, &
84 : mixed_cdft_diagonalize_blocks, mixed_cdft_get_blocks, mixed_cdft_init_structures, &
85 : mixed_cdft_parse_settings, mixed_cdft_print_couplings, mixed_cdft_read_block_diag, &
86 : mixed_cdft_redistribute_arrays, mixed_cdft_release_work, mixed_cdft_transfer_settings
87 : USE mixed_environment_types, ONLY: get_mixed_env,&
88 : mixed_environment_type,&
89 : set_mixed_env
90 : USE parallel_gemm_api, ONLY: parallel_gemm
91 : USE particle_list_types, ONLY: particle_list_type
92 : USE particle_types, ONLY: particle_type
93 : USE pw_env_types, ONLY: pw_env_get,&
94 : pw_env_type
95 : USE pw_methods, ONLY: pw_copy,&
96 : pw_scale,&
97 : pw_zero
98 : USE pw_pool_types, ONLY: pw_pool_type
99 : USE qs_cdft_types, ONLY: cdft_control_type
100 : USE qs_energy_types, ONLY: qs_energy_type
101 : USE qs_environment_types, ONLY: get_qs_env,&
102 : qs_environment_type
103 : USE qs_kind_types, ONLY: qs_kind_type
104 : USE qs_mo_io, ONLY: read_mo_set_from_restart,&
105 : wfn_restart_file_name
106 : USE qs_mo_methods, ONLY: make_basis_simple,&
107 : make_basis_sm
108 : USE qs_mo_types, ONLY: allocate_mo_set,&
109 : deallocate_mo_set,&
110 : mo_set_type,&
111 : set_mo_set
112 : USE realspace_grid_types, ONLY: realspace_grid_type,&
113 : rs_grid_zero,&
114 : transfer_rs2pw
115 : USE util, ONLY: sort
116 : #include "./base/base_uses.f90"
117 :
118 : IMPLICIT NONE
119 :
120 : PRIVATE
121 :
122 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_methods'
123 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
124 :
125 : TYPE buffers_idx_irr
126 : INTEGER :: imap(6) = 0
127 : INTEGER, DIMENSION(:), &
128 : POINTER :: iv => null()
129 : REAL(KIND=dp), POINTER, &
130 : DIMENSION(:, :, :) :: r3 => null()
131 : REAL(KIND=dp), POINTER, &
132 : DIMENSION(:, :, :, :) :: r4 => null()
133 : END TYPE buffers_idx_irr
134 :
135 : TYPE buffers_bi
136 : LOGICAL, POINTER, DIMENSION(:) :: bv => NULL()
137 : INTEGER, POINTER, DIMENSION(:) :: iv => NULL()
138 : END TYPE buffers_bi
139 :
140 : PUBLIC :: mixed_cdft_init, &
141 : mixed_cdft_build_weight, &
142 : mixed_cdft_calculate_coupling
143 :
144 : CONTAINS
145 :
146 : ! **************************************************************************************************
147 : !> \brief Initialize a mixed CDFT calculation
148 : !> \param force_env the force_env that holds the CDFT states
149 : !> \param calculate_forces determines if forces should be calculated
150 : !> \par History
151 : !> 01.2016 created [Nico Holmberg]
152 : ! **************************************************************************************************
153 530 : SUBROUTINE mixed_cdft_init(force_env, calculate_forces)
154 : TYPE(force_env_type), POINTER :: force_env
155 : LOGICAL, INTENT(IN) :: calculate_forces
156 :
157 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_init'
158 :
159 : INTEGER :: et_freq, handle, iforce_eval, iounit, &
160 : mixing_type, nforce_eval
161 : LOGICAL :: explicit, is_parallel, is_qmmm
162 : TYPE(cp_logger_type), POINTER :: logger
163 : TYPE(cp_subsys_type), POINTER :: subsys_mix
164 : TYPE(force_env_type), POINTER :: force_env_qs
165 : TYPE(mixed_cdft_settings_type) :: settings
166 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
167 : TYPE(mixed_environment_type), POINTER :: mixed_env
168 : TYPE(particle_list_type), POINTER :: particles_mix
169 : TYPE(section_vals_type), POINTER :: force_env_section, mapping_section, &
170 : md_section, mixed_section, &
171 : print_section, root_section
172 :
173 530 : NULLIFY (subsys_mix, force_env_qs, force_env_section, print_section, &
174 530 : root_section, mixed_section, md_section, mixed_env, mixed_cdft, &
175 530 : mapping_section)
176 :
177 : NULLIFY (settings%grid_span, settings%npts, settings%cutoff, settings%rel_cutoff, &
178 : settings%spherical, settings%rs_dims, settings%odd, settings%atoms, &
179 : settings%coeffs, settings%si, settings%sr, &
180 : settings%cutoffs, settings%radii)
181 :
182 530 : is_qmmm = .FALSE.
183 1060 : logger => cp_get_default_logger()
184 530 : CPASSERT(ASSOCIATED(force_env))
185 530 : nforce_eval = SIZE(force_env%sub_force_env)
186 530 : CALL timeset(routineN, handle)
187 530 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
188 530 : mixed_env => force_env%mixed_env
189 530 : mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
190 530 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
191 530 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
192 : ! Check if a mixed CDFT calculation is requested
193 530 : CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
194 530 : IF (mixing_type == mix_cdft .AND. .NOT. ASSOCIATED(mixed_env%cdft_control)) THEN
195 78 : mixed_env%do_mixed_cdft = .TRUE.
196 156 : IF (mixed_env%do_mixed_cdft) THEN
197 : ! Sanity check
198 78 : IF (nforce_eval < 2) THEN
199 : CALL cp_abort(__LOCATION__, &
200 0 : "Mixed CDFT calculation requires at least 2 force_evals.")
201 : END IF
202 78 : mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
203 78 : CALL section_vals_get(mapping_section, explicit=explicit)
204 : ! The sub_force_envs must share the same geometrical structure
205 78 : IF (explicit) THEN
206 0 : CPABORT("Please disable section &MAPPING for mixed CDFT calculations")
207 : END IF
208 78 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%COUPLING", i_val=et_freq)
209 78 : IF (et_freq < 0) THEN
210 0 : mixed_env%do_mixed_et = .FALSE.
211 : ELSE
212 78 : mixed_env%do_mixed_et = .TRUE.
213 78 : IF (et_freq == 0) THEN
214 0 : mixed_env%et_freq = 1
215 : ELSE
216 78 : mixed_env%et_freq = et_freq
217 : END IF
218 : END IF
219 : ! Start initializing the mixed_cdft type
220 : ! First determine if the calculation is pure DFT or QMMM and find the qs force_env
221 258 : DO iforce_eval = 1, nforce_eval
222 180 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
223 296 : SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
224 : CASE (use_qs_force)
225 142 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
226 : CASE (use_qmmm)
227 12 : is_qmmm = .TRUE.
228 : ! This is really the container for QMMM
229 12 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
230 : CASE (use_qmmmx)
231 0 : CPABORT("No force mixing allowed for mixed CDFT QM/MM")
232 : CASE DEFAULT
233 : CALL cp_abort(__LOCATION__, &
234 : "Only use_qs_force and use_qmmm are "// &
235 154 : "supported for mixed_cdft_init")
236 : END SELECT
237 232 : CPASSERT(ASSOCIATED(force_env_qs))
238 : END DO
239 : ! Get infos about the mixed subsys
240 78 : IF (.NOT. is_qmmm) THEN
241 : CALL force_env_get(force_env=force_env, &
242 70 : subsys=subsys_mix)
243 : CALL cp_subsys_get(subsys=subsys_mix, &
244 70 : particles=particles_mix)
245 : ELSE
246 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
247 8 : cp_subsys=subsys_mix)
248 : CALL cp_subsys_get(subsys=subsys_mix, &
249 8 : particles=particles_mix)
250 : END IF
251 : ! Init mixed_cdft_type
252 78 : ALLOCATE (mixed_cdft)
253 78 : CALL mixed_cdft_type_create(mixed_cdft)
254 78 : mixed_cdft%first_iteration = .TRUE.
255 : ! Determine what run type to use
256 78 : IF (mixed_env%ngroups == 1) THEN
257 : ! States treated in serial, possibly copying CDFT weight function and gradients from state to state
258 52 : mixed_cdft%run_type = mixed_cdft_serial
259 26 : ELSE IF (mixed_env%ngroups == 2) THEN
260 26 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%PARALLEL_BUILD", l_val=is_parallel)
261 26 : IF (is_parallel) THEN
262 : ! Treat states in parallel, build weight function and gradients in parallel before SCF process
263 24 : mixed_cdft%run_type = mixed_cdft_parallel
264 24 : IF (.NOT. nforce_eval == 2) THEN
265 : CALL cp_abort(__LOCATION__, &
266 0 : "Parallel mode mixed CDFT calculation supports only 2 force_evals.")
267 : END IF
268 : ELSE
269 : ! Treat states in parallel, but each states builds its own weight function and gradients
270 2 : mixed_cdft%run_type = mixed_cdft_parallel_nobuild
271 : END IF
272 : ELSE
273 0 : mixed_cdft%run_type = mixed_cdft_parallel_nobuild
274 : END IF
275 : ! Store QMMM flag
276 78 : mixed_env%do_mixed_qmmm_cdft = is_qmmm
277 : ! Setup dynamic load balancing
278 78 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%DLB", l_val=mixed_cdft%dlb)
279 78 : mixed_cdft%dlb = mixed_cdft%dlb .AND. calculate_forces ! disable if forces are not needed
280 78 : mixed_cdft%dlb = mixed_cdft%dlb .AND. (mixed_cdft%run_type == mixed_cdft_parallel) ! disable if not parallel
281 78 : IF (mixed_cdft%dlb) THEN
282 40 : ALLOCATE (mixed_cdft%dlb_control)
283 4 : NULLIFY (mixed_cdft%dlb_control%weight, mixed_cdft%dlb_control%gradients, &
284 4 : mixed_cdft%dlb_control%cavity, mixed_cdft%dlb_control%target_list, &
285 4 : mixed_cdft%dlb_control%bo, mixed_cdft%dlb_control%expected_work, &
286 4 : mixed_cdft%dlb_control%prediction_error, mixed_cdft%dlb_control%sendbuff, &
287 4 : mixed_cdft%dlb_control%recvbuff, mixed_cdft%dlb_control%recv_work_repl, &
288 4 : mixed_cdft%dlb_control%recv_info)
289 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOAD_SCALE", &
290 4 : r_val=mixed_cdft%dlb_control%load_scale)
291 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%VERY_OVERLOADED", &
292 4 : r_val=mixed_cdft%dlb_control%very_overloaded)
293 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%MORE_WORK", &
294 4 : i_val=mixed_cdft%dlb_control%more_work)
295 : END IF
296 : ! Metric/Wavefunction overlap method/Lowdin orthogonalization/CDFT-CI
297 78 : mixed_cdft%calculate_metric = .FALSE.
298 78 : mixed_cdft%wfn_overlap_method = .FALSE.
299 78 : mixed_cdft%use_lowdin = .FALSE.
300 78 : mixed_cdft%do_ci = .FALSE.
301 78 : mixed_cdft%nonortho_coupling = .FALSE.
302 78 : mixed_cdft%identical_constraints = .TRUE.
303 78 : mixed_cdft%block_diagonalize = .FALSE.
304 78 : IF (mixed_env%do_mixed_et) THEN
305 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%METRIC", &
306 78 : l_val=mixed_cdft%calculate_metric)
307 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%WFN_OVERLAP", &
308 78 : l_val=mixed_cdft%wfn_overlap_method)
309 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOWDIN", &
310 78 : l_val=mixed_cdft%use_lowdin)
311 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%CI", &
312 78 : l_val=mixed_cdft%do_ci)
313 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%NONORTHOGONAL_COUPLING", &
314 78 : l_val=mixed_cdft%nonortho_coupling)
315 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%BLOCK_DIAGONALIZE", &
316 78 : l_val=mixed_cdft%block_diagonalize)
317 : END IF
318 : ! Inversion method
319 78 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%EPS_SVD", r_val=mixed_cdft%eps_svd)
320 78 : IF (mixed_cdft%eps_svd < 0.0_dp .OR. mixed_cdft%eps_svd > 1.0_dp) THEN
321 0 : CPABORT("Illegal value for EPS_SVD. Value must be between 0.0 and 1.0.")
322 : END IF
323 : ! MD related settings
324 78 : CALL force_env_get(force_env, root_section=root_section)
325 78 : md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
326 78 : CALL section_vals_val_get(md_section, "TIMESTEP", r_val=mixed_cdft%sim_dt)
327 78 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=mixed_cdft%sim_step)
328 78 : mixed_cdft%sim_step = mixed_cdft%sim_step - 1 ! to get the first step correct
329 78 : mixed_cdft%sim_dt = cp_unit_from_cp2k(mixed_cdft%sim_dt, "fs")
330 : ! Parse constraint settings from the individual force_evals and check consistency
331 : CALL mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, &
332 78 : settings, natom=SIZE(particles_mix%els))
333 : ! Transfer settings to mixed_cdft
334 78 : CALL mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
335 : ! Initilize necessary structures
336 78 : CALL mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
337 : ! Write information about the mixed CDFT calculation
338 78 : IF (iounit > 0) THEN
339 39 : WRITE (iounit, *) ""
340 : WRITE (iounit, FMT="(T2,A,T71)") &
341 39 : "MIXED_CDFT| Activating mixed CDFT calculation"
342 : WRITE (iounit, FMT="(T2,A,T71,I10)") &
343 39 : "MIXED_CDFT| Number of CDFT states: ", nforce_eval
344 51 : SELECT CASE (mixed_cdft%run_type)
345 : CASE (mixed_cdft_parallel)
346 : WRITE (iounit, FMT="(T2,A,T71)") &
347 12 : "MIXED_CDFT| CDFT states calculation mode: parallel with build"
348 : WRITE (iounit, FMT="(T2,A,T71)") &
349 12 : "MIXED_CDFT| Becke constraint is first built using all available processors"
350 : WRITE (iounit, FMT="(T2,A,T71)") &
351 12 : " and then copied to both states with their own processor groups"
352 : CASE (mixed_cdft_serial)
353 : WRITE (iounit, FMT="(T2,A,T71)") &
354 26 : "MIXED_CDFT| CDFT states calculation mode: serial"
355 26 : IF (mixed_cdft%identical_constraints) THEN
356 : WRITE (iounit, FMT="(T2,A,T71)") &
357 25 : "MIXED_CDFT| The constraints are built before the SCF procedure of the first"
358 : WRITE (iounit, FMT="(T2,A,T71)") &
359 25 : " CDFT state and subsequently copied to the other states"
360 : ELSE
361 : WRITE (iounit, FMT="(T2,A,T71)") &
362 1 : "MIXED_CDFT| The constraints are separately built for all CDFT states"
363 : END IF
364 : CASE (mixed_cdft_parallel_nobuild)
365 : WRITE (iounit, FMT="(T2,A,T71)") &
366 1 : "MIXED_CDFT| CDFT states calculation mode: parallel without build"
367 : WRITE (iounit, FMT="(T2,A,T71)") &
368 1 : "MIXED_CDFT| The constraints are separately built for all CDFT states"
369 : CASE DEFAULT
370 39 : CPABORT("Unknown mixed CDFT run type.")
371 : END SELECT
372 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
373 39 : "MIXED_CDFT| Calculating electronic coupling between states: ", mixed_env%do_mixed_et
374 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
375 39 : "MIXED_CDFT| Calculating electronic coupling reliability metric: ", mixed_cdft%calculate_metric
376 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
377 39 : "MIXED_CDFT| Configuration interaction (CDFT-CI) was requested: ", mixed_cdft%do_ci
378 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
379 39 : "MIXED_CDFT| Block diagonalizing the mixed CDFT Hamiltonian: ", mixed_cdft%block_diagonalize
380 39 : IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
381 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
382 12 : "MIXED_CDFT| Dynamic load balancing enabled: ", mixed_cdft%dlb
383 12 : IF (mixed_cdft%dlb) THEN
384 2 : WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Dynamic load balancing parameters:"
385 : WRITE (iounit, FMT="(T2,A,T71,F10.2)") &
386 2 : "MIXED_CDFT| load_scale:", mixed_cdft%dlb_control%load_scale
387 : WRITE (iounit, FMT="(T2,A,T71,F10.2)") &
388 2 : "MIXED_CDFT| very_overloaded:", mixed_cdft%dlb_control%very_overloaded
389 : WRITE (iounit, FMT="(T2,A,T71,I10)") &
390 2 : "MIXED_CDFT| more_work:", mixed_cdft%dlb_control%more_work
391 : END IF
392 : END IF
393 39 : IF (mixed_env%do_mixed_et) THEN
394 39 : IF (mixed_cdft%eps_svd == 0.0_dp) THEN
395 31 : WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with LU decomposition."
396 : ELSE
397 8 : WRITE (iounit, FMT="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with SVD decomposition."
398 8 : WRITE (iounit, FMT="(T2,A,T71,ES10.2)") "MIXED_CDFT| EPS_SVD:", mixed_cdft%eps_svd
399 : END IF
400 : END IF
401 : END IF
402 78 : CALL set_mixed_env(mixed_env, cdft_control=mixed_cdft)
403 : END IF
404 : END IF
405 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
406 530 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
407 530 : CALL timestop(handle)
408 :
409 5300 : END SUBROUTINE mixed_cdft_init
410 :
411 : ! **************************************************************************************************
412 : !> \brief Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes
413 : !> \param force_env the force_env that holds the CDFT states
414 : !> \param calculate_forces if forces should be calculated
415 : !> \param iforce_eval index of the currently active CDFT state (serial mode only)
416 : !> \par History
417 : !> 01.2017 created [Nico Holmberg]
418 : ! **************************************************************************************************
419 286 : SUBROUTINE mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
420 : TYPE(force_env_type), POINTER :: force_env
421 : LOGICAL, INTENT(IN) :: calculate_forces
422 : INTEGER, INTENT(IN), OPTIONAL :: iforce_eval
423 :
424 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
425 :
426 286 : NULLIFY (mixed_cdft)
427 286 : CPASSERT(ASSOCIATED(force_env))
428 286 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
429 286 : CPASSERT(ASSOCIATED(mixed_cdft))
430 286 : IF (.NOT. PRESENT(iforce_eval)) THEN
431 136 : SELECT CASE (mixed_cdft%run_type)
432 : CASE (mixed_cdft_parallel)
433 36 : CALL mixed_cdft_build_weight_parallel(force_env, calculate_forces)
434 : CASE (mixed_cdft_parallel_nobuild)
435 100 : CALL mixed_cdft_set_flags(force_env)
436 : CASE DEFAULT
437 : ! Do nothing
438 : END SELECT
439 : ELSE
440 334 : SELECT CASE (mixed_cdft%run_type)
441 : CASE (mixed_cdft_serial)
442 186 : CALL mixed_cdft_transfer_weight(force_env, calculate_forces, iforce_eval)
443 : CASE DEFAULT
444 : ! Do nothing
445 : END SELECT
446 : END IF
447 :
448 286 : END SUBROUTINE mixed_cdft_build_weight
449 :
450 : ! **************************************************************************************************
451 : !> \brief Build CDFT weight and gradient on 2N processors and copy it to two N processor subgroups
452 : !> \param force_env the force_env that holds the CDFT states
453 : !> \param calculate_forces if forces should be calculted
454 : !> \par History
455 : !> 01.2016 created [Nico Holmberg]
456 : ! **************************************************************************************************
457 36 : SUBROUTINE mixed_cdft_build_weight_parallel(force_env, calculate_forces)
458 : TYPE(force_env_type), POINTER :: force_env
459 : LOGICAL, INTENT(IN) :: calculate_forces
460 :
461 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_build_weight_parallel'
462 :
463 : INTEGER :: handle, handle2, i, iforce_eval, ind, INDEX(6), iounit, j, lb_min, &
464 : my_special_work, natom, nforce_eval, recv_offset, ub_max
465 : INTEGER, DIMENSION(2, 3) :: bo
466 36 : INTEGER, DIMENSION(:), POINTER :: lb, sendbuffer_i, ub
467 : REAL(KIND=dp) :: t1, t2
468 36 : TYPE(buffers_idx_irr), DIMENSION(:), POINTER :: recvbuffer
469 : TYPE(cdft_control_type), POINTER :: cdft_control, cdft_control_target
470 : TYPE(cp_logger_type), POINTER :: logger
471 : TYPE(cp_subsys_type), POINTER :: subsys_mix
472 : TYPE(dft_control_type), POINTER :: dft_control
473 : TYPE(force_env_type), POINTER :: force_env_qs
474 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
475 : TYPE(mixed_environment_type), POINTER :: mixed_env
476 36 : TYPE(mp_request_type), DIMENSION(:), POINTER :: req_total
477 : TYPE(particle_list_type), POINTER :: particles_mix
478 : TYPE(pw_env_type), POINTER :: pw_env
479 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, mixed_auxbas_pw_pool
480 : TYPE(qs_environment_type), POINTER :: qs_env
481 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
482 :
483 36 : NULLIFY (subsys_mix, force_env_qs, particles_mix, force_env_section, print_section, &
484 36 : mixed_env, mixed_cdft, pw_env, auxbas_pw_pool, mixed_auxbas_pw_pool, &
485 36 : qs_env, dft_control, sendbuffer_i, lb, ub, req_total, recvbuffer, &
486 36 : cdft_control, cdft_control_target)
487 :
488 72 : logger => cp_get_default_logger()
489 36 : CPASSERT(ASSOCIATED(force_env))
490 36 : nforce_eval = SIZE(force_env%sub_force_env)
491 36 : CALL timeset(routineN, handle)
492 36 : t1 = m_walltime()
493 : ! Get infos about the mixed subsys
494 : CALL force_env_get(force_env=force_env, &
495 : subsys=subsys_mix, &
496 36 : force_env_section=force_env_section)
497 : CALL cp_subsys_get(subsys=subsys_mix, &
498 36 : particles=particles_mix)
499 108 : DO iforce_eval = 1, nforce_eval
500 72 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
501 36 : SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
502 : CASE (use_qs_force)
503 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
504 : CASE (use_qmmm)
505 0 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
506 : CASE DEFAULT
507 : CALL cp_abort(__LOCATION__, &
508 : "Only use_qs_force and use_qmmm are "// &
509 72 : "supported for mixed_cdft_build_weight_parallel")
510 : END SELECT
511 : END DO
512 36 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
513 : CALL force_env_get(force_env=force_env_qs, &
514 : qs_env=qs_env, &
515 28 : subsys=subsys_mix)
516 : CALL cp_subsys_get(subsys=subsys_mix, &
517 28 : particles=particles_mix)
518 : ELSE
519 8 : qs_env => force_env_qs%qmmm_env%qs_env
520 8 : CALL get_qs_env(qs_env, cp_subsys=subsys_mix)
521 : CALL cp_subsys_get(subsys=subsys_mix, &
522 8 : particles=particles_mix)
523 : END IF
524 36 : mixed_env => force_env%mixed_env
525 36 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
526 36 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
527 36 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
528 36 : CPASSERT(ASSOCIATED(mixed_cdft))
529 36 : cdft_control => mixed_cdft%cdft_control
530 36 : CPASSERT(ASSOCIATED(cdft_control))
531 : ! Calculate the Becke weight function and gradient on all np processors
532 36 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=mixed_auxbas_pw_pool)
533 36 : natom = SIZE(particles_mix%els)
534 36 : CALL mixed_becke_constraint(force_env, calculate_forces)
535 : ! Start replicating the working arrays on both np/2 processor groups
536 36 : mixed_cdft%sim_step = mixed_cdft%sim_step + 1
537 36 : CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
538 36 : cdft_control_target => dft_control%qs_control%cdft_control
539 36 : CPASSERT(dft_control%qs_control%cdft)
540 36 : CPASSERT(ASSOCIATED(cdft_control_target))
541 36 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
542 360 : bo = auxbas_pw_pool%pw_grid%bounds_local
543 : ! First determine the size of the arrays along the confinement dir
544 36 : IF (mixed_cdft%is_special) THEN
545 : my_special_work = 2
546 : ELSE
547 36 : my_special_work = 1
548 : END IF
549 432 : ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)))
550 252 : ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)))
551 144 : ALLOCATE (lb(SIZE(mixed_cdft%source_list)), ub(SIZE(mixed_cdft%source_list)))
552 36 : IF (cdft_control%becke_control%cavity_confine) THEN
553 : ! Gaussian confinement => the bounds depend on the processor and need to be communicated
554 34 : ALLOCATE (sendbuffer_i(2))
555 204 : sendbuffer_i = cdft_control%becke_control%confine_bounds
556 102 : DO i = 1, SIZE(mixed_cdft%source_list)
557 68 : ALLOCATE (recvbuffer(i)%iv(2))
558 : CALL force_env%para_env%irecv(msgout=recvbuffer(i)%iv, source=mixed_cdft%source_list(i), &
559 102 : request=req_total(i))
560 : END DO
561 68 : DO i = 1, my_special_work
562 136 : DO j = 1, SIZE(mixed_cdft%dest_list)
563 68 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
564 : CALL force_env%para_env%isend(msgin=sendbuffer_i, &
565 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
566 102 : request=req_total(ind))
567 : END DO
568 : END DO
569 34 : CALL mp_waitall(req_total)
570 : ! Find the largest/smallest value of ub/lb
571 34 : DEALLOCATE (sendbuffer_i)
572 34 : lb_min = HUGE(0)
573 34 : ub_max = -HUGE(0)
574 102 : DO i = 1, SIZE(mixed_cdft%source_list)
575 68 : lb(i) = recvbuffer(i)%iv(1)
576 68 : ub(i) = recvbuffer(i)%iv(2)
577 68 : IF (lb(i) < lb_min) lb_min = lb(i)
578 68 : IF (ub(i) > ub_max) ub_max = ub(i)
579 102 : DEALLOCATE (recvbuffer(i)%iv)
580 : END DO
581 : ! Take into account the grids already communicated during dlb
582 34 : IF (mixed_cdft%dlb) THEN
583 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
584 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
585 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
586 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
587 16 : IF (LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
588 0 : < lb_min) lb_min = LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)
589 16 : IF (UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
590 8 : > ub_max) ub_max = UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)
591 : END DO
592 : END IF
593 : END DO
594 : END IF
595 : END IF
596 : ELSE
597 : ! No confinement
598 2 : ub_max = bo(2, 3)
599 2 : lb_min = bo(1, 3)
600 6 : lb = lb_min
601 6 : ub = ub_max
602 : END IF
603 : ! Determine the sender specific indices of grid slices that are to be received
604 36 : CALL timeset(routineN//"_comm", handle2)
605 108 : DO j = 1, SIZE(recvbuffer)
606 72 : ind = j + (j/2)
607 108 : IF (mixed_cdft%is_special) THEN
608 : recvbuffer(j)%imap = [mixed_cdft%source_list_bo(1, j), mixed_cdft%source_list_bo(2, j), &
609 : mixed_cdft%source_list_bo(3, j), mixed_cdft%source_list_bo(4, j), &
610 0 : lb(j), ub(j)]
611 72 : ELSE IF (mixed_cdft%is_pencil) THEN
612 0 : recvbuffer(j)%imap = [bo(1, 1), bo(2, 1), mixed_cdft%recv_bo(ind), mixed_cdft%recv_bo(ind + 1), lb(j), ub(j)]
613 : ELSE
614 504 : recvbuffer(j)%imap = [mixed_cdft%recv_bo(ind), mixed_cdft%recv_bo(ind + 1), bo(1, 2), bo(2, 2), lb(j), ub(j)]
615 : END IF
616 : END DO
617 36 : IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_special) THEN
618 8 : IF (mixed_cdft%dlb_control%recv_work_repl(1) .OR. mixed_cdft%dlb_control%recv_work_repl(2)) THEN
619 24 : DO j = 1, 2
620 16 : recv_offset = 0
621 16 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
622 16 : recv_offset = SUM(mixed_cdft%dlb_control%recv_info(j)%target_list(2, :))
623 : END IF
624 24 : IF (mixed_cdft%is_pencil) THEN
625 0 : recvbuffer(j)%imap(1) = recvbuffer(j)%imap(1) + recv_offset
626 : ELSE
627 16 : recvbuffer(j)%imap(3) = recvbuffer(j)%imap(3) + recv_offset
628 : END IF
629 : END DO
630 : END IF
631 : END IF
632 : ! Transfer the arrays one-by-one and deallocate shared storage
633 : ! Start with the weight function
634 108 : DO j = 1, SIZE(mixed_cdft%source_list)
635 : ALLOCATE (recvbuffer(j)%r3(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
636 : recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
637 360 : recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)))
638 :
639 : CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r3, source=mixed_cdft%source_list(j), &
640 108 : request=req_total(j))
641 : END DO
642 72 : DO i = 1, my_special_work
643 144 : DO j = 1, SIZE(mixed_cdft%dest_list)
644 72 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
645 108 : IF (mixed_cdft%is_special) THEN
646 : CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%weight, &
647 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
648 0 : request=req_total(ind))
649 : ELSE
650 : CALL force_env%para_env%isend(msgin=mixed_cdft%weight, dest=mixed_cdft%dest_list(j), &
651 72 : request=req_total(ind))
652 : END IF
653 : END DO
654 : END DO
655 36 : CALL mp_waitall(req_total)
656 36 : IF (mixed_cdft%is_special) THEN
657 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
658 0 : DEALLOCATE (mixed_cdft%sendbuff(j)%weight)
659 : END DO
660 : ELSE
661 36 : DEALLOCATE (mixed_cdft%weight)
662 : END IF
663 : ! In principle, we could reduce the memory footprint of becke_pot by only transferring the nonzero portion
664 : ! of the potential, but this would require a custom integrate_v_rspace
665 36 : ALLOCATE (cdft_control_target%group(1)%weight)
666 36 : CALL auxbas_pw_pool%create_pw(cdft_control_target%group(1)%weight)
667 36 : CALL pw_zero(cdft_control_target%group(1)%weight)
668 : ! Assemble the recved slices
669 108 : DO j = 1, SIZE(mixed_cdft%source_list)
670 : cdft_control_target%group(1)%weight%array(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
671 : recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
672 7517624 : recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
673 : END DO
674 : ! Do the same for slices sent during dlb
675 36 : IF (mixed_cdft%dlb) THEN
676 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
677 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
678 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
679 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
680 : index = [LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 1), &
681 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 1), &
682 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 2), &
683 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 2), &
684 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3), &
685 104 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)]
686 : cdft_control_target%group(1)%weight%array(INDEX(1):INDEX(2), &
687 : INDEX(3):INDEX(4), &
688 : INDEX(5):INDEX(6)) = &
689 14096 : mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight
690 16 : DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight)
691 : END DO
692 : END IF
693 : END DO
694 : END IF
695 : END IF
696 : ! Gaussian confinement cavity
697 36 : IF (cdft_control%becke_control%cavity_confine) THEN
698 102 : DO j = 1, SIZE(mixed_cdft%source_list)
699 : CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r3, source=mixed_cdft%source_list(j), &
700 102 : request=req_total(j))
701 : END DO
702 68 : DO i = 1, my_special_work
703 136 : DO j = 1, SIZE(mixed_cdft%dest_list)
704 68 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
705 102 : IF (mixed_cdft%is_special) THEN
706 : CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%cavity, &
707 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
708 0 : request=req_total(ind))
709 : ELSE
710 : CALL force_env%para_env%isend(msgin=mixed_cdft%cavity, dest=mixed_cdft%dest_list(j), &
711 68 : request=req_total(ind))
712 : END IF
713 : END DO
714 : END DO
715 34 : CALL mp_waitall(req_total)
716 34 : IF (mixed_cdft%is_special) THEN
717 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
718 0 : DEALLOCATE (mixed_cdft%sendbuff(j)%cavity)
719 : END DO
720 : ELSE
721 34 : DEALLOCATE (mixed_cdft%cavity)
722 : END IF
723 : ! We only need the nonzero part of the confinement cavity
724 : ALLOCATE (cdft_control_target%becke_control%cavity_mat(bo(1, 1):bo(2, 1), &
725 : bo(1, 2):bo(2, 2), &
726 170 : lb_min:ub_max))
727 3504226 : cdft_control_target%becke_control%cavity_mat = 0.0_dp
728 :
729 102 : DO j = 1, SIZE(mixed_cdft%source_list)
730 : cdft_control_target%becke_control%cavity_mat(recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
731 : recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
732 7136554 : recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
733 : END DO
734 34 : IF (mixed_cdft%dlb) THEN
735 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
736 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
737 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
738 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
739 : index = [LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 1), &
740 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 1), &
741 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 2), &
742 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 2), &
743 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 3), &
744 104 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, 3)]
745 : cdft_control_target%becke_control%cavity_mat(INDEX(1):INDEX(2), &
746 : INDEX(3):INDEX(4), &
747 : INDEX(5):INDEX(6)) = &
748 14096 : mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity
749 16 : DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity)
750 : END DO
751 : END IF
752 : END DO
753 : END IF
754 : END IF
755 : END IF
756 108 : DO j = 1, SIZE(mixed_cdft%source_list)
757 108 : DEALLOCATE (recvbuffer(j)%r3)
758 : END DO
759 36 : IF (calculate_forces) THEN
760 : ! Gradients of the weight function
761 72 : DO j = 1, SIZE(mixed_cdft%source_list)
762 : ALLOCATE (recvbuffer(j)%r4(3*natom, recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
763 : recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
764 288 : recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)))
765 : CALL force_env%para_env%irecv(msgout=recvbuffer(j)%r4, source=mixed_cdft%source_list(j), &
766 72 : request=req_total(j))
767 : END DO
768 48 : DO i = 1, my_special_work
769 96 : DO j = 1, SIZE(mixed_cdft%dest_list)
770 48 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
771 72 : IF (mixed_cdft%is_special) THEN
772 : CALL force_env%para_env%isend(msgin=mixed_cdft%sendbuff(j)%gradients, &
773 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
774 0 : request=req_total(ind))
775 : ELSE
776 : CALL force_env%para_env%isend(msgin=cdft_control%group(1)%gradients, dest=mixed_cdft%dest_list(j), &
777 48 : request=req_total(ind))
778 : END IF
779 : END DO
780 : END DO
781 24 : CALL mp_waitall(req_total)
782 24 : IF (mixed_cdft%is_special) THEN
783 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
784 0 : DEALLOCATE (mixed_cdft%sendbuff(j)%gradients)
785 : END DO
786 0 : DEALLOCATE (mixed_cdft%sendbuff)
787 : ELSE
788 24 : DEALLOCATE (cdft_control%group(1)%gradients)
789 : END IF
790 : ALLOCATE (cdft_control_target%group(1)%gradients(3*natom, bo(1, 1):bo(2, 1), &
791 144 : bo(1, 2):bo(2, 2), lb_min:ub_max))
792 72 : DO j = 1, SIZE(mixed_cdft%source_list)
793 : cdft_control_target%group(1)%gradients(:, recvbuffer(j)%imap(1):recvbuffer(j)%imap(2), &
794 : recvbuffer(j)%imap(3):recvbuffer(j)%imap(4), &
795 39235744 : recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r4
796 72 : DEALLOCATE (recvbuffer(j)%r4)
797 : END DO
798 24 : IF (mixed_cdft%dlb) THEN
799 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
800 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
801 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
802 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
803 : index = [LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 2), &
804 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 2), &
805 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 3), &
806 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 3), &
807 : LBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 4), &
808 104 : UBOUND(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, 4)]
809 : cdft_control_target%group(1)%gradients(:, INDEX(1):INDEX(2), &
810 : INDEX(3):INDEX(4), &
811 : INDEX(5):INDEX(6)) = &
812 90896 : mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients
813 16 : DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients)
814 : END DO
815 : END IF
816 : END DO
817 : END IF
818 : END IF
819 : END IF
820 : ! Clean up remaining temporaries
821 36 : IF (mixed_cdft%dlb) THEN
822 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
823 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
824 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
825 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%recv_info(j)%target_list)) THEN
826 8 : DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%target_list)
827 : END IF
828 8 : DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs)
829 : END IF
830 : END DO
831 8 : DEALLOCATE (mixed_cdft%dlb_control%recv_info, mixed_cdft%dlb_control%recvbuff)
832 : END IF
833 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%target_list)) THEN
834 4 : DEALLOCATE (mixed_cdft%dlb_control%target_list)
835 : END IF
836 8 : DEALLOCATE (mixed_cdft%dlb_control%recv_work_repl)
837 : END IF
838 36 : DEALLOCATE (recvbuffer)
839 36 : DEALLOCATE (req_total)
840 36 : DEALLOCATE (lb)
841 36 : DEALLOCATE (ub)
842 36 : CALL timestop(handle2)
843 : ! Set some flags so the weight is not rebuilt during SCF
844 36 : cdft_control_target%external_control = .TRUE.
845 36 : cdft_control_target%need_pot = .FALSE.
846 36 : cdft_control_target%transfer_pot = .FALSE.
847 : ! Store the bound indices for force calculation
848 36 : IF (calculate_forces) THEN
849 24 : cdft_control_target%becke_control%confine_bounds(2) = ub_max
850 24 : cdft_control_target%becke_control%confine_bounds(1) = lb_min
851 : END IF
852 : CALL pw_scale(cdft_control_target%group(1)%weight, &
853 36 : cdft_control_target%group(1)%weight%pw_grid%dvol)
854 : ! Set flags for ET coupling calculation
855 36 : IF (mixed_env%do_mixed_et) THEN
856 36 : IF (MODULO(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
857 36 : dft_control%qs_control%cdft_control%do_et = .TRUE.
858 36 : dft_control%qs_control%cdft_control%calculate_metric = mixed_cdft%calculate_metric
859 : ELSE
860 0 : dft_control%qs_control%cdft_control%do_et = .FALSE.
861 0 : dft_control%qs_control%cdft_control%calculate_metric = .FALSE.
862 : END IF
863 : END IF
864 36 : t2 = m_walltime()
865 36 : IF (iounit > 0) THEN
866 18 : WRITE (iounit, '(A)') ' '
867 18 : WRITE (iounit, '(T2,A,F6.1,A)') 'MIXED_CDFT| Becke constraint built in ', t2 - t1, ' seconds'
868 18 : WRITE (iounit, '(A)') ' '
869 : END IF
870 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
871 36 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
872 36 : CALL timestop(handle)
873 :
874 108 : END SUBROUTINE mixed_cdft_build_weight_parallel
875 :
876 : ! **************************************************************************************************
877 : !> \brief Transfer CDFT weight/gradient between force_evals
878 : !> \param force_env the force_env that holds the CDFT sub_force_envs
879 : !> \param calculate_forces if forces should be computed
880 : !> \param iforce_eval index of the currently active CDFT state
881 : !> \par History
882 : !> 01.2017 created [Nico Holmberg]
883 : ! **************************************************************************************************
884 296 : SUBROUTINE mixed_cdft_transfer_weight(force_env, calculate_forces, iforce_eval)
885 : TYPE(force_env_type), POINTER :: force_env
886 : LOGICAL, INTENT(IN) :: calculate_forces
887 : INTEGER, INTENT(IN) :: iforce_eval
888 :
889 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_transfer_weight'
890 :
891 : INTEGER :: bounds_of(8), handle, iatom, igroup, &
892 : jforce_eval, nforce_eval
893 : LOGICAL, SAVE :: first_call = .TRUE.
894 : TYPE(cdft_control_type), POINTER :: cdft_control_source, cdft_control_target
895 : TYPE(dft_control_type), POINTER :: dft_control_source, dft_control_target
896 : TYPE(force_env_type), POINTER :: force_env_qs_source, force_env_qs_target
897 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
898 : TYPE(mixed_environment_type), POINTER :: mixed_env
899 : TYPE(pw_env_type), POINTER :: pw_env_source, pw_env_target
900 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool_source, &
901 : auxbas_pw_pool_target
902 : TYPE(qs_environment_type), POINTER :: qs_env_source, qs_env_target
903 :
904 148 : NULLIFY (mixed_cdft, dft_control_source, dft_control_target, force_env_qs_source, &
905 148 : force_env_qs_target, pw_env_source, pw_env_target, auxbas_pw_pool_source, &
906 148 : auxbas_pw_pool_target, qs_env_source, qs_env_target, mixed_env, &
907 148 : cdft_control_source, cdft_control_target)
908 148 : mixed_env => force_env%mixed_env
909 148 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
910 148 : CALL timeset(routineN, handle)
911 148 : IF (iforce_eval == 1) THEN
912 : jforce_eval = 1
913 : ELSE
914 86 : jforce_eval = iforce_eval - 1
915 : END IF
916 148 : nforce_eval = SIZE(force_env%sub_force_env)
917 296 : SELECT CASE (force_env%sub_force_env(jforce_eval)%force_env%in_use)
918 : CASE (use_qs_force, use_qmmm)
919 148 : force_env_qs_source => force_env%sub_force_env(jforce_eval)%force_env
920 148 : force_env_qs_target => force_env%sub_force_env(iforce_eval)%force_env
921 : CASE DEFAULT
922 : CALL cp_abort(__LOCATION__, &
923 : "Only use_qs_force and use_qmmm are "// &
924 148 : "supported for mixed_cdft_transfer_weight")
925 : END SELECT
926 148 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
927 : CALL force_env_get(force_env=force_env_qs_source, &
928 132 : qs_env=qs_env_source)
929 : CALL force_env_get(force_env=force_env_qs_target, &
930 132 : qs_env=qs_env_target)
931 : ELSE
932 16 : qs_env_source => force_env_qs_source%qmmm_env%qs_env
933 16 : qs_env_target => force_env_qs_target%qmmm_env%qs_env
934 : END IF
935 148 : IF (iforce_eval == 1) THEN
936 : ! The first force_eval builds the weight function and gradients in qs_cdft_methods.F
937 : ! Set some flags so the constraint is saved if the constraint definitions are identical in all CDFT states
938 62 : CALL get_qs_env(qs_env_source, dft_control=dft_control_source)
939 62 : cdft_control_source => dft_control_source%qs_control%cdft_control
940 62 : cdft_control_source%external_control = .FALSE.
941 62 : cdft_control_source%need_pot = .TRUE.
942 62 : IF (mixed_cdft%identical_constraints) THEN
943 60 : cdft_control_source%transfer_pot = .TRUE.
944 : ELSE
945 2 : cdft_control_source%transfer_pot = .FALSE.
946 : END IF
947 62 : mixed_cdft%sim_step = mixed_cdft%sim_step + 1
948 : ELSE
949 : ! Transfer the constraint from the ith force_eval to the i+1th
950 : CALL get_qs_env(qs_env_source, dft_control=dft_control_source, &
951 86 : pw_env=pw_env_source)
952 86 : CALL pw_env_get(pw_env_source, auxbas_pw_pool=auxbas_pw_pool_source)
953 86 : cdft_control_source => dft_control_source%qs_control%cdft_control
954 : CALL get_qs_env(qs_env_target, dft_control=dft_control_target, &
955 86 : pw_env=pw_env_target)
956 86 : CALL pw_env_get(pw_env_target, auxbas_pw_pool=auxbas_pw_pool_target)
957 86 : cdft_control_target => dft_control_target%qs_control%cdft_control
958 : ! The constraint can be transferred only when the constraint defitions are identical in all CDFT states
959 86 : IF (mixed_cdft%identical_constraints) THEN
960 : ! Weight function
961 170 : DO igroup = 1, SIZE(cdft_control_target%group)
962 86 : ALLOCATE (cdft_control_target%group(igroup)%weight)
963 86 : CALL auxbas_pw_pool_target%create_pw(cdft_control_target%group(igroup)%weight)
964 : ! We have ensured that the grids are consistent => no danger in using explicit copy
965 86 : CALL pw_copy(cdft_control_source%group(igroup)%weight, cdft_control_target%group(igroup)%weight)
966 86 : CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%group(igroup)%weight)
967 170 : DEALLOCATE (cdft_control_source%group(igroup)%weight)
968 : END DO
969 : ! Cavity
970 84 : IF (cdft_control_source%type == outer_scf_becke_constraint) THEN
971 78 : IF (cdft_control_source%becke_control%cavity_confine) THEN
972 72 : CALL auxbas_pw_pool_target%create_pw(cdft_control_target%becke_control%cavity)
973 72 : CALL pw_copy(cdft_control_source%becke_control%cavity, cdft_control_target%becke_control%cavity)
974 72 : CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%becke_control%cavity)
975 : END IF
976 : END IF
977 : ! Gradients
978 84 : IF (calculate_forces) THEN
979 40 : DO igroup = 1, SIZE(cdft_control_source%group)
980 : bounds_of = [LBOUND(cdft_control_source%group(igroup)%gradients, 1), &
981 : UBOUND(cdft_control_source%group(igroup)%gradients, 1), &
982 : LBOUND(cdft_control_source%group(igroup)%gradients, 2), &
983 : UBOUND(cdft_control_source%group(igroup)%gradients, 2), &
984 : LBOUND(cdft_control_source%group(igroup)%gradients, 3), &
985 : UBOUND(cdft_control_source%group(igroup)%gradients, 3), &
986 : LBOUND(cdft_control_source%group(igroup)%gradients, 4), &
987 340 : UBOUND(cdft_control_source%group(igroup)%gradients, 4)]
988 : ALLOCATE (cdft_control_target%group(igroup)% &
989 : gradients(bounds_of(1):bounds_of(2), bounds_of(3):bounds_of(4), &
990 120 : bounds_of(5):bounds_of(6), bounds_of(7):bounds_of(8)))
991 17857864 : cdft_control_target%group(igroup)%gradients = cdft_control_source%group(igroup)%gradients
992 40 : DEALLOCATE (cdft_control_source%group(igroup)%gradients)
993 : END DO
994 : END IF
995 : ! Atomic weight functions needed for CDFT charges
996 84 : IF (cdft_control_source%atomic_charges) THEN
997 18 : IF (.NOT. ASSOCIATED(cdft_control_target%charge)) THEN
998 10 : ALLOCATE (cdft_control_target%charge(cdft_control_target%natoms))
999 : END IF
1000 54 : DO iatom = 1, cdft_control_target%natoms
1001 36 : CALL auxbas_pw_pool_target%create_pw(cdft_control_target%charge(iatom))
1002 36 : CALL pw_copy(cdft_control_source%charge(iatom), cdft_control_target%charge(iatom))
1003 54 : CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%charge(iatom))
1004 : END DO
1005 : END IF
1006 : ! Set some flags so the weight is not rebuilt during SCF
1007 84 : cdft_control_target%external_control = .FALSE.
1008 84 : cdft_control_target%need_pot = .FALSE.
1009 : ! For states i+1 < nforce_eval, prevent deallocation of constraint
1010 84 : IF (iforce_eval == nforce_eval) THEN
1011 60 : cdft_control_target%transfer_pot = .FALSE.
1012 : ELSE
1013 24 : cdft_control_target%transfer_pot = .TRUE.
1014 : END IF
1015 84 : cdft_control_target%first_iteration = .FALSE.
1016 : ELSE
1017 : ! Force rebuild of constraint and dont save it
1018 2 : cdft_control_target%external_control = .FALSE.
1019 2 : cdft_control_target%need_pot = .TRUE.
1020 2 : cdft_control_target%transfer_pot = .FALSE.
1021 2 : IF (first_call) THEN
1022 2 : cdft_control_target%first_iteration = .TRUE.
1023 : ELSE
1024 0 : cdft_control_target%first_iteration = .FALSE.
1025 : END IF
1026 : END IF
1027 : END IF
1028 : ! Set flags for ET coupling calculation
1029 148 : IF (mixed_env%do_mixed_et) THEN
1030 148 : IF (MODULO(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
1031 148 : IF (iforce_eval == 1) THEN
1032 62 : cdft_control_source%do_et = .TRUE.
1033 62 : cdft_control_source%calculate_metric = mixed_cdft%calculate_metric
1034 : ELSE
1035 86 : cdft_control_target%do_et = .TRUE.
1036 86 : cdft_control_target%calculate_metric = mixed_cdft%calculate_metric
1037 : END IF
1038 : ELSE
1039 0 : IF (iforce_eval == 1) THEN
1040 0 : cdft_control_source%do_et = .FALSE.
1041 0 : cdft_control_source%calculate_metric = .FALSE.
1042 : ELSE
1043 0 : cdft_control_target%do_et = .FALSE.
1044 0 : cdft_control_target%calculate_metric = .FALSE.
1045 : END IF
1046 : END IF
1047 : END IF
1048 148 : IF (iforce_eval == nforce_eval .AND. first_call) first_call = .FALSE.
1049 148 : CALL timestop(handle)
1050 :
1051 148 : END SUBROUTINE mixed_cdft_transfer_weight
1052 :
1053 : ! **************************************************************************************************
1054 : !> \brief In case CDFT states are treated in parallel, sets flags so that each CDFT state
1055 : !> builds their own weight functions and gradients
1056 : !> \param force_env the force_env that holds the CDFT sub_force_envs
1057 : !> \par History
1058 : !> 09.2018 created [Nico Holmberg]
1059 : ! **************************************************************************************************
1060 4 : SUBROUTINE mixed_cdft_set_flags(force_env)
1061 : TYPE(force_env_type), POINTER :: force_env
1062 :
1063 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_set_flags'
1064 :
1065 : INTEGER :: handle, iforce_eval, nforce_eval
1066 : LOGICAL, SAVE :: first_call = .TRUE.
1067 : TYPE(cdft_control_type), POINTER :: cdft_control
1068 : TYPE(dft_control_type), POINTER :: dft_control
1069 : TYPE(force_env_type), POINTER :: force_env_qs
1070 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1071 : TYPE(mixed_environment_type), POINTER :: mixed_env
1072 : TYPE(qs_environment_type), POINTER :: qs_env
1073 :
1074 2 : NULLIFY (mixed_cdft, dft_control, force_env_qs, qs_env, mixed_env, cdft_control)
1075 2 : mixed_env => force_env%mixed_env
1076 2 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1077 2 : CALL timeset(routineN, handle)
1078 2 : nforce_eval = SIZE(force_env%sub_force_env)
1079 6 : DO iforce_eval = 1, nforce_eval
1080 4 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
1081 2 : SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
1082 : CASE (use_qs_force, use_qmmm)
1083 0 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
1084 : CASE DEFAULT
1085 : CALL cp_abort(__LOCATION__, &
1086 : "Only use_qs_force and use_qmmm are "// &
1087 2 : "supported for mixed_cdft_set_flags")
1088 : END SELECT
1089 2 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1090 2 : CALL force_env_get(force_env=force_env_qs, qs_env=qs_env)
1091 : ELSE
1092 0 : qs_env => force_env_qs%qmmm_env%qs_env
1093 : END IF
1094 : ! All force_evals build the weight function and gradients in qs_cdft_methods.F
1095 : ! Update flags to match run type
1096 2 : CALL get_qs_env(qs_env, dft_control=dft_control)
1097 2 : cdft_control => dft_control%qs_control%cdft_control
1098 2 : cdft_control%external_control = .FALSE.
1099 2 : cdft_control%need_pot = .TRUE.
1100 2 : cdft_control%transfer_pot = .FALSE.
1101 2 : IF (first_call) THEN
1102 2 : cdft_control%first_iteration = .TRUE.
1103 : ELSE
1104 0 : cdft_control%first_iteration = .FALSE.
1105 : END IF
1106 : ! Set flags for ET coupling calculation
1107 4 : IF (mixed_env%do_mixed_et) THEN
1108 2 : IF (MODULO(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
1109 2 : cdft_control%do_et = .TRUE.
1110 2 : cdft_control%calculate_metric = mixed_cdft%calculate_metric
1111 : ELSE
1112 0 : cdft_control%do_et = .FALSE.
1113 0 : cdft_control%calculate_metric = .FALSE.
1114 : END IF
1115 : END IF
1116 : END DO
1117 2 : mixed_cdft%sim_step = mixed_cdft%sim_step + 1
1118 2 : IF (first_call) first_call = .FALSE.
1119 2 : CALL timestop(handle)
1120 :
1121 2 : END SUBROUTINE mixed_cdft_set_flags
1122 :
1123 : ! **************************************************************************************************
1124 : !> \brief Driver routine to calculate the electronic coupling(s) between CDFT states.
1125 : !> \param force_env the force_env that holds the CDFT states
1126 : !> \par History
1127 : !> 02.15 created [Nico Holmberg]
1128 : ! **************************************************************************************************
1129 200 : SUBROUTINE mixed_cdft_calculate_coupling(force_env)
1130 : TYPE(force_env_type), POINTER :: force_env
1131 :
1132 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_calculate_coupling'
1133 :
1134 : INTEGER :: handle
1135 :
1136 100 : CPASSERT(ASSOCIATED(force_env))
1137 100 : CALL timeset(routineN, handle)
1138 : ! Move needed arrays from individual CDFT states to the mixed CDFT env
1139 100 : CALL mixed_cdft_redistribute_arrays(force_env)
1140 : ! Calculate the mixed CDFT Hamiltonian and overlap matrices.
1141 : ! All work matrices defined in the wavefunction basis get deallocated on exit.
1142 : ! Any analyses which depend on these work matrices are performed within.
1143 100 : CALL mixed_cdft_interaction_matrices(force_env)
1144 : ! Calculate eletronic couplings between states (Lowdin/rotation)
1145 100 : CALL mixed_cdft_calculate_coupling_low(force_env)
1146 : ! Print out couplings
1147 100 : CALL mixed_cdft_print_couplings(force_env)
1148 : ! Block diagonalize the mixed CDFT Hamiltonian matrix
1149 100 : CALL mixed_cdft_block_diag(force_env)
1150 : ! CDFT Configuration Interaction
1151 100 : CALL mixed_cdft_configuration_interaction(force_env)
1152 : ! Clean up
1153 100 : CALL mixed_cdft_release_work(force_env)
1154 100 : CALL timestop(handle)
1155 :
1156 100 : END SUBROUTINE mixed_cdft_calculate_coupling
1157 :
1158 : ! **************************************************************************************************
1159 : !> \brief Routine to calculate the mixed CDFT Hamiltonian and overlap matrices.
1160 : !> \param force_env the force_env that holds the CDFT states
1161 : !> \par History
1162 : !> 11.17 created [Nico Holmberg]
1163 : ! **************************************************************************************************
1164 100 : SUBROUTINE mixed_cdft_interaction_matrices(force_env)
1165 : TYPE(force_env_type), POINTER :: force_env
1166 :
1167 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_interaction_matrices'
1168 :
1169 : INTEGER :: check_ao(2), check_mo(2), handle, iforce_eval, ipermutation, ispin, istate, ivar, &
1170 : j, jstate, k, moeigvalunit, mounit, nao, ncol_local, nforce_eval, nmo, npermutations, &
1171 : nrow_local, nspins, nvar
1172 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ncol_mo, nrow_mo
1173 100 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: homo
1174 : LOGICAL :: nelectron_mismatch, print_mo, &
1175 : print_mo_eigval, should_scale, &
1176 : uniform_occupation
1177 : REAL(KIND=dp) :: c(2), eps_occupied, nelectron_tot, &
1178 : sum_a(2), sum_b(2)
1179 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_nonortho, eigenv, energy, Sda
1180 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, S_det, S_mat, strength, tmp_mat, &
1181 100 : W_diagonal, Wad, Wda
1182 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: a, b
1183 100 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigval
1184 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo, mo_mo_fmstruct
1185 : TYPE(cp_fm_type) :: inverse_mat, Tinverse, tmp2
1186 100 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_overlap
1187 100 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: w_matrix_mo
1188 100 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
1189 : TYPE(cp_logger_type), POINTER :: logger
1190 100 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix, density_matrix_diff, &
1191 100 : w_matrix
1192 : TYPE(dbcsr_type), POINTER :: mixed_matrix_s
1193 : TYPE(dft_control_type), POINTER :: dft_control
1194 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1195 : TYPE(mixed_environment_type), POINTER :: mixed_env
1196 : TYPE(qs_energy_type), POINTER :: energy_qs
1197 : TYPE(qs_environment_type), POINTER :: qs_env
1198 : TYPE(section_vals_type), POINTER :: force_env_section, mixed_cdft_section, &
1199 : print_section
1200 :
1201 100 : NULLIFY (force_env_section, print_section, mixed_cdft_section, &
1202 100 : mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
1203 100 : density_matrix_diff, mo_mo_fmstruct, &
1204 100 : mixed_mo_coeff, mixed_matrix_s, &
1205 100 : density_matrix, energy_qs, w_matrix, mo_eigval)
1206 200 : logger => cp_get_default_logger()
1207 100 : CPASSERT(ASSOCIATED(force_env))
1208 100 : CALL timeset(routineN, handle)
1209 : CALL force_env_get(force_env=force_env, &
1210 100 : force_env_section=force_env_section)
1211 100 : mixed_env => force_env%mixed_env
1212 100 : nforce_eval = SIZE(force_env%sub_force_env)
1213 100 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1214 100 : IF (section_get_lval(print_section, "MO_OVERLAP_MATRIX")) THEN
1215 2 : print_mo = .TRUE.
1216 2 : mounit = cp_print_key_unit_nr(logger, print_section, extension='.moOverlap', on_file=.TRUE.)
1217 : ELSE
1218 : print_mo = .FALSE.
1219 : END IF
1220 100 : IF (section_get_lval(print_section, "MO_OVERLAP_EIGENVALUES")) THEN
1221 14 : print_mo_eigval = .TRUE.
1222 14 : moeigvalunit = cp_print_key_unit_nr(logger, print_section, extension='.moOverlapEigval', on_file=.TRUE.)
1223 : ELSE
1224 : print_mo_eigval = .FALSE.
1225 : END IF
1226 100 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1227 : ! Get redistributed work matrices
1228 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1229 100 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%mixed_mo_coeff))
1230 100 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%w_matrix))
1231 100 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%mixed_matrix_s))
1232 100 : mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
1233 100 : w_matrix => mixed_cdft%matrix%w_matrix
1234 100 : mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
1235 100 : IF (mixed_cdft%calculate_metric) THEN
1236 14 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%density_matrix))
1237 14 : density_matrix => mixed_cdft%matrix%density_matrix
1238 : END IF
1239 : ! Get number of weight functions per state
1240 100 : nvar = SIZE(w_matrix, 2)
1241 100 : nspins = SIZE(mixed_mo_coeff, 2)
1242 : ! Check that the number of MOs/AOs is equal in every CDFT state
1243 400 : ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
1244 294 : DO ispin = 1, nspins
1245 194 : CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=check_mo(1), nrow_global=check_ao(1))
1246 536 : DO iforce_eval = 2, nforce_eval
1247 242 : CALL cp_fm_get_info(mixed_mo_coeff(iforce_eval, ispin), ncol_global=check_mo(2), nrow_global=check_ao(2))
1248 242 : IF (check_ao(1) /= check_ao(2)) THEN
1249 : CALL cp_abort(__LOCATION__, &
1250 0 : "The number of atomic orbitals must be the same in every CDFT state.")
1251 : END IF
1252 436 : IF (check_mo(1) /= check_mo(2)) THEN
1253 : CALL cp_abort(__LOCATION__, &
1254 0 : "The number of molecular orbitals must be the same in every CDFT state.")
1255 : END IF
1256 : END DO
1257 : END DO
1258 : ! Allocate work
1259 100 : npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
1260 1426 : ALLOCATE (w_matrix_mo(nforce_eval, nforce_eval, nvar))
1261 782 : ALLOCATE (mo_overlap(npermutations), S_det(npermutations, nspins))
1262 800 : ALLOCATE (a(nspins, nvar, npermutations), b(nspins, nvar, npermutations))
1263 100 : a = 0.0_dp
1264 100 : b = 0.0_dp
1265 100 : IF (mixed_cdft%calculate_metric) THEN
1266 106 : ALLOCATE (density_matrix_diff(npermutations, nspins))
1267 42 : DO ispin = 1, nspins
1268 78 : DO ipermutation = 1, npermutations
1269 36 : NULLIFY (density_matrix_diff(ipermutation, ispin)%matrix)
1270 36 : CALL dbcsr_init_p(density_matrix_diff(ipermutation, ispin)%matrix)
1271 36 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1272 : CALL dbcsr_copy(density_matrix_diff(ipermutation, ispin)%matrix, &
1273 64 : density_matrix(istate, ispin)%matrix, name="DENSITY_MATRIX")
1274 : END DO
1275 : END DO
1276 : END IF
1277 : ! Check for uniform occupations
1278 100 : uniform_occupation = .NOT. ALLOCATED(mixed_cdft%occupations)
1279 100 : should_scale = .FALSE.
1280 100 : IF (.NOT. uniform_occupation) THEN
1281 56 : ALLOCATE (homo(nforce_eval, nspins))
1282 14 : mixed_cdft_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT")
1283 14 : CALL section_vals_val_get(mixed_cdft_section, "EPS_OCCUPIED", r_val=eps_occupied)
1284 14 : IF (eps_occupied > 1.0_dp .OR. eps_occupied < 0.0_dp) THEN
1285 : CALL cp_abort(__LOCATION__, &
1286 0 : "Keyword EPS_OCCUPIED only accepts values between 0.0 and 1.0")
1287 : END IF
1288 14 : IF (mixed_cdft%eps_svd == 0.0_dp) THEN
1289 : CALL cp_warn(__LOCATION__, &
1290 : "The usage of SVD based matrix inversions with fractionally occupied "// &
1291 0 : "orbitals is strongly recommended to screen nearly orthogonal states.")
1292 : END IF
1293 28 : CALL section_vals_val_get(mixed_cdft_section, "SCALE_WITH_OCCUPATION_NUMBERS", l_val=should_scale)
1294 : END IF
1295 : ! Start the actual calculation
1296 294 : DO ispin = 1, nspins
1297 : ! Create the MOxMO fm struct (mo_mo_fm_pools%struct)
1298 : ! The number of MOs/AOs is equal to the number of columns/rows of mo_coeff(:,:)%matrix
1299 194 : NULLIFY (fm_struct_mo, mo_mo_fmstruct)
1300 194 : CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=ncol_mo(ispin), nrow_global=nrow_mo(ispin))
1301 194 : nao = nrow_mo(ispin)
1302 194 : IF (uniform_occupation) THEN
1303 166 : nmo = ncol_mo(ispin)
1304 : ELSE
1305 28 : nmo = ncol_mo(ispin)
1306 : ! Find indices of highest (fractionally) occupied molecular orbital
1307 84 : homo(:, ispin) = nmo
1308 84 : DO istate = 1, nforce_eval
1309 152 : DO j = nmo, 1, -1
1310 124 : IF (mixed_cdft%occupations(istate, ispin)%array(j) >= eps_occupied) THEN
1311 56 : homo(istate, ispin) = j
1312 56 : EXIT
1313 : END IF
1314 : END DO
1315 : END DO
1316 : ! Make matrices square by using the largest homo and emit warning if a state has fewer occupied MOs
1317 : ! Although it would be possible to handle the nonsquare situation as well,
1318 : ! all CDFT states should be in the same spin state for meaningful results
1319 84 : nmo = MAXVAL(homo(:, ispin))
1320 : ! Also check that the number of electrons is conserved (using a fixed sensible threshold)
1321 28 : nelectron_mismatch = .FALSE.
1322 92 : nelectron_tot = SUM(mixed_cdft%occupations(1, ispin)%array(1:nmo))
1323 56 : DO istate = 2, nforce_eval
1324 120 : IF (ABS(SUM(mixed_cdft%occupations(istate, ispin)%array(1:nmo)) - nelectron_tot) > 1.0E-4_dp) THEN
1325 0 : nelectron_mismatch = .TRUE.
1326 : END IF
1327 : END DO
1328 84 : IF (ANY(homo(:, ispin) /= nmo)) THEN
1329 0 : IF (ispin == 1) THEN
1330 : CALL cp_warn(__LOCATION__, &
1331 : "The number of occupied alpha MOs is not constant across all CDFT states. "// &
1332 0 : "Calculation proceeds but the results will likely be meaningless.")
1333 : ELSE
1334 : CALL cp_warn(__LOCATION__, &
1335 : "The number of occupied beta MOs is not constant across all CDFT states. "// &
1336 0 : "Calculation proceeds but the results will likely be meaningless.")
1337 : END IF
1338 28 : ELSE IF (nelectron_mismatch) THEN
1339 0 : IF (ispin == 1) THEN
1340 : CALL cp_warn(__LOCATION__, &
1341 : "The number of alpha electrons is not constant across all CDFT states. "// &
1342 0 : "Calculation proceeds but the results will likely be meaningless.")
1343 : ELSE
1344 : CALL cp_warn(__LOCATION__, &
1345 : "The number of beta electrons is not constant across all CDFT states. "// &
1346 0 : "Calculation proceeds but the results will likely be meaningless.")
1347 : END IF
1348 : END IF
1349 : END IF
1350 : CALL cp_fm_struct_create(fm_struct_mo, nrow_global=nao, ncol_global=nmo, &
1351 194 : context=mixed_cdft%blacs_env, para_env=force_env%para_env)
1352 : CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, &
1353 194 : context=mixed_cdft%blacs_env, para_env=force_env%para_env)
1354 : ! Allocate work
1355 : CALL cp_fm_create(matrix=tmp2, matrix_struct=fm_struct_mo, &
1356 194 : name="ET_TMP_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1357 194 : CALL cp_fm_struct_release(fm_struct_mo)
1358 : CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
1359 194 : name="INVERSE_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1360 : CALL cp_fm_create(matrix=Tinverse, matrix_struct=mo_mo_fmstruct, &
1361 194 : name="T_INVERSE_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1362 552 : DO istate = 1, npermutations
1363 : CALL cp_fm_create(matrix=mo_overlap(istate), matrix_struct=mo_mo_fmstruct, &
1364 : name="MO_OVERLAP_"//TRIM(ADJUSTL(cp_to_string(istate)))//"_"// &
1365 552 : TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1366 : END DO
1367 392 : DO ivar = 1, nvar
1368 836 : DO istate = 1, nforce_eval
1369 1810 : DO jstate = 1, nforce_eval
1370 1168 : IF (istate == jstate) CYCLE
1371 : CALL cp_fm_create(matrix=w_matrix_mo(istate, jstate, ivar), matrix_struct=mo_mo_fmstruct, &
1372 : name="W_"//TRIM(ADJUSTL(cp_to_string(istate)))//"_"// &
1373 : TRIM(ADJUSTL(cp_to_string(jstate)))//"_"// &
1374 1612 : TRIM(ADJUSTL(cp_to_string(ivar)))//"_MATRIX")
1375 : END DO
1376 : END DO
1377 : END DO
1378 194 : CALL cp_fm_struct_release(mo_mo_fmstruct)
1379 : ! Remove empty MOs and (possibly) scale rest with occupation numbers
1380 194 : IF (.NOT. uniform_occupation) THEN
1381 84 : DO iforce_eval = 1, nforce_eval
1382 56 : CALL cp_fm_to_fm(mixed_mo_coeff(iforce_eval, ispin), tmp2, nmo, 1, 1)
1383 56 : CALL cp_fm_release(mixed_mo_coeff(iforce_eval, ispin))
1384 : CALL cp_fm_create(mixed_mo_coeff(iforce_eval, ispin), &
1385 : matrix_struct=tmp2%matrix_struct, &
1386 : name="MO_COEFF_"//TRIM(ADJUSTL(cp_to_string(iforce_eval)))//"_" &
1387 56 : //TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1388 56 : CALL cp_fm_to_fm(tmp2, mixed_mo_coeff(iforce_eval, ispin))
1389 56 : IF (should_scale) THEN
1390 : CALL cp_fm_column_scale(mixed_mo_coeff(iforce_eval, ispin), &
1391 40 : mixed_cdft%occupations(iforce_eval, ispin)%array(1:nmo))
1392 : END IF
1393 84 : DEALLOCATE (mixed_cdft%occupations(iforce_eval, ispin)%array)
1394 : END DO
1395 : END IF
1396 : ! calculate the MO overlaps (C_j)^T S C_i
1397 194 : ipermutation = 0
1398 630 : DO istate = 1, nforce_eval
1399 988 : DO jstate = istate + 1, nforce_eval
1400 358 : ipermutation = ipermutation + 1
1401 : CALL cp_dbcsr_sm_fm_multiply(mixed_matrix_s, mixed_mo_coeff(istate, ispin), &
1402 358 : tmp2, nmo, 1.0_dp, 0.0_dp)
1403 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
1404 : mixed_mo_coeff(jstate, ispin), &
1405 358 : tmp2, 0.0_dp, mo_overlap(ipermutation))
1406 794 : IF (print_mo) THEN
1407 : CALL cp_fm_write_formatted(mo_overlap(ipermutation), mounit, &
1408 : "# MO overlap matrix (step "//TRIM(ADJUSTL(cp_to_string(mixed_cdft%sim_step)))// &
1409 : "): CDFT states "//TRIM(ADJUSTL(cp_to_string(istate)))//" and "// &
1410 : TRIM(ADJUSTL(cp_to_string(jstate)))//" (spin "// &
1411 4 : TRIM(ADJUSTL(cp_to_string(ispin)))//")")
1412 : END IF
1413 : END DO
1414 : END DO
1415 : ! calculate the MO-representations of the restraint matrices of all CDFT states
1416 392 : DO ivar = 1, nvar
1417 836 : DO jstate = 1, nforce_eval
1418 1810 : DO istate = 1, nforce_eval
1419 1168 : IF (istate == jstate) CYCLE
1420 : ! State i: (C_j)^T W_i C_i
1421 : CALL cp_dbcsr_sm_fm_multiply(w_matrix(istate, ivar)%matrix, &
1422 : mixed_mo_coeff(istate, ispin), &
1423 724 : tmp2, nmo, 1.0_dp, 0.0_dp)
1424 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
1425 : mixed_mo_coeff(jstate, ispin), &
1426 1612 : tmp2, 0.0_dp, w_matrix_mo(istate, jstate, ivar))
1427 : END DO
1428 : END DO
1429 : END DO
1430 552 : DO ipermutation = 1, npermutations
1431 : ! Invert and calculate determinant of MO overlaps
1432 358 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1433 358 : IF (print_mo_eigval) THEN
1434 28 : NULLIFY (mo_eigval)
1435 : CALL cp_fm_invert(mo_overlap(ipermutation), inverse_mat, &
1436 : S_det(ipermutation, ispin), eps_svd=mixed_cdft%eps_svd, &
1437 28 : eigval=mo_eigval)
1438 28 : IF (moeigvalunit > 0) THEN
1439 14 : IF (mixed_cdft%eps_svd == 0.0_dp) THEN
1440 : WRITE (moeigvalunit, '(A,I2,A,I2,A,I1,A)') &
1441 0 : "# MO Overlap matrix eigenvalues for CDFT states ", istate, " and ", jstate, &
1442 0 : " (spin ", ispin, ")"
1443 : ELSE
1444 : WRITE (moeigvalunit, '(A,I2,A,I2,A,I1,A)') &
1445 14 : "# MO Overlap matrix singular values for CDFT states ", istate, " and ", jstate, &
1446 28 : " (spin ", ispin, ")"
1447 : END IF
1448 14 : WRITE (moeigvalunit, '(A1, A9, A12)') "#", "Index", ADJUSTL("Value")
1449 46 : DO j = 1, SIZE(mo_eigval)
1450 46 : WRITE (moeigvalunit, '(I10, F12.8)') j, mo_eigval(j)
1451 : END DO
1452 : END IF
1453 28 : DEALLOCATE (mo_eigval)
1454 : ELSE
1455 : CALL cp_fm_invert(mo_overlap(ipermutation), inverse_mat, &
1456 330 : S_det(ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
1457 : END IF
1458 358 : CALL cp_fm_get_info(inverse_mat, nrow_local=nrow_local, ncol_local=ncol_local)
1459 : ! Calculate <Psi_i | w_j(r) | Psi_j> for ivar:th constraint
1460 924 : DO j = 1, ncol_local
1461 1437 : DO k = 1, nrow_local
1462 1597 : DO ivar = 1, nvar
1463 : b(ispin, ivar, ipermutation) = b(ispin, ivar, ipermutation) + &
1464 : w_matrix_mo(jstate, istate, ivar)%local_data(k, j)* &
1465 1031 : inverse_mat%local_data(k, j)
1466 : END DO
1467 : END DO
1468 : END DO
1469 : ! Calculate <Psi_j | w_i(r) | Psi_i> for ivar:th constraint
1470 358 : CALL cp_fm_transpose(inverse_mat, Tinverse)
1471 924 : DO j = 1, ncol_local
1472 1437 : DO k = 1, nrow_local
1473 1597 : DO ivar = 1, nvar
1474 : a(ispin, ivar, ipermutation) = a(ispin, ivar, ipermutation) + &
1475 : w_matrix_mo(istate, jstate, ivar)%local_data(k, j)* &
1476 1031 : Tinverse%local_data(k, j)
1477 : END DO
1478 : END DO
1479 : END DO
1480 : ! Handle different constraint types
1481 720 : DO ivar = 1, nvar
1482 362 : SELECT CASE (mixed_cdft%constraint_type(ivar, istate))
1483 : CASE (cdft_charge_constraint)
1484 : ! No action needed
1485 : CASE (cdft_magnetization_constraint)
1486 0 : IF (ispin == 2) a(ispin, ivar, ipermutation) = -a(ispin, ivar, ipermutation)
1487 : CASE (cdft_alpha_constraint)
1488 : ! Constraint applied to alpha electrons only, set integrals involving beta to zero
1489 4 : IF (ispin == 2) a(ispin, ivar, ipermutation) = 0.0_dp
1490 : CASE (cdft_beta_constraint)
1491 : ! Constraint applied to beta electrons only, set integrals involving alpha to zero
1492 4 : IF (ispin == 1) a(ispin, ivar, ipermutation) = 0.0_dp
1493 : CASE DEFAULT
1494 362 : CPABORT("Unknown constraint type.")
1495 : END SELECT
1496 358 : SELECT CASE (mixed_cdft%constraint_type(ivar, jstate))
1497 : CASE (cdft_charge_constraint)
1498 : ! No action needed
1499 : CASE (cdft_magnetization_constraint)
1500 0 : IF (ispin == 2) b(ispin, ivar, ipermutation) = -b(ispin, ivar, ipermutation)
1501 : CASE (cdft_alpha_constraint)
1502 : ! Constraint applied to alpha electrons only, set integrals involving beta to zero
1503 4 : IF (ispin == 2) b(ispin, ivar, ipermutation) = 0.0_dp
1504 : CASE (cdft_beta_constraint)
1505 : ! Constraint applied to beta electrons only, set integrals involving alpha to zero
1506 4 : IF (ispin == 1) b(ispin, ivar, ipermutation) = 0.0_dp
1507 : CASE DEFAULT
1508 362 : CPABORT("Unknown constraint type.")
1509 : END SELECT
1510 : END DO
1511 : ! Compute density matrix difference P = P_j - P_i
1512 358 : IF (mixed_cdft%calculate_metric) THEN
1513 : CALL dbcsr_add(density_matrix_diff(ipermutation, ispin)%matrix, &
1514 36 : density_matrix(jstate, ispin)%matrix, -1.0_dp, 1.0_dp)
1515 : END IF
1516 : !
1517 1082 : CALL force_env%para_env%sum(a(ispin, :, ipermutation))
1518 1992 : CALL force_env%para_env%sum(b(ispin, :, ipermutation))
1519 : END DO
1520 : ! Release work
1521 194 : CALL cp_fm_release(tmp2)
1522 392 : DO ivar = 1, nvar
1523 836 : DO jstate = 1, nforce_eval
1524 1810 : DO istate = 1, nforce_eval
1525 1168 : IF (istate == jstate) CYCLE
1526 1612 : CALL cp_fm_release(w_matrix_mo(istate, jstate, ivar))
1527 : END DO
1528 : END DO
1529 : END DO
1530 552 : DO ipermutation = 1, npermutations
1531 552 : CALL cp_fm_release(mo_overlap(ipermutation))
1532 : END DO
1533 194 : CALL cp_fm_release(Tinverse)
1534 294 : CALL cp_fm_release(inverse_mat)
1535 : END DO
1536 100 : DEALLOCATE (mo_overlap)
1537 100 : DEALLOCATE (w_matrix_mo)
1538 100 : IF (.NOT. uniform_occupation) THEN
1539 14 : DEALLOCATE (homo)
1540 14 : DEALLOCATE (mixed_cdft%occupations)
1541 : END IF
1542 100 : IF (print_mo) THEN
1543 : CALL cp_print_key_finished_output(mounit, logger, force_env_section, &
1544 2 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO", on_file=.TRUE.)
1545 : END IF
1546 100 : IF (print_mo_eigval) THEN
1547 : CALL cp_print_key_finished_output(moeigvalunit, logger, force_env_section, &
1548 14 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO", on_file=.TRUE.)
1549 : END IF
1550 : ! solve eigenstates for the projector matrix
1551 400 : ALLOCATE (Wda(nvar, npermutations))
1552 300 : ALLOCATE (Sda(npermutations))
1553 104 : IF (.NOT. mixed_cdft%identical_constraints) ALLOCATE (Wad(nvar, npermutations))
1554 282 : DO ipermutation = 1, npermutations
1555 182 : IF (nspins == 2) THEN
1556 176 : Sda(ipermutation) = ABS(S_det(ipermutation, 1)*S_det(ipermutation, 2))
1557 : ELSE
1558 6 : Sda(ipermutation) = S_det(ipermutation, 1)**2
1559 : END IF
1560 : ! Finalize <Psi_j | w_i(r) | Psi_i> by multiplication with Sda
1561 466 : DO ivar = 1, nvar
1562 366 : IF (mixed_cdft%identical_constraints) THEN
1563 : Wda(ivar, ipermutation) = (SUM(a(:, ivar, ipermutation)) + SUM(b(:, ivar, ipermutation)))* &
1564 898 : Sda(ipermutation)/2.0_dp
1565 : ELSE
1566 6 : Wda(ivar, ipermutation) = SUM(a(:, ivar, ipermutation))*Sda(ipermutation)
1567 6 : Wad(ivar, ipermutation) = SUM(b(:, ivar, ipermutation))*Sda(ipermutation)
1568 : END IF
1569 : END DO
1570 : END DO
1571 100 : DEALLOCATE (a, b, S_det)
1572 : ! Transfer info about the constraint calculations
1573 800 : ALLOCATE (W_diagonal(nvar, nforce_eval), strength(nvar, nforce_eval), energy(nforce_eval))
1574 100 : W_diagonal = 0.0_dp
1575 324 : DO iforce_eval = 1, nforce_eval
1576 552 : strength(:, iforce_eval) = mixed_env%strength(iforce_eval, :)
1577 : END DO
1578 100 : energy = 0.0_dp
1579 324 : DO iforce_eval = 1, nforce_eval
1580 224 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
1581 186 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1582 24 : qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1583 : ELSE
1584 162 : CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1585 : END IF
1586 186 : CALL get_qs_env(qs_env, energy=energy_qs, dft_control=dft_control)
1587 286 : IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1588 226 : W_diagonal(:, iforce_eval) = dft_control%qs_control%cdft_control%value(:)
1589 150 : energy(iforce_eval) = energy_qs%total
1590 : END IF
1591 : END DO
1592 100 : CALL force_env%para_env%sum(W_diagonal)
1593 100 : CALL force_env%para_env%sum(energy)
1594 : CALL mixed_cdft_result_type_set(mixed_cdft%results, Wda=Wda, W_diagonal=W_diagonal, &
1595 100 : energy=energy, strength=strength)
1596 100 : IF (.NOT. mixed_cdft%identical_constraints) CALL mixed_cdft_result_type_set(mixed_cdft%results, Wad=Wad)
1597 : ! Construct S
1598 400 : ALLOCATE (S_mat(nforce_eval, nforce_eval))
1599 324 : DO istate = 1, nforce_eval
1600 324 : S_mat(istate, istate) = 1.0_dp
1601 : END DO
1602 282 : DO ipermutation = 1, npermutations
1603 182 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1604 182 : S_mat(istate, jstate) = Sda(ipermutation)
1605 282 : S_mat(jstate, istate) = Sda(ipermutation)
1606 : END DO
1607 100 : CALL mixed_cdft_result_type_set(mixed_cdft%results, S=S_mat)
1608 : ! Invert S via eigendecomposition and compute S^-(1/2)
1609 400 : ALLOCATE (eigenv(nforce_eval), tmp_mat(nforce_eval, nforce_eval))
1610 100 : CALL diamat_all(S_mat, eigenv, .TRUE.)
1611 100 : tmp_mat = 0.0_dp
1612 324 : DO istate = 1, nforce_eval
1613 224 : IF (eigenv(istate) < 1.0e-14_dp) THEN
1614 : ! Safeguard against division with 0 and negative numbers
1615 10 : eigenv(istate) = 1.0e-14_dp
1616 : CALL cp_warn(__LOCATION__, &
1617 : "The overlap matrix is numerically nearly singular. "// &
1618 10 : "Calculation proceeds but the results might be meaningless.")
1619 : END IF
1620 324 : tmp_mat(istate, istate) = 1.0_dp/SQRT(eigenv(istate))
1621 : END DO
1622 5940 : tmp_mat(:, :) = MATMUL(tmp_mat, TRANSPOSE(S_mat))
1623 10356 : S_mat(:, :) = MATMUL(S_mat, tmp_mat) ! S^(-1/2)
1624 100 : CALL mixed_cdft_result_type_set(mixed_cdft%results, S_minushalf=S_mat)
1625 100 : DEALLOCATE (eigenv, tmp_mat, S_mat)
1626 : ! Construct nonorthogonal diabatic Hamiltonian matrix H''
1627 300 : ALLOCATE (H_mat(nforce_eval, nforce_eval))
1628 118 : IF (mixed_cdft%nonortho_coupling) ALLOCATE (coupling_nonortho(npermutations))
1629 324 : DO istate = 1, nforce_eval
1630 324 : H_mat(istate, istate) = energy(istate)
1631 : END DO
1632 282 : DO ipermutation = 1, npermutations
1633 182 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1634 182 : sum_a = 0.0_dp
1635 182 : sum_b = 0.0_dp
1636 366 : DO ivar = 1, nvar
1637 : ! V_J * <Psi_J | w_J(r) | Psi_J>
1638 184 : sum_b(1) = sum_b(1) + strength(ivar, jstate)*W_diagonal(ivar, jstate)
1639 : ! V_I * <Psi_I | w_I(r) | Psi_I>
1640 184 : sum_a(1) = sum_a(1) + strength(ivar, istate)*W_diagonal(ivar, istate)
1641 366 : IF (mixed_cdft%identical_constraints) THEN
1642 : ! V_J * W_IJ
1643 182 : sum_b(2) = sum_b(2) + strength(ivar, jstate)*Wda(ivar, ipermutation)
1644 : ! V_I * W_JI
1645 182 : sum_a(2) = sum_a(2) + strength(ivar, istate)*Wda(ivar, ipermutation)
1646 : ELSE
1647 : ! V_J * W_IJ
1648 2 : sum_b(2) = sum_b(2) + strength(ivar, jstate)*Wad(ivar, ipermutation)
1649 : ! V_I * W_JI
1650 2 : sum_a(2) = sum_a(2) + strength(ivar, istate)*Wda(ivar, ipermutation)
1651 : END IF
1652 : END DO
1653 : ! Denote F_X = <Psi_X | H_X + V_X*w_X(r) | Psi_X> = E_X + V_X*<Psi_X | w_X(r) | Psi_X>
1654 : ! H_IJ = F_J*S_IJ - V_J * W_IJ
1655 182 : c(1) = (energy(jstate) + sum_b(1))*Sda(ipermutation) - sum_b(2)
1656 : ! H_JI = F_I*S_JI - V_I * W_JI
1657 182 : c(2) = (energy(istate) + sum_a(1))*Sda(ipermutation) - sum_a(2)
1658 : ! H''(I,J) = 0.5*(H_IJ+H_JI) = H''(J,I)
1659 182 : H_mat(istate, jstate) = (c(1) + c(2))*0.5_dp
1660 182 : H_mat(jstate, istate) = H_mat(istate, jstate)
1661 464 : IF (mixed_cdft%nonortho_coupling) coupling_nonortho(ipermutation) = H_mat(istate, jstate)
1662 : END DO
1663 100 : CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat)
1664 100 : DEALLOCATE (H_mat, W_diagonal, Wda, strength, energy, Sda)
1665 100 : IF (ALLOCATED(Wad)) DEALLOCATE (Wad)
1666 100 : IF (mixed_cdft%nonortho_coupling) THEN
1667 18 : CALL mixed_cdft_result_type_set(mixed_cdft%results, nonortho=coupling_nonortho)
1668 18 : DEALLOCATE (coupling_nonortho)
1669 : END IF
1670 : ! Compute metric to assess reliability of coupling
1671 100 : IF (mixed_cdft%calculate_metric) CALL mixed_cdft_calculate_metric(force_env, mixed_cdft, density_matrix_diff, ncol_mo)
1672 : ! Compute coupling also with the wavefunction overlap method, see Migliore2009
1673 : ! Requires the unconstrained KS ground state wavefunction as input
1674 100 : IF (mixed_cdft%wfn_overlap_method) THEN
1675 8 : IF (.NOT. uniform_occupation) THEN
1676 : CALL cp_abort(__LOCATION__, &
1677 0 : "Wavefunction overlap method supports only uniformly occupied MOs.")
1678 : END IF
1679 8 : CALL mixed_cdft_wfn_overlap_method(force_env, mixed_cdft, ncol_mo, nrow_mo)
1680 : END IF
1681 : ! Release remaining work
1682 100 : DEALLOCATE (nrow_mo, ncol_mo)
1683 100 : CALL mixed_cdft_work_type_release(mixed_cdft%matrix)
1684 100 : CALL timestop(handle)
1685 :
1686 300 : END SUBROUTINE mixed_cdft_interaction_matrices
1687 :
1688 : ! **************************************************************************************************
1689 : !> \brief Routine to calculate the CDFT electronic couplings.
1690 : !> \param force_env the force_env that holds the CDFT states
1691 : !> \par History
1692 : !> 11.17 created [Nico Holmberg]
1693 : ! **************************************************************************************************
1694 100 : SUBROUTINE mixed_cdft_calculate_coupling_low(force_env)
1695 : TYPE(force_env_type), POINTER :: force_env
1696 :
1697 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_calculate_coupling_low'
1698 :
1699 : INTEGER :: handle, ipermutation, istate, jstate, &
1700 : nforce_eval, npermutations, nvar
1701 : LOGICAL :: use_lowdin, use_rotation
1702 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_lowdin, coupling_rotation, &
1703 100 : eigenv
1704 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, W_mat
1705 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1706 :
1707 100 : NULLIFY (mixed_cdft)
1708 100 : CPASSERT(ASSOCIATED(force_env))
1709 100 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1710 100 : CALL timeset(routineN, handle)
1711 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1712 100 : CPASSERT(ALLOCATED(mixed_cdft%results%W_diagonal))
1713 100 : CPASSERT(ALLOCATED(mixed_cdft%results%Wda))
1714 100 : CPASSERT(ALLOCATED(mixed_cdft%results%S_minushalf))
1715 100 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1716 : ! Decide which methods to use for computing the coupling
1717 : ! Default behavior is to use rotation when a single constraint is active, otherwise uses Lowdin orthogonalization
1718 : ! The latter can also be explicitly requested when a single constraint is active
1719 : ! Possibly computes the coupling additionally with the wavefunction overlap method
1720 100 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
1721 100 : nvar = SIZE(mixed_cdft%results%Wda, 1)
1722 100 : npermutations = nforce_eval*(nforce_eval - 1)/2
1723 400 : ALLOCATE (tmp_mat(nforce_eval, nforce_eval))
1724 100 : IF (nvar == 1 .AND. mixed_cdft%identical_constraints) THEN
1725 96 : use_rotation = .TRUE.
1726 96 : use_lowdin = mixed_cdft%use_lowdin
1727 : ELSE
1728 : use_rotation = .FALSE.
1729 : use_lowdin = .TRUE.
1730 : END IF
1731 : ! Calculate coupling by rotating the CDFT states to eigenstates of the weight matrix W (single constraint only)
1732 : IF (use_rotation) THEN
1733 : ! Construct W
1734 480 : ALLOCATE (W_mat(nforce_eval, nforce_eval), coupling_rotation(npermutations))
1735 288 : ALLOCATE (eigenv(nforce_eval))
1736 : ! W_mat(i, i) = N_i where N_i is the value of the constraint in state i
1737 312 : DO istate = 1, nforce_eval
1738 528 : W_mat(istate, istate) = SUM(mixed_cdft%results%W_diagonal(:, istate))
1739 : END DO
1740 : ! W_mat(i, j) = <Psi_i|w(r)|Psi_j>
1741 274 : DO ipermutation = 1, npermutations
1742 178 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1743 356 : W_mat(istate, jstate) = SUM(mixed_cdft%results%Wda(:, ipermutation))
1744 452 : W_mat(jstate, istate) = W_mat(istate, jstate)
1745 : END DO
1746 : ! Solve generalized eigenvalue equation WV = SVL
1747 : ! Convert to standard eigenvalue problem via symmetric orthogonalisation
1748 5036 : tmp_mat(:, :) = MATMUL(W_mat, mixed_cdft%results%S_minushalf) ! W * S^(-1/2)
1749 5036 : W_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, tmp_mat) ! W' = S^(-1/2) * W * S^(-1/2)
1750 96 : CALL diamat_all(W_mat, eigenv, .TRUE.) ! Solve W'V' = AV'
1751 9188 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, W_mat) ! Reverse transformation V = S^(-1/2) V'
1752 : ! Construct final, orthogonal diabatic Hamiltonian matrix H
1753 9188 : W_mat(:, :) = MATMUL(mixed_cdft%results%H, tmp_mat) ! H'' * V
1754 10168 : W_mat(:, :) = MATMUL(TRANSPOSE(tmp_mat), W_mat) ! H = V^T * H'' * V
1755 274 : DO ipermutation = 1, npermutations
1756 178 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1757 274 : coupling_rotation(ipermutation) = W_mat(istate, jstate)
1758 : END DO
1759 96 : CALL mixed_cdft_result_type_set(mixed_cdft%results, rotation=coupling_rotation)
1760 96 : DEALLOCATE (W_mat, coupling_rotation, eigenv)
1761 : END IF
1762 : ! Calculate coupling by Lowdin orthogonalization
1763 100 : IF (use_lowdin) THEN
1764 60 : ALLOCATE (coupling_lowdin(npermutations))
1765 780 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%H, mixed_cdft%results%S_minushalf) ! H'' * S^(-1/2)
1766 : ! Final orthogonal diabatic Hamiltonian matrix H
1767 580 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, tmp_mat) ! H = S^(-1/2) * H'' * S^(-1/2)
1768 40 : DO ipermutation = 1, npermutations
1769 20 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1770 40 : coupling_lowdin(ipermutation) = tmp_mat(istate, jstate)
1771 : END DO
1772 20 : CALL mixed_cdft_result_type_set(mixed_cdft%results, lowdin=coupling_lowdin)
1773 20 : DEALLOCATE (coupling_lowdin)
1774 : END IF
1775 100 : DEALLOCATE (tmp_mat)
1776 100 : CALL timestop(handle)
1777 :
1778 200 : END SUBROUTINE mixed_cdft_calculate_coupling_low
1779 :
1780 : ! **************************************************************************************************
1781 : !> \brief Performs a configuration interaction calculation in the basis spanned by the CDFT states.
1782 : !> \param force_env the force_env that holds the CDFT states
1783 : !> \par History
1784 : !> 11.17 created [Nico Holmberg]
1785 : ! **************************************************************************************************
1786 100 : SUBROUTINE mixed_cdft_configuration_interaction(force_env)
1787 : TYPE(force_env_type), POINTER :: force_env
1788 :
1789 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_configuration_interaction'
1790 :
1791 : INTEGER :: handle, info, iounit, istate, ivar, &
1792 : nforce_eval, work_array_size
1793 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenv, work
1794 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, H_mat_copy, S_mat, S_mat_copy
1795 : REAL(KIND=dp), EXTERNAL :: dnrm2
1796 : TYPE(cp_logger_type), POINTER :: logger
1797 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1798 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
1799 :
1800 : EXTERNAL :: dsygv
1801 :
1802 100 : NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1803 :
1804 100 : CPASSERT(ASSOCIATED(force_env))
1805 100 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1806 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1807 :
1808 100 : IF (.NOT. mixed_cdft%do_ci) RETURN
1809 :
1810 20 : logger => cp_get_default_logger()
1811 20 : CALL timeset(routineN, handle)
1812 : CALL force_env_get(force_env=force_env, &
1813 20 : force_env_section=force_env_section)
1814 20 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1815 20 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1816 :
1817 20 : CPASSERT(ALLOCATED(mixed_cdft%results%S))
1818 20 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1819 20 : nforce_eval = SIZE(mixed_cdft%results%S, 1)
1820 120 : ALLOCATE (S_mat(nforce_eval, nforce_eval), H_mat(nforce_eval, nforce_eval))
1821 60 : ALLOCATE (eigenv(nforce_eval))
1822 168 : S_mat(:, :) = mixed_cdft%results%S(:, :)
1823 168 : H_mat(:, :) = mixed_cdft%results%H(:, :)
1824 : ! Workspace query
1825 20 : ALLOCATE (work(1))
1826 20 : info = 0
1827 100 : ALLOCATE (H_mat_copy(nforce_eval, nforce_eval), S_mat_copy(nforce_eval, nforce_eval))
1828 168 : H_mat_copy(:, :) = H_mat(:, :) ! Need explicit copies because dsygv destroys original values
1829 168 : S_mat_copy(:, :) = S_mat(:, :)
1830 20 : CALL dsygv(1, 'V', 'U', nforce_eval, H_mat_copy, nforce_eval, S_mat_copy, nforce_eval, eigenv, work, -1, info)
1831 20 : work_array_size = NINT(work(1))
1832 20 : DEALLOCATE (H_mat_copy, S_mat_copy)
1833 : ! Allocate work array
1834 20 : DEALLOCATE (work)
1835 60 : ALLOCATE (work(work_array_size))
1836 20 : work = 0.0_dp
1837 : ! Solve Hc = eSc
1838 20 : info = 0
1839 20 : CALL dsygv(1, 'V', 'U', nforce_eval, H_mat, nforce_eval, S_mat, nforce_eval, eigenv, work, work_array_size, info)
1840 20 : IF (info /= 0) THEN
1841 0 : IF (info > nforce_eval) THEN
1842 0 : CPABORT("Matrix S is not positive definite")
1843 : ELSE
1844 0 : CPABORT("Diagonalization of H matrix failed.")
1845 : END IF
1846 : END IF
1847 : ! dsygv returns eigenvectors (stored in columns of H_mat) that are normalized to H^T * S * H = I
1848 : ! Renormalize eigenvectors to H^T * H = I
1849 64 : DO ivar = 1, nforce_eval
1850 168 : H_mat(:, ivar) = H_mat(:, ivar)/dnrm2(nforce_eval, H_mat(:, ivar), 1)
1851 : END DO
1852 20 : DEALLOCATE (work)
1853 20 : IF (iounit > 0) THEN
1854 10 : WRITE (iounit, '(/,T3,A)') '------------------ CDFT Configuration Interaction (CDFT-CI) ------------------'
1855 32 : DO ivar = 1, nforce_eval
1856 22 : IF (ivar == 1) THEN
1857 10 : WRITE (iounit, '(T3,A,T58,(3X,F20.14))') 'Ground state energy:', eigenv(ivar)
1858 : ELSE
1859 12 : WRITE (iounit, '(/,T3,A,I2,A,T58,(3X,F20.14))') 'Excited state (', ivar - 1, ' ) energy:', eigenv(ivar)
1860 : END IF
1861 58 : DO istate = 1, nforce_eval, 2
1862 48 : IF (istate == 1) THEN
1863 : WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
1864 22 : 'Expansion coefficients:', H_mat(istate, ivar), H_mat(istate + 1, ivar)
1865 4 : ELSE IF (istate < nforce_eval) THEN
1866 4 : WRITE (iounit, '(T54,(3X,2F12.6))') H_mat(istate, ivar), H_mat(istate + 1, ivar)
1867 : ELSE
1868 0 : WRITE (iounit, '(T54,(3X,F12.6))') H_mat(istate, ivar)
1869 : END IF
1870 : END DO
1871 : END DO
1872 : WRITE (iounit, '(T3,A)') &
1873 10 : '------------------------------------------------------------------------------'
1874 : END IF
1875 20 : DEALLOCATE (S_mat, H_mat, eigenv)
1876 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1877 20 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1878 20 : CALL timestop(handle)
1879 :
1880 20 : END SUBROUTINE mixed_cdft_configuration_interaction
1881 : ! **************************************************************************************************
1882 : !> \brief Block diagonalizes the mixed CDFT Hamiltonian matrix.
1883 : !> \param force_env the force_env that holds the CDFT states
1884 : !> \par History
1885 : !> 11.17 created [Nico Holmberg]
1886 : !> 01.18 added recursive diagonalization
1887 : !> split to subroutines [Nico Holmberg]
1888 : ! **************************************************************************************************
1889 100 : SUBROUTINE mixed_cdft_block_diag(force_env)
1890 : TYPE(force_env_type), POINTER :: force_env
1891 :
1892 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_block_diag'
1893 :
1894 : INTEGER :: handle, i, iounit, irecursion, j, n, &
1895 : nblk, nforce_eval, nrecursion
1896 : LOGICAL :: ignore_excited
1897 100 : TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1898 100 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1899 100 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: H_block, S_block
1900 : TYPE(cp_logger_type), POINTER :: logger
1901 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1902 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
1903 :
1904 : EXTERNAL :: dsygv
1905 :
1906 100 : NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1907 :
1908 100 : CPASSERT(ASSOCIATED(force_env))
1909 100 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1910 100 : CPASSERT(ASSOCIATED(mixed_cdft))
1911 :
1912 100 : IF (.NOT. mixed_cdft%block_diagonalize) RETURN
1913 :
1914 8 : logger => cp_get_default_logger()
1915 8 : CALL timeset(routineN, handle)
1916 :
1917 8 : CPASSERT(ALLOCATED(mixed_cdft%results%S))
1918 8 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1919 8 : nforce_eval = SIZE(mixed_cdft%results%S, 1)
1920 :
1921 : CALL force_env_get(force_env=force_env, &
1922 8 : force_env_section=force_env_section)
1923 8 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1924 8 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1925 : ! Read block definitions from input
1926 8 : CALL mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
1927 8 : nblk = SIZE(blocks)
1928 : ! Start block diagonalization
1929 18 : DO irecursion = 1, nrecursion
1930 : ! Print block definitions
1931 10 : IF (iounit > 0 .AND. irecursion == 1) THEN
1932 4 : WRITE (iounit, '(/,T3,A)') '-------------------------- CDFT BLOCK DIAGONALIZATION ------------------------'
1933 4 : WRITE (iounit, '(T3,A)') 'Block diagonalizing the mixed CDFT Hamiltonian'
1934 4 : WRITE (iounit, '(T3,A,I3)') 'Number of blocks:', nblk
1935 4 : WRITE (iounit, '(T3,A,L3)') 'Ignoring excited states within blocks:', ignore_excited
1936 4 : WRITE (iounit, '(/,T3,A)') 'List of CDFT states for each block'
1937 14 : DO i = 1, nblk
1938 33 : WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array))
1939 : END DO
1940 : END IF
1941 : ! Recursive diagonalization: update counters and references
1942 10 : IF (irecursion > 1) THEN
1943 2 : nblk = nblk/2
1944 10 : ALLOCATE (blocks(nblk))
1945 2 : j = 1
1946 6 : DO i = 1, nblk
1947 4 : NULLIFY (blocks(i)%array)
1948 4 : ALLOCATE (blocks(i)%array(2))
1949 12 : blocks(i)%array = [j, j + 1]
1950 6 : j = j + 2
1951 : END DO
1952 : ! Print info
1953 2 : IF (iounit > 0) THEN
1954 1 : WRITE (iounit, '(/, T3,A)') 'Recursive block diagonalization of the mixed CDFT Hamiltonian'
1955 1 : WRITE (iounit, '(T6,A)') 'Block diagonalization is continued until only two matrix blocks remain.'
1956 1 : WRITE (iounit, '(T6,A)') 'The new blocks are formed by collecting pairs of blocks from the previous'
1957 1 : WRITE (iounit, '(T6,A)') 'block diagonalized matrix in ascending order.'
1958 1 : WRITE (iounit, '(/,T3,A,I3,A,I3)') 'Recursion step:', irecursion - 1, ' of ', nrecursion - 1
1959 1 : WRITE (iounit, '(/,T3,A)') 'List of old block indices for each new block'
1960 3 : DO i = 1, nblk
1961 7 : WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array))
1962 : END DO
1963 : END IF
1964 : END IF
1965 : ! Get the Hamiltonian and overlap matrices of each block
1966 10 : CALL mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block)
1967 : ! Diagonalize blocks
1968 10 : CALL mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues)
1969 : ! Assemble the block diagonalized matrices
1970 10 : IF (ignore_excited) THEN
1971 8 : n = nblk
1972 : ELSE
1973 2 : n = nforce_eval
1974 : END IF
1975 10 : CALL mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, n, iounit)
1976 : ! Deallocate work
1977 34 : DO i = 1, nblk
1978 24 : DEALLOCATE (H_block(i)%array)
1979 24 : DEALLOCATE (S_block(i)%array)
1980 24 : DEALLOCATE (eigenvalues(i)%array)
1981 34 : DEALLOCATE (blocks(i)%array)
1982 : END DO
1983 18 : DEALLOCATE (H_block, S_block, eigenvalues, blocks)
1984 : END DO ! recursion
1985 8 : IF (iounit > 0) THEN
1986 : WRITE (iounit, '(T3,A)') &
1987 4 : '------------------------------------------------------------------------------'
1988 : END IF
1989 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1990 8 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1991 8 : CALL timestop(handle)
1992 :
1993 100 : END SUBROUTINE mixed_cdft_block_diag
1994 : ! **************************************************************************************************
1995 : !> \brief Routine to calculate the CDFT electronic coupling reliability metric
1996 : !> \param force_env the force_env that holds the CDFT states
1997 : !> \param mixed_cdft the mixed_cdft env
1998 : !> \param density_matrix_diff array holding difference density matrices (P_j - P_i) for every CDFT
1999 : !> state permutation
2000 : !> \param ncol_mo the number of MOs per spin
2001 : !> \par History
2002 : !> 11.17 created [Nico Holmberg]
2003 : ! **************************************************************************************************
2004 14 : SUBROUTINE mixed_cdft_calculate_metric(force_env, mixed_cdft, density_matrix_diff, ncol_mo)
2005 : TYPE(force_env_type), POINTER :: force_env
2006 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2007 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix_diff
2008 : INTEGER, DIMENSION(:) :: ncol_mo
2009 :
2010 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_calculate_metric'
2011 :
2012 : INTEGER :: handle, ipermutation, ispin, j, &
2013 : nforce_eval, npermutations, nspins
2014 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
2015 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: metric
2016 : TYPE(dbcsr_type) :: e_vectors
2017 :
2018 14 : CALL timeset(routineN, handle)
2019 14 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
2020 14 : npermutations = nforce_eval*(nforce_eval - 1)/2
2021 14 : nspins = SIZE(density_matrix_diff, 2)
2022 56 : ALLOCATE (metric(npermutations, nspins))
2023 14 : metric = 0.0_dp
2024 14 : CALL dbcsr_create(e_vectors, template=density_matrix_diff(1, 1)%matrix)
2025 42 : DO ispin = 1, nspins
2026 84 : ALLOCATE (evals(ncol_mo(ispin)))
2027 64 : DO ipermutation = 1, npermutations
2028 : ! Take into account doubly occupied orbitals without LSD
2029 36 : IF (nspins == 1) THEN
2030 0 : CALL dbcsr_scale(density_matrix_diff(ipermutation, 1)%matrix, alpha_scalar=0.5_dp)
2031 : END IF
2032 : ! Diagonalize difference density matrix
2033 : CALL cp_dbcsr_syevd(density_matrix_diff(ipermutation, ispin)%matrix, e_vectors, evals, &
2034 36 : para_env=force_env%para_env, blacs_env=mixed_cdft%blacs_env)
2035 36 : CALL dbcsr_release_p(density_matrix_diff(ipermutation, ispin)%matrix)
2036 100 : DO j = 1, ncol_mo(ispin)
2037 72 : metric(ipermutation, ispin) = metric(ipermutation, ispin) + (evals(j)**2 - evals(j)**4)
2038 : END DO
2039 : END DO
2040 42 : DEALLOCATE (evals)
2041 : END DO
2042 14 : CALL dbcsr_release(e_vectors)
2043 14 : DEALLOCATE (density_matrix_diff)
2044 78 : metric(:, :) = metric(:, :)/4.0_dp
2045 14 : CALL mixed_cdft_result_type_set(mixed_cdft%results, metric=metric)
2046 14 : DEALLOCATE (metric)
2047 14 : CALL timestop(handle)
2048 :
2049 28 : END SUBROUTINE mixed_cdft_calculate_metric
2050 :
2051 : ! **************************************************************************************************
2052 : !> \brief Routine to calculate the electronic coupling according to the wavefunction overlap method
2053 : !> \param force_env the force_env that holds the CDFT states
2054 : !> \param mixed_cdft the mixed_cdft env
2055 : !> \param ncol_mo the number of MOs per spin
2056 : !> \param nrow_mo the number of AOs per spin
2057 : !> \par History
2058 : !> 11.17 created [Nico Holmberg]
2059 : ! **************************************************************************************************
2060 8 : SUBROUTINE mixed_cdft_wfn_overlap_method(force_env, mixed_cdft, ncol_mo, nrow_mo)
2061 : TYPE(force_env_type), POINTER :: force_env
2062 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2063 : INTEGER, DIMENSION(:) :: ncol_mo, nrow_mo
2064 :
2065 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_wfn_overlap_method'
2066 :
2067 : CHARACTER(LEN=default_path_length) :: file_name
2068 : INTEGER :: handle, ipermutation, ispin, istate, &
2069 : jstate, nao, nforce_eval, nmo, &
2070 : npermutations, nspins
2071 : LOGICAL :: exist, natom_mismatch
2072 : REAL(KIND=dp) :: energy_diff, maxocc, Sda
2073 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_wfn
2074 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: overlaps
2075 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2076 : TYPE(cp_fm_struct_type), POINTER :: mo_mo_fmstruct
2077 : TYPE(cp_fm_type) :: inverse_mat, mo_overlap_wfn, mo_tmp
2078 8 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
2079 : TYPE(cp_logger_type), POINTER :: logger
2080 : TYPE(cp_subsys_type), POINTER :: subsys_mix
2081 : TYPE(dbcsr_type), POINTER :: mixed_matrix_s
2082 8 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_set
2083 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2084 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2085 : TYPE(section_vals_type), POINTER :: force_env_section, mixed_cdft_section
2086 :
2087 8 : NULLIFY (mixed_cdft_section, subsys_mix, particle_set, qs_kind_set, atomic_kind_set, &
2088 8 : mixed_mo_coeff, mixed_matrix_s, force_env_section)
2089 16 : logger => cp_get_default_logger()
2090 :
2091 8 : CALL timeset(routineN, handle)
2092 8 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
2093 8 : npermutations = nforce_eval*(nforce_eval - 1)/2
2094 8 : nspins = SIZE(nrow_mo)
2095 8 : mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
2096 8 : mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
2097 : CALL force_env_get(force_env=force_env, &
2098 8 : force_env_section=force_env_section)
2099 : ! Create mo_set for input wfn
2100 40 : ALLOCATE (mo_set(nspins))
2101 8 : IF (nspins == 2) THEN
2102 8 : maxocc = 1.0_dp
2103 : ELSE
2104 0 : maxocc = 2.0_dp
2105 : END IF
2106 24 : DO ispin = 1, nspins
2107 16 : nao = nrow_mo(ispin)
2108 16 : nmo = ncol_mo(ispin)
2109 : ! Only OT with fully occupied orbitals is implicitly supported
2110 : CALL allocate_mo_set(mo_set(ispin), nao=nao, nmo=nmo, nelectron=INT(maxocc*nmo), &
2111 : n_el_f=REAL(maxocc*nmo, dp), maxocc=maxocc, &
2112 16 : flexible_electron_count=0.0_dp)
2113 16 : CALL set_mo_set(mo_set(ispin), uniform_occupation=.TRUE., homo=nmo)
2114 16 : ALLOCATE (mo_set(ispin)%mo_coeff)
2115 : CALL cp_fm_create(matrix=mo_set(ispin)%mo_coeff, &
2116 : matrix_struct=mixed_mo_coeff(1, ispin)%matrix_struct, &
2117 16 : name="GS_MO_COEFF"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX")
2118 48 : ALLOCATE (mo_set(ispin)%eigenvalues(nmo))
2119 40 : ALLOCATE (mo_set(ispin)%occupation_numbers(nmo))
2120 : END DO
2121 : ! Read wfn file (note we assume that the basis set is the same)
2122 8 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
2123 : ! This really shouldnt be a problem?
2124 : CALL cp_abort(__LOCATION__, &
2125 0 : "QMMM + wavefunction overlap method not supported.")
2126 : END IF
2127 8 : CALL force_env_get(force_env=force_env, subsys=subsys_mix)
2128 8 : mixed_cdft_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT")
2129 8 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
2130 8 : CPASSERT(ASSOCIATED(mixed_cdft%qs_kind_set))
2131 8 : IF (force_env%para_env%is_source()) THEN
2132 4 : CALL wfn_restart_file_name(file_name, exist, mixed_cdft_section, logger)
2133 : END IF
2134 8 : CALL force_env%para_env%bcast(exist)
2135 8 : CALL force_env%para_env%bcast(file_name)
2136 8 : IF (.NOT. exist) THEN
2137 : CALL cp_abort(__LOCATION__, &
2138 : "User requested to restart the wavefunction from the file named: "// &
2139 : TRIM(file_name)//". This file does not exist. Please check the existence of"// &
2140 : " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME in"// &
2141 0 : " section FORCE_EVAL\MIXED\MIXED_CDFT.")
2142 : END IF
2143 : CALL read_mo_set_from_restart(mo_array=mo_set, qs_kind_set=mixed_cdft%qs_kind_set, particle_set=particle_set, &
2144 : para_env=force_env%para_env, id_nr=0, multiplicity=mixed_cdft%multiplicity, &
2145 : dft_section=mixed_cdft_section, natom_mismatch=natom_mismatch, &
2146 8 : cdft=.TRUE.)
2147 8 : IF (natom_mismatch) THEN
2148 : CALL cp_abort(__LOCATION__, &
2149 0 : "Restart wfn file has a wrong number of atoms")
2150 : END IF
2151 : ! Orthonormalize wfn
2152 24 : DO ispin = 1, nspins
2153 24 : IF (mixed_cdft%has_unit_metric) THEN
2154 0 : CALL make_basis_simple(mo_set(ispin)%mo_coeff, ncol_mo(ispin))
2155 : ELSE
2156 16 : CALL make_basis_sm(mo_set(ispin)%mo_coeff, ncol_mo(ispin), mixed_matrix_s)
2157 : END IF
2158 : END DO
2159 : ! Calculate MO overlaps between reference state (R) and CDFT state pairs I/J
2160 24 : ALLOCATE (coupling_wfn(npermutations))
2161 32 : ALLOCATE (overlaps(2, npermutations, nspins))
2162 8 : overlaps = 0.0_dp
2163 24 : DO ispin = 1, nspins
2164 : ! Allocate work
2165 16 : nao = nrow_mo(ispin)
2166 16 : nmo = ncol_mo(ispin)
2167 : CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, &
2168 16 : context=mixed_cdft%blacs_env, para_env=force_env%para_env)
2169 : CALL cp_fm_create(matrix=mo_overlap_wfn, matrix_struct=mo_mo_fmstruct, &
2170 16 : name="MO_OVERLAP_MATRIX_WFN")
2171 : CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
2172 16 : name="INVERSE_MO_OVERLAP_MATRIX_WFN")
2173 16 : CALL cp_fm_struct_release(mo_mo_fmstruct)
2174 : CALL cp_fm_create(matrix=mo_tmp, &
2175 : matrix_struct=mixed_mo_coeff(1, ispin)%matrix_struct, &
2176 16 : name="OVERLAP_MO_COEFF_WFN")
2177 40 : DO ipermutation = 1, npermutations
2178 24 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
2179 : ! S*C_r
2180 : CALL cp_dbcsr_sm_fm_multiply(mixed_matrix_s, mo_set(ispin)%mo_coeff, &
2181 24 : mo_tmp, nmo, 1.0_dp, 0.0_dp)
2182 : ! C_i^T * (S*C_r)
2183 24 : CALL cp_fm_set_all(mo_overlap_wfn, alpha=0.0_dp)
2184 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
2185 : mixed_mo_coeff(istate, ispin), &
2186 24 : mo_tmp, 0.0_dp, mo_overlap_wfn)
2187 24 : CALL cp_fm_invert(mo_overlap_wfn, inverse_mat, overlaps(1, ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
2188 : ! C_j^T * (S*C_r)
2189 24 : CALL cp_fm_set_all(mo_overlap_wfn, alpha=0.0_dp)
2190 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
2191 : mixed_mo_coeff(jstate, ispin), &
2192 24 : mo_tmp, 0.0_dp, mo_overlap_wfn)
2193 64 : CALL cp_fm_invert(mo_overlap_wfn, inverse_mat, overlaps(2, ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
2194 : END DO
2195 16 : CALL cp_fm_release(mo_overlap_wfn)
2196 16 : CALL cp_fm_release(inverse_mat)
2197 16 : CALL cp_fm_release(mo_tmp)
2198 56 : CALL deallocate_mo_set(mo_set(ispin))
2199 : END DO
2200 8 : DEALLOCATE (mo_set)
2201 20 : DO ipermutation = 1, npermutations
2202 12 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
2203 12 : IF (nspins == 2) THEN
2204 12 : overlaps(1, ipermutation, 1) = ABS(overlaps(1, ipermutation, 1)*overlaps(1, ipermutation, 2)) ! A in eq. 12c
2205 12 : overlaps(2, ipermutation, 1) = ABS(overlaps(2, ipermutation, 1)*overlaps(2, ipermutation, 2)) ! B in eq. 12c
2206 : ELSE
2207 0 : overlaps(1, ipermutation, 1) = overlaps(1, ipermutation, 1)**2
2208 0 : overlaps(2, ipermutation, 1) = overlaps(2, ipermutation, 1)**2
2209 : END IF
2210 : ! Calculate coupling using eq. 12c
2211 : ! The coupling is singular if A = B (i.e. states I/J are identical or charge in ground state is fully delocalized)
2212 32 : IF (ABS(overlaps(1, ipermutation, 1) - overlaps(2, ipermutation, 1)) <= 1.0e-14_dp) THEN
2213 : CALL cp_warn(__LOCATION__, &
2214 : "Coupling between states is singular and set to zero. "// &
2215 : "Potential causes: coupling is computed between identical CDFT states or the spin/charge "// &
2216 2 : "density is fully delocalized in the unconstrained ground state.")
2217 2 : coupling_wfn(ipermutation) = 0.0_dp
2218 : ELSE
2219 10 : energy_diff = mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate)
2220 10 : Sda = mixed_cdft%results%S(istate, jstate)
2221 : coupling_wfn(ipermutation) = ABS((overlaps(1, ipermutation, 1)*overlaps(2, ipermutation, 1)/ &
2222 : (overlaps(1, ipermutation, 1)**2 - overlaps(2, ipermutation, 1)**2))* &
2223 : (energy_diff)/(1.0_dp - Sda**2)* &
2224 : (1.0_dp - (overlaps(1, ipermutation, 1)**2 + overlaps(2, ipermutation, 1)**2)/ &
2225 : (2.0_dp*overlaps(1, ipermutation, 1)*overlaps(2, ipermutation, 1))* &
2226 10 : Sda))
2227 : END IF
2228 : END DO
2229 8 : DEALLOCATE (overlaps)
2230 8 : CALL mixed_cdft_result_type_set(mixed_cdft%results, wfn=coupling_wfn)
2231 8 : DEALLOCATE (coupling_wfn)
2232 8 : CALL timestop(handle)
2233 :
2234 16 : END SUBROUTINE mixed_cdft_wfn_overlap_method
2235 :
2236 : ! **************************************************************************************************
2237 : !> \brief Becke constraint adapted to mixed calculations, details in qs_cdft_methods.F
2238 : !> \param force_env the force_env that holds the CDFT states
2239 : !> \param calculate_forces determines if forces should be calculted
2240 : !> \par History
2241 : !> 02.2016 created [Nico Holmberg]
2242 : !> 03.2016 added dynamic load balancing (dlb)
2243 : !> changed pw_p_type data types to rank-3 reals to accommodate dlb
2244 : !> and to reduce overall memory footprint
2245 : !> split to subroutines [Nico Holmberg]
2246 : !> 04.2016 introduced mixed grid mapping [Nico Holmberg]
2247 : ! **************************************************************************************************
2248 36 : SUBROUTINE mixed_becke_constraint(force_env, calculate_forces)
2249 : TYPE(force_env_type), POINTER :: force_env
2250 : LOGICAL, INTENT(IN) :: calculate_forces
2251 :
2252 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint'
2253 :
2254 : INTEGER :: handle
2255 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
2256 : LOGICAL :: in_memory, store_vectors
2257 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint
2258 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients
2259 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: position_vecs, R12
2260 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pair_dist_vecs
2261 : TYPE(cp_logger_type), POINTER :: logger
2262 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2263 : TYPE(mixed_environment_type), POINTER :: mixed_env
2264 :
2265 36 : NULLIFY (mixed_env, mixed_cdft)
2266 36 : store_vectors = .TRUE.
2267 36 : logger => cp_get_default_logger()
2268 36 : CALL timeset(routineN, handle)
2269 36 : mixed_env => force_env%mixed_env
2270 36 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
2271 : CALL mixed_becke_constraint_init(force_env, mixed_cdft, calculate_forces, &
2272 : is_constraint, in_memory, store_vectors, &
2273 : R12, position_vecs, pair_dist_vecs, &
2274 36 : coefficients, catom)
2275 : CALL mixed_becke_constraint_low(force_env, mixed_cdft, in_memory, &
2276 : is_constraint, store_vectors, R12, &
2277 : position_vecs, pair_dist_vecs, &
2278 36 : coefficients, catom)
2279 36 : CALL timestop(handle)
2280 :
2281 36 : END SUBROUTINE mixed_becke_constraint
2282 : ! **************************************************************************************************
2283 : !> \brief Initialize the mixed Becke constraint calculation
2284 : !> \param force_env the force_env that holds the CDFT states
2285 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
2286 : !> \param calculate_forces determines if forces should be calculted
2287 : !> \param is_constraint a list used to determine which atoms in the system define the constraint
2288 : !> \param in_memory decides whether to build the weight function gradients in parallel before solving
2289 : !> the CDFT states or later during the SCF procedure of the individual states
2290 : !> \param store_vectors should temporary arrays be stored in memory to accelerate the calculation
2291 : !> \param R12 temporary array holding the pairwise atomic distances
2292 : !> \param position_vecs temporary array holding the pbc corrected atomic position vectors
2293 : !> \param pair_dist_vecs temporary array holding the pairwise displament vectors
2294 : !> \param coefficients array that determines how atoms should be summed to form the constraint
2295 : !> \param catom temporary array to map the global index of constraint atoms to their position
2296 : !> in a list that holds only constraint atoms
2297 : !> \par History
2298 : !> 03.2016 created [Nico Holmberg]
2299 : ! **************************************************************************************************
2300 36 : SUBROUTINE mixed_becke_constraint_init(force_env, mixed_cdft, calculate_forces, &
2301 : is_constraint, in_memory, store_vectors, &
2302 : R12, position_vecs, pair_dist_vecs, coefficients, &
2303 : catom)
2304 : TYPE(force_env_type), POINTER :: force_env
2305 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2306 : LOGICAL, INTENT(IN) :: calculate_forces
2307 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: is_constraint
2308 : LOGICAL, INTENT(OUT) :: in_memory
2309 : LOGICAL, INTENT(IN) :: store_vectors
2310 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
2311 : INTENT(out) :: R12, position_vecs
2312 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
2313 : INTENT(out) :: pair_dist_vecs
2314 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
2315 : INTENT(OUT) :: coefficients
2316 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(out) :: catom
2317 :
2318 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_init'
2319 :
2320 : CHARACTER(len=2) :: element_symbol
2321 : INTEGER :: atom_a, bounds(2), handle, i, iatom, iex, iforce_eval, ikind, iounit, ithread, j, &
2322 : jatom, katom, my_work, my_work_size, natom, nforce_eval, nkind, np(3), npme, nthread, &
2323 : numexp, offset_dlb, unit_nr
2324 : INTEGER, DIMENSION(2, 3) :: bo, bo_conf
2325 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores, stride
2326 : LOGICAL :: build, mpi_io
2327 : REAL(kind=dp) :: alpha, chi, coef, ircov, jrcov, ra(3), &
2328 : radius, uij
2329 : REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dr, r, r1, shift
2330 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list
2331 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
2332 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2333 : TYPE(cdft_control_type), POINTER :: cdft_control
2334 : TYPE(cell_type), POINTER :: cell
2335 : TYPE(cp_logger_type), POINTER :: logger
2336 : TYPE(cp_subsys_type), POINTER :: subsys_mix
2337 : TYPE(force_env_type), POINTER :: force_env_qs
2338 : TYPE(hirshfeld_type), POINTER :: cavity_env
2339 : TYPE(particle_list_type), POINTER :: particles
2340 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2341 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2342 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2343 : TYPE(realspace_grid_type), POINTER :: rs_cavity
2344 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
2345 :
2346 36 : NULLIFY (pab, cell, force_env_qs, particle_set, force_env_section, print_section, &
2347 36 : qs_kind_set, particles, subsys_mix, rs_cavity, cavity_env, auxbas_pw_pool, &
2348 36 : atomic_kind_set, radii_list, cdft_control)
2349 72 : logger => cp_get_default_logger()
2350 36 : nforce_eval = SIZE(force_env%sub_force_env)
2351 36 : CALL timeset(routineN, handle)
2352 36 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2353 36 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
2354 : CALL force_env_get(force_env=force_env, &
2355 : subsys=subsys_mix, &
2356 28 : cell=cell)
2357 : CALL cp_subsys_get(subsys=subsys_mix, &
2358 : particles=particles, &
2359 28 : particle_set=particle_set)
2360 : ELSE
2361 24 : DO iforce_eval = 1, nforce_eval
2362 16 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
2363 24 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
2364 : END DO
2365 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
2366 : cp_subsys=subsys_mix, &
2367 8 : cell=cell)
2368 : CALL cp_subsys_get(subsys=subsys_mix, &
2369 : particles=particles, &
2370 8 : particle_set=particle_set)
2371 : END IF
2372 36 : natom = SIZE(particles%els)
2373 36 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2374 36 : cdft_control => mixed_cdft%cdft_control
2375 36 : CPASSERT(ASSOCIATED(cdft_control))
2376 36 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cutoffs)) THEN
2377 24 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2378 72 : ALLOCATE (cdft_control%becke_control%cutoffs(natom))
2379 30 : SELECT CASE (cdft_control%becke_control%cutoff_type)
2380 : CASE (becke_cutoff_global)
2381 18 : cdft_control%becke_control%cutoffs(:) = cdft_control%becke_control%rglobal
2382 : CASE (becke_cutoff_element)
2383 18 : IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%becke_control%cutoffs_tmp)) THEN
2384 : CALL cp_abort(__LOCATION__, &
2385 : "Size of keyword BECKE_CONSTRAINT\ELEMENT_CUTOFFS does "// &
2386 0 : "not match number of atomic kinds in the input coordinate file.")
2387 : END IF
2388 54 : DO ikind = 1, SIZE(atomic_kind_set)
2389 36 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2390 90 : DO iatom = 1, katom
2391 36 : atom_a = atom_list(iatom)
2392 72 : cdft_control%becke_control%cutoffs(atom_a) = cdft_control%becke_control%cutoffs_tmp(ikind)
2393 : END DO
2394 : END DO
2395 42 : DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
2396 : END SELECT
2397 : END IF
2398 36 : build = .FALSE.
2399 36 : IF (cdft_control%becke_control%adjust .AND. &
2400 : .NOT. ASSOCIATED(cdft_control%becke_control%aij)) THEN
2401 72 : ALLOCATE (cdft_control%becke_control%aij(natom, natom))
2402 18 : build = .TRUE.
2403 : END IF
2404 108 : ALLOCATE (catom(cdft_control%natoms))
2405 : IF (cdft_control%save_pot .OR. &
2406 : cdft_control%becke_control%cavity_confine .OR. &
2407 36 : cdft_control%becke_control%should_skip .OR. &
2408 : mixed_cdft%first_iteration) THEN
2409 108 : ALLOCATE (is_constraint(natom))
2410 36 : is_constraint = .FALSE.
2411 : END IF
2412 36 : in_memory = calculate_forces .AND. cdft_control%becke_control%in_memory
2413 36 : IF (in_memory .NEQV. calculate_forces) THEN
2414 : CALL cp_abort(__LOCATION__, &
2415 : "The flag BECKE_CONSTRAINT\IN_MEMORY must be activated "// &
2416 0 : "for the calculation of mixed CDFT forces")
2417 : END IF
2418 108 : IF (in_memory .OR. mixed_cdft%first_iteration) ALLOCATE (coefficients(natom))
2419 108 : DO i = 1, cdft_control%natoms
2420 72 : catom(i) = cdft_control%atoms(i)
2421 : IF (cdft_control%save_pot .OR. &
2422 : cdft_control%becke_control%cavity_confine .OR. &
2423 72 : cdft_control%becke_control%should_skip .OR. &
2424 : mixed_cdft%first_iteration) THEN
2425 72 : is_constraint(catom(i)) = .TRUE.
2426 : END IF
2427 108 : IF (in_memory .OR. mixed_cdft%first_iteration) THEN
2428 72 : coefficients(catom(i)) = cdft_control%group(1)%coeff(i)
2429 : END IF
2430 : END DO
2431 36 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
2432 360 : bo = auxbas_pw_pool%pw_grid%bounds_local
2433 144 : np = auxbas_pw_pool%pw_grid%npts
2434 144 : dr = auxbas_pw_pool%pw_grid%dr
2435 144 : shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp
2436 36 : IF (store_vectors) THEN
2437 108 : IF (in_memory) ALLOCATE (pair_dist_vecs(3, natom, natom))
2438 108 : ALLOCATE (position_vecs(3, natom))
2439 : END IF
2440 144 : DO i = 1, 3
2441 144 : cell_v(i) = cell%hmat(i, i)
2442 : END DO
2443 144 : ALLOCATE (R12(natom, natom))
2444 72 : DO iatom = 1, natom - 1
2445 108 : DO jatom = iatom + 1, natom
2446 144 : r = particle_set(iatom)%r
2447 144 : r1 = particle_set(jatom)%r
2448 144 : DO i = 1, 3
2449 108 : r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2450 144 : r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2451 : END DO
2452 144 : dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
2453 36 : IF (store_vectors) THEN
2454 144 : position_vecs(:, iatom) = r(:)
2455 144 : IF (iatom == 1 .AND. jatom == natom) position_vecs(:, jatom) = r1(:)
2456 36 : IF (in_memory) THEN
2457 96 : pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
2458 96 : pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
2459 : END IF
2460 : END IF
2461 144 : R12(iatom, jatom) = NORM2(dist_vec)
2462 36 : R12(jatom, iatom) = R12(iatom, jatom)
2463 72 : IF (build) THEN
2464 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2465 18 : kind_number=ikind)
2466 18 : ircov = cdft_control%becke_control%radii(ikind)
2467 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
2468 18 : kind_number=ikind)
2469 18 : jrcov = cdft_control%becke_control%radii(ikind)
2470 18 : IF (ircov /= jrcov) THEN
2471 18 : chi = ircov/jrcov
2472 18 : uij = (chi - 1.0_dp)/(chi + 1.0_dp)
2473 18 : cdft_control%becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp)
2474 18 : IF (cdft_control%becke_control%aij(iatom, jatom) &
2475 : > 0.5_dp) THEN
2476 0 : cdft_control%becke_control%aij(iatom, jatom) = 0.5_dp
2477 18 : ELSE IF (cdft_control%becke_control%aij(iatom, jatom) &
2478 : < -0.5_dp) THEN
2479 0 : cdft_control%becke_control%aij(iatom, jatom) = -0.5_dp
2480 : END IF
2481 : ELSE
2482 0 : cdft_control%becke_control%aij(iatom, jatom) = 0.0_dp
2483 : END IF
2484 : cdft_control%becke_control%aij(jatom, iatom) = &
2485 18 : -cdft_control%becke_control%aij(iatom, jatom)
2486 : END IF
2487 : END DO
2488 : END DO
2489 : ! Dump some additional information about the calculation
2490 36 : IF (mixed_cdft%first_iteration) THEN
2491 24 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2492 24 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2493 24 : IF (iounit > 0) THEN
2494 : WRITE (iounit, '(/,T3,A,T66)') &
2495 12 : '-------------------------- Becke atomic parameters ---------------------------'
2496 12 : IF (cdft_control%becke_control%adjust) THEN
2497 : WRITE (iounit, '(T3,A,A)') &
2498 9 : 'Atom Element Coefficient', ' Cutoff (angstrom) CDFT Radius (angstrom)'
2499 27 : DO iatom = 1, natom
2500 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2501 : element_symbol=element_symbol, &
2502 18 : kind_number=ikind)
2503 18 : ircov = cp_unit_from_cp2k(cdft_control%becke_control%radii(ikind), "angstrom")
2504 18 : IF (is_constraint(iatom)) THEN
2505 18 : coef = coefficients(iatom)
2506 : ELSE
2507 0 : coef = 0.0_dp
2508 : END IF
2509 : WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3,T73,F8.3)") &
2510 18 : iatom, ADJUSTR(element_symbol), coef, &
2511 18 : cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom"), &
2512 63 : ircov
2513 : END DO
2514 : ELSE
2515 : WRITE (iounit, '(T3,A,A)') &
2516 3 : 'Atom Element Coefficient', ' Cutoff (angstrom)'
2517 9 : DO iatom = 1, natom
2518 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2519 6 : element_symbol=element_symbol)
2520 6 : IF (is_constraint(iatom)) THEN
2521 6 : coef = coefficients(iatom)
2522 : ELSE
2523 0 : coef = 0.0_dp
2524 : END IF
2525 : WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3)") &
2526 6 : iatom, ADJUSTR(element_symbol), coef, &
2527 15 : cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom")
2528 : END DO
2529 : END IF
2530 : WRITE (iounit, '(T3,A)') &
2531 12 : '------------------------------------------------------------------------------'
2532 : END IF
2533 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
2534 24 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2535 24 : mixed_cdft%first_iteration = .FALSE.
2536 : END IF
2537 :
2538 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2539 34 : CPASSERT(ASSOCIATED(mixed_cdft%qs_kind_set))
2540 34 : cavity_env => cdft_control%becke_control%cavity_env
2541 34 : qs_kind_set => mixed_cdft%qs_kind_set
2542 34 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2543 34 : nkind = SIZE(qs_kind_set)
2544 34 : IF (.NOT. ASSOCIATED(cavity_env%kind_shape_fn)) THEN
2545 22 : IF (ASSOCIATED(cdft_control%becke_control%radii)) THEN
2546 54 : ALLOCATE (radii_list(SIZE(cdft_control%becke_control%radii)))
2547 54 : DO ikind = 1, SIZE(cdft_control%becke_control%radii)
2548 54 : IF (cavity_env%use_bohr) THEN
2549 0 : radii_list(ikind) = cdft_control%becke_control%radii(ikind)
2550 : ELSE
2551 36 : radii_list(ikind) = cp_unit_from_cp2k(cdft_control%becke_control%radii(ikind), "angstrom")
2552 : END IF
2553 : END DO
2554 : END IF
2555 : CALL create_shape_function(cavity_env, qs_kind_set, atomic_kind_set, &
2556 : radius=cdft_control%becke_control%rcavity, &
2557 22 : radii_list=radii_list)
2558 22 : IF (ASSOCIATED(radii_list)) THEN
2559 18 : DEALLOCATE (radii_list)
2560 : END IF
2561 : END IF
2562 34 : NULLIFY (rs_cavity)
2563 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_rs_grid=rs_cavity, &
2564 34 : auxbas_pw_pool=auxbas_pw_pool)
2565 : ! be careful in parallel nsmax is chosen with multigrid in mind!
2566 34 : CALL rs_grid_zero(rs_cavity)
2567 34 : ALLOCATE (pab(1, 1))
2568 34 : nthread = 1
2569 34 : ithread = 0
2570 102 : DO ikind = 1, SIZE(atomic_kind_set)
2571 68 : numexp = cavity_env%kind_shape_fn(ikind)%numexp
2572 68 : IF (numexp <= 0) CYCLE
2573 68 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2574 204 : ALLOCATE (cores(katom))
2575 136 : DO iex = 1, numexp
2576 68 : alpha = cavity_env%kind_shape_fn(ikind)%zet(iex)
2577 68 : coef = cavity_env%kind_shape_fn(ikind)%coef(iex)
2578 68 : npme = 0
2579 136 : cores = 0
2580 136 : DO iatom = 1, katom
2581 136 : IF (rs_cavity%desc%parallel .AND. .NOT. rs_cavity%desc%distributed) THEN
2582 : ! replicated realspace grid, split the atoms up between procs
2583 68 : IF (MODULO(iatom, rs_cavity%desc%group_size) == rs_cavity%desc%my_pos) THEN
2584 34 : npme = npme + 1
2585 34 : cores(npme) = iatom
2586 : END IF
2587 : ELSE
2588 0 : npme = npme + 1
2589 0 : cores(npme) = iatom
2590 : END IF
2591 : END DO
2592 170 : DO j = 1, npme
2593 34 : iatom = cores(j)
2594 34 : atom_a = atom_list(iatom)
2595 34 : pab(1, 1) = coef
2596 34 : IF (store_vectors) THEN
2597 136 : ra(:) = position_vecs(:, atom_a) + cell_v(:)/2._dp
2598 : ELSE
2599 0 : ra(:) = pbc(particle_set(atom_a)%r, cell)
2600 : END IF
2601 102 : IF (is_constraint(atom_a)) THEN
2602 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
2603 : ra=ra, rb=ra, rp=ra, &
2604 : zetp=alpha, eps=mixed_cdft%eps_rho_rspace, &
2605 : pab=pab, o1=0, o2=0, & ! without map_consistent
2606 34 : prefactor=1.0_dp, cutoff=0.0_dp)
2607 :
2608 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
2609 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
2610 : rs_cavity, &
2611 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
2612 : use_subpatch=.TRUE., &
2613 34 : subpatch_pattern=0)
2614 : END IF
2615 : END DO
2616 : END DO
2617 170 : DEALLOCATE (cores)
2618 : END DO
2619 34 : DEALLOCATE (pab)
2620 34 : CALL auxbas_pw_pool%create_pw(cdft_control%becke_control%cavity)
2621 34 : CALL transfer_rs2pw(rs_cavity, cdft_control%becke_control%cavity)
2622 : CALL hfun_zero(cdft_control%becke_control%cavity%array, &
2623 : cdft_control%becke_control%eps_cavity, &
2624 34 : just_zero=.FALSE., bounds=bounds, work=my_work)
2625 34 : IF (bounds(2) < bo(2, 3)) THEN
2626 8 : bounds(2) = bounds(2) - 1
2627 : ELSE
2628 26 : bounds(2) = bo(2, 3)
2629 : END IF
2630 34 : IF (bounds(1) > bo(1, 3)) THEN
2631 : ! In the special case bounds(1) == bounds(2) == bo(2, 3), after this check
2632 : ! bounds(1) > bounds(2) and the subsequent gradient allocation (:, :, :, bounds(1):bounds(2))
2633 : ! will correctly allocate a 0-sized array
2634 8 : bounds(1) = bounds(1) + 1
2635 : ELSE
2636 26 : bounds(1) = bo(1, 3)
2637 : END IF
2638 34 : IF (bounds(1) > bounds(2)) THEN
2639 0 : my_work_size = 0
2640 : ELSE
2641 34 : my_work_size = (bounds(2) - bounds(1) + 1)
2642 34 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2643 0 : my_work_size = my_work_size*(bo(2, 2) - bo(1, 2) + 1)
2644 : ELSE
2645 34 : my_work_size = my_work_size*(bo(2, 1) - bo(1, 1) + 1)
2646 : END IF
2647 : END IF
2648 102 : cdft_control%becke_control%confine_bounds = bounds
2649 34 : IF (cdft_control%becke_control%print_cavity) THEN
2650 : CALL hfun_zero(cdft_control%becke_control%cavity%array, &
2651 0 : cdft_control%becke_control%eps_cavity, just_zero=.TRUE.)
2652 : NULLIFY (stride)
2653 0 : ALLOCATE (stride(3))
2654 0 : stride = [2, 2, 2]
2655 0 : mpi_io = .TRUE.
2656 : unit_nr = cp_print_key_unit_nr(logger, print_section, "", &
2657 : middle_name="BECKE_CAVITY", &
2658 : extension=".cube", file_position="REWIND", &
2659 0 : log_filename=.FALSE., mpi_io=mpi_io)
2660 0 : IF (force_env%para_env%is_source() .AND. unit_nr < 1) THEN
2661 : CALL cp_abort(__LOCATION__, &
2662 0 : "Please turn on PROGRAM_RUN_INFO to print cavity")
2663 : END IF
2664 : CALL cp_pw_to_cube(cdft_control%becke_control%cavity, &
2665 : unit_nr, "CAVITY", particles=particles, &
2666 0 : stride=stride, mpi_io=mpi_io)
2667 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', mpi_io=mpi_io)
2668 0 : DEALLOCATE (stride)
2669 : END IF
2670 : END IF
2671 36 : bo_conf = bo
2672 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2673 102 : bo_conf(:, 3) = cdft_control%becke_control%confine_bounds
2674 : END IF
2675 : ! Load balance
2676 36 : IF (mixed_cdft%dlb) THEN
2677 : CALL mixed_becke_constraint_dlb(force_env, mixed_cdft, my_work, &
2678 8 : my_work_size, natom, bo, bo_conf)
2679 : END IF
2680 : ! The bounds have been finalized => time to allocate storage for working matrices
2681 36 : offset_dlb = 0
2682 36 : IF (mixed_cdft%dlb) THEN
2683 8 : IF (mixed_cdft%dlb_control%send_work .AND. .NOT. mixed_cdft%is_special) THEN
2684 8 : offset_dlb = SUM(mixed_cdft%dlb_control%target_list(2, :))
2685 : END IF
2686 : END IF
2687 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2688 : ! Get rid of the zero part of the confinement cavity (cr3d -> real(:,:,:))
2689 34 : IF (mixed_cdft%is_special) THEN
2690 0 : ALLOCATE (mixed_cdft%sendbuff(SIZE(mixed_cdft%dest_list)))
2691 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2692 : ALLOCATE (mixed_cdft%sendbuff(i)%cavity(mixed_cdft%dest_list_bo(1, i):mixed_cdft%dest_list_bo(2, i), &
2693 0 : bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2694 : mixed_cdft%sendbuff(i)%cavity = cdft_control%becke_control%cavity%array(mixed_cdft%dest_list_bo(1, i): &
2695 : mixed_cdft%dest_list_bo(2, i), &
2696 : bo(1, 2):bo(2, 2), &
2697 0 : bo_conf(1, 3):bo_conf(2, 3))
2698 : END DO
2699 34 : ELSE IF (mixed_cdft%is_pencil) THEN
2700 0 : ALLOCATE (mixed_cdft%cavity(bo(1, 1) + offset_dlb:bo(2, 1), bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2701 : mixed_cdft%cavity = cdft_control%becke_control%cavity%array(bo(1, 1) + offset_dlb:bo(2, 1), &
2702 : bo(1, 2):bo(2, 2), &
2703 0 : bo_conf(1, 3):bo_conf(2, 3))
2704 : ELSE
2705 170 : ALLOCATE (mixed_cdft%cavity(bo(1, 1):bo(2, 1), bo(1, 2) + offset_dlb:bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2706 : mixed_cdft%cavity = cdft_control%becke_control%cavity%array(bo(1, 1):bo(2, 1), &
2707 : bo(1, 2) + offset_dlb:bo(2, 2), &
2708 3568226 : bo_conf(1, 3):bo_conf(2, 3))
2709 : END IF
2710 34 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
2711 : END IF
2712 36 : IF (mixed_cdft%is_special) THEN
2713 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2714 : ALLOCATE (mixed_cdft%sendbuff(i)%weight(mixed_cdft%dest_list_bo(1, i):mixed_cdft%dest_list_bo(2, i), &
2715 0 : bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2716 0 : mixed_cdft%sendbuff(i)%weight = 0.0_dp
2717 : END DO
2718 36 : ELSE IF (mixed_cdft%is_pencil) THEN
2719 0 : ALLOCATE (mixed_cdft%weight(bo(1, 1) + offset_dlb:bo(2, 1), bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2720 0 : mixed_cdft%weight = 0.0_dp
2721 : ELSE
2722 180 : ALLOCATE (mixed_cdft%weight(bo(1, 1):bo(2, 1), bo(1, 2) + offset_dlb:bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2723 1879397 : mixed_cdft%weight = 0.0_dp
2724 : END IF
2725 36 : IF (in_memory) THEN
2726 24 : IF (mixed_cdft%is_special) THEN
2727 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2728 : ALLOCATE (mixed_cdft%sendbuff(i)%gradients(3*natom, mixed_cdft%dest_list_bo(1, i): &
2729 : mixed_cdft%dest_list_bo(2, i), &
2730 : bo(1, 2):bo(2, 2), &
2731 0 : bo_conf(1, 3):bo_conf(2, 3)))
2732 0 : mixed_cdft%sendbuff(i)%gradients = 0.0_dp
2733 : END DO
2734 24 : ELSE IF (mixed_cdft%is_pencil) THEN
2735 : ALLOCATE (cdft_control%group(1)%gradients(3*natom, bo(1, 1) + offset_dlb:bo(2, 1), &
2736 : bo(1, 2):bo(2, 2), &
2737 0 : bo_conf(1, 3):bo_conf(2, 3)))
2738 0 : cdft_control%group(1)%gradients = 0.0_dp
2739 : ELSE
2740 : ALLOCATE (cdft_control%group(1)%gradients(3*natom, bo(1, 1):bo(2, 1), &
2741 : bo(1, 2) + offset_dlb:bo(2, 2), &
2742 144 : bo_conf(1, 3):bo_conf(2, 3)))
2743 9808936 : cdft_control%group(1)%gradients = 0.0_dp
2744 : END IF
2745 : END IF
2746 :
2747 36 : CALL timestop(handle)
2748 :
2749 36 : END SUBROUTINE mixed_becke_constraint_init
2750 :
2751 : ! **************************************************************************************************
2752 : !> \brief Setup load balancing for mixed Becke calculation
2753 : !> \param force_env the force_env that holds the CDFT states
2754 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
2755 : !> \param my_work an estimate of the work per processor
2756 : !> \param my_work_size size of the smallest array slice per processor. overloaded processors will
2757 : !> redistribute works as integer multiples of this value.
2758 : !> \param natom the total number of atoms
2759 : !> \param bo bounds of the realspace grid that holds the electron density
2760 : !> \param bo_conf same as bo, but bounds along z-direction have been compacted with confinement
2761 : !> \par History
2762 : !> 03.2016 created [Nico Holmberg]
2763 : ! **************************************************************************************************
2764 8 : SUBROUTINE mixed_becke_constraint_dlb(force_env, mixed_cdft, my_work, &
2765 : my_work_size, natom, bo, bo_conf)
2766 : TYPE(force_env_type), POINTER :: force_env
2767 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2768 : INTEGER, INTENT(IN) :: my_work, my_work_size, natom
2769 : INTEGER, DIMENSION(2, 3) :: bo, bo_conf
2770 :
2771 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_dlb'
2772 : INTEGER, PARAMETER :: should_deallocate = 7000, &
2773 : uninitialized = -7000
2774 :
2775 : CHARACTER(len=2) :: dummy
2776 : INTEGER :: actually_sent, exhausted_work, handle, i, ind, iounit, ispecial, j, max_targets, &
2777 : more_work, my_pos, my_special_work, my_target, no_overloaded, no_underloaded, nsend, &
2778 : nsend_limit, nsend_max, offset, offset_proc, offset_special, send_total, tags(2)
2779 8 : INTEGER, DIMENSION(:), POINTER :: buffsize, cumulative_work, expected_work, load_imbalance, &
2780 16 : nrecv, nsend_proc, sendbuffer, should_warn, tmp, work_index, work_size
2781 8 : INTEGER, DIMENSION(:, :), POINTER :: targets, tmp_bo
2782 : LOGICAL :: consistent
2783 16 : LOGICAL, DIMENSION(:), POINTER :: mask_recv, mask_send, touched
2784 : REAL(kind=dp) :: average_work, load_scale, &
2785 : very_overloaded, work_factor
2786 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cavity
2787 16 : TYPE(buffers_bi), DIMENSION(:), POINTER :: recvbuffer, sbuff
2788 : TYPE(cdft_control_type), POINTER :: cdft_control
2789 : TYPE(cp_logger_type), POINTER :: logger
2790 40 : TYPE(mp_request_type), DIMENSION(4) :: req
2791 8 : TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
2792 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
2793 :
2794 16 : logger => cp_get_default_logger()
2795 8 : CALL timeset(routineN, handle)
2796 8 : mixed_cdft%dlb_control%recv_work = .FALSE.
2797 8 : mixed_cdft%dlb_control%send_work = .FALSE.
2798 8 : NULLIFY (expected_work, work_index, load_imbalance, work_size, &
2799 8 : cumulative_work, sendbuffer, buffsize, req_recv, req_total, &
2800 8 : tmp, nrecv, nsend_proc, targets, tmp_bo, touched, &
2801 8 : mask_recv, mask_send, cavity, recvbuffer, sbuff, force_env_section, &
2802 8 : print_section, cdft_control)
2803 8 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2804 8 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2805 8 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2806 8 : cdft_control => mixed_cdft%cdft_control
2807 : ! These numerical values control data redistribution and are system sensitive
2808 : ! Currently they are not refined during run time which may cause crashes
2809 : ! However, using too many processors or a confinement cavity that is too large relative to the
2810 : ! total system volume are more likely culprits.
2811 8 : load_scale = mixed_cdft%dlb_control%load_scale
2812 8 : very_overloaded = mixed_cdft%dlb_control%very_overloaded
2813 8 : more_work = mixed_cdft%dlb_control%more_work
2814 8 : max_targets = 40
2815 8 : work_factor = 0.8_dp
2816 : ! Reset targets/sources
2817 8 : IF (mixed_cdft%is_special) THEN
2818 0 : DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo, &
2819 0 : mixed_cdft%source_list, mixed_cdft%source_list_bo)
2820 : ALLOCATE (mixed_cdft%dest_list(SIZE(mixed_cdft%dest_list_save)), &
2821 : mixed_cdft%dest_list_bo(SIZE(mixed_cdft%dest_bo_save, 1), SIZE(mixed_cdft%dest_bo_save, 2)), &
2822 : mixed_cdft%source_list(SIZE(mixed_cdft%source_list_save)), &
2823 0 : mixed_cdft%source_list_bo(SIZE(mixed_cdft%source_bo_save, 1), SIZE(mixed_cdft%source_bo_save, 2)))
2824 0 : mixed_cdft%dest_list = mixed_cdft%dest_list_save
2825 0 : mixed_cdft%source_list = mixed_cdft%source_list_save
2826 0 : mixed_cdft%dest_list_bo = mixed_cdft%dest_bo_save
2827 0 : mixed_cdft%source_list_bo = mixed_cdft%source_bo_save
2828 : END IF
2829 : ALLOCATE (mixed_cdft%dlb_control%expected_work(force_env%para_env%num_pe), &
2830 : expected_work(force_env%para_env%num_pe), &
2831 48 : work_size(force_env%para_env%num_pe))
2832 : IF (debug_this_module) THEN
2833 : ALLOCATE (should_warn(force_env%para_env%num_pe))
2834 : should_warn = 0
2835 : END IF
2836 24 : expected_work = 0
2837 8 : expected_work(force_env%para_env%mepos + 1) = my_work
2838 24 : work_size = 0
2839 8 : work_size(force_env%para_env%mepos + 1) = my_work_size
2840 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2841 4 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2842 : work_size(force_env%para_env%mepos + 1) = work_size(force_env%para_env%mepos + 1) - &
2843 : NINT(REAL(mixed_cdft%dlb_control% &
2844 : prediction_error(force_env%para_env%mepos + 1), dp)/ &
2845 0 : REAL(bo(2, 1) - bo(1, 1) + 1, dp))
2846 : ELSE
2847 : work_size(force_env%para_env%mepos + 1) = work_size(force_env%para_env%mepos + 1) - &
2848 : NINT(REAL(mixed_cdft%dlb_control% &
2849 : prediction_error(force_env%para_env%mepos + 1), dp)/ &
2850 4 : REAL(bo(2, 2) - bo(1, 2) + 1, dp))
2851 : END IF
2852 : END IF
2853 40 : CALL force_env%para_env%sum(expected_work)
2854 40 : CALL force_env%para_env%sum(work_size)
2855 : ! We store the unsorted expected work to refine the estimate on subsequent calls to this routine
2856 40 : mixed_cdft%dlb_control%expected_work = expected_work
2857 : ! Take into account the prediction error of the last step
2858 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2859 20 : expected_work = expected_work - mixed_cdft%dlb_control%prediction_error
2860 : END IF
2861 : !
2862 24 : average_work = REAL(SUM(expected_work), dp)/REAL(force_env%para_env%num_pe, dp)
2863 : ALLOCATE (work_index(force_env%para_env%num_pe), &
2864 : load_imbalance(force_env%para_env%num_pe), &
2865 48 : targets(2, force_env%para_env%num_pe))
2866 40 : load_imbalance = expected_work - NINT(average_work)
2867 8 : no_overloaded = 0
2868 8 : no_underloaded = 0
2869 56 : targets = 0
2870 : ! Convert the load imbalance to a multiple of the actual work size
2871 24 : DO i = 1, force_env%para_env%num_pe
2872 24 : IF (load_imbalance(i) > 0) THEN
2873 8 : no_overloaded = no_overloaded + 1
2874 : ! Allow heavily overloaded processors to dump more data since most likely they have a lot of 'real' work
2875 8 : IF (expected_work(i) > NINT(very_overloaded*average_work)) THEN
2876 0 : load_imbalance(i) = (CEILING(REAL(load_imbalance(i), dp)/REAL(work_size(i), dp)) + more_work)*work_size(i)
2877 : ELSE
2878 8 : load_imbalance(i) = CEILING(REAL(load_imbalance(i), dp)/REAL(work_size(i), dp))*work_size(i)
2879 : END IF
2880 : ELSE
2881 : ! Allow the underloaded processors to take load_scale amount of additional work
2882 : ! otherwise we may be unable to exhaust all overloaded processors
2883 8 : load_imbalance(i) = NINT(load_imbalance(i)*load_scale)
2884 8 : no_underloaded = no_underloaded + 1
2885 : END IF
2886 : END DO
2887 8 : CALL sort(expected_work, force_env%para_env%num_pe, indices=work_index)
2888 : ! Redistribute work in order from the most overloaded processors to the most underloaded processors
2889 : ! Each underloaded processor is limited to one overloaded processor
2890 8 : IF (load_imbalance(force_env%para_env%mepos + 1) > 0) THEN
2891 4 : offset = 0
2892 4 : mixed_cdft%dlb_control%send_work = .TRUE.
2893 : ! Build up the total amount of work that needs redistribution
2894 12 : ALLOCATE (cumulative_work(force_env%para_env%num_pe))
2895 12 : cumulative_work = 0
2896 4 : DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
2897 4 : IF (work_index(i) == force_env%para_env%mepos + 1) THEN
2898 : EXIT
2899 : ELSE
2900 0 : offset = offset + load_imbalance(work_index(i))
2901 0 : IF (i == force_env%para_env%num_pe) THEN
2902 0 : cumulative_work(i) = load_imbalance(work_index(i))
2903 : ELSE
2904 0 : cumulative_work(i) = cumulative_work(i + 1) + load_imbalance(work_index(i))
2905 : END IF
2906 : END IF
2907 : END DO
2908 4 : my_pos = i
2909 4 : j = force_env%para_env%num_pe
2910 4 : nsend_max = load_imbalance(work_index(j))/work_size(work_index(j))
2911 4 : exhausted_work = 0
2912 : ! Determine send offset by going through all processors that are more overloaded than my_pos
2913 4 : DO i = 1, no_underloaded
2914 4 : IF (my_pos == force_env%para_env%num_pe) EXIT
2915 0 : nsend = -load_imbalance(work_index(i))/work_size(work_index(j))
2916 0 : IF (nsend < 1) nsend = 1
2917 0 : nsend_max = nsend_max - nsend
2918 0 : IF (nsend_max < 0) nsend = nsend + nsend_max
2919 0 : exhausted_work = exhausted_work + nsend*work_size(work_index(j))
2920 0 : offset = offset - nsend*work_size(work_index(j))
2921 0 : IF (offset < 0) EXIT
2922 4 : IF (exhausted_work == cumulative_work(j)) THEN
2923 0 : j = j - 1
2924 0 : nsend_max = load_imbalance(work_index(j))/work_size(work_index(j))
2925 : END IF
2926 : END DO
2927 : ! Underloaded processors were fully exhausted: rewind index
2928 : ! Load balancing will fail if this happens on multiple processors
2929 4 : IF (i > no_underloaded) THEN
2930 0 : i = no_underloaded
2931 : END IF
2932 4 : my_target = i
2933 4 : DEALLOCATE (cumulative_work)
2934 : ! Determine how much and who to send slices of my grid points
2935 4 : nsend_max = load_imbalance(force_env%para_env%mepos + 1)/work_size(force_env%para_env%mepos + 1)
2936 : ! This the actual number of available array slices
2937 4 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2938 0 : nsend_limit = bo(2, 1) - bo(1, 1) + 1
2939 : ELSE
2940 4 : nsend_limit = bo(2, 2) - bo(1, 2) + 1
2941 : END IF
2942 4 : IF (.NOT. mixed_cdft%is_special) THEN
2943 4 : ALLOCATE (mixed_cdft%dlb_control%target_list(3, max_targets))
2944 : ELSE
2945 0 : ALLOCATE (mixed_cdft%dlb_control%target_list(3 + 2*SIZE(mixed_cdft%dest_list), max_targets))
2946 0 : ALLOCATE (touched(SIZE(mixed_cdft%dest_list)))
2947 0 : touched = .FALSE.
2948 : END IF
2949 644 : mixed_cdft%dlb_control%target_list = uninitialized
2950 4 : i = 1
2951 4 : ispecial = 1
2952 4 : offset_special = 0
2953 4 : targets(1, my_pos) = my_target
2954 4 : send_total = 0
2955 : ! Main loop. Note, we actually allow my_pos to offload more slices than nsend_max
2956 : DO
2957 4 : nsend = -load_imbalance(work_index(my_target))/work_size(force_env%para_env%mepos + 1)
2958 4 : IF (nsend < 1) nsend = 1 ! send at least one block
2959 : ! Prevent over redistribution: leave at least (1-work_factor)*nsend_limit slices to my_pos
2960 4 : IF (nsend > NINT(work_factor*nsend_limit - send_total)) THEN
2961 : nsend = NINT(work_factor*nsend_limit - send_total)
2962 : IF (debug_this_module) THEN
2963 : should_warn(force_env%para_env%mepos + 1) = 1
2964 : END IF
2965 : END IF
2966 4 : mixed_cdft%dlb_control%target_list(1, i) = work_index(my_target) - 1 ! This is the actual processor rank
2967 4 : IF (mixed_cdft%is_special) THEN
2968 0 : mixed_cdft%dlb_control%target_list(2, i) = 0
2969 0 : actually_sent = nsend
2970 0 : DO j = ispecial, SIZE(mixed_cdft%dest_list)
2971 0 : mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + 1
2972 0 : touched(j) = .TRUE.
2973 0 : IF (nsend < mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1) THEN
2974 0 : mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2975 0 : mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(1, j) + nsend - 1
2976 0 : mixed_cdft%dest_list_bo(1, j) = mixed_cdft%dest_list_bo(1, j) + nsend
2977 0 : nsend = 0
2978 0 : EXIT
2979 : ELSE
2980 0 : mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2981 0 : mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(2, j)
2982 0 : nsend = nsend - (mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1)
2983 0 : mixed_cdft%dest_list_bo(1:2, j) = should_deallocate
2984 : END IF
2985 0 : IF (nsend <= 0) EXIT
2986 : END DO
2987 0 : IF (mixed_cdft%dest_list_bo(1, ispecial) == should_deallocate) ispecial = j + 1
2988 0 : actually_sent = actually_sent - nsend
2989 0 : nsend_max = nsend_max - actually_sent
2990 0 : send_total = send_total + actually_sent
2991 : ELSE
2992 4 : mixed_cdft%dlb_control%target_list(2, i) = nsend
2993 4 : nsend_max = nsend_max - nsend
2994 4 : send_total = send_total + nsend
2995 : END IF
2996 4 : IF (nsend_max < 0) nsend_max = 0
2997 4 : IF (nsend_max == 0) EXIT
2998 0 : IF (my_target /= no_underloaded) THEN
2999 0 : my_target = my_target + 1
3000 : ELSE
3001 : ! If multiple processors execute this block load balancing will fail
3002 0 : mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + nsend_max
3003 0 : nsend_max = 0
3004 0 : EXIT
3005 : END IF
3006 0 : i = i + 1
3007 4 : IF (i > max_targets) THEN
3008 : CALL cp_abort(__LOCATION__, &
3009 0 : "Load balancing error: increase max_targets")
3010 : END IF
3011 : END DO
3012 4 : IF (.NOT. mixed_cdft%is_special) THEN
3013 4 : CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3, 1, i)
3014 : ELSE
3015 0 : CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3 + 2*SIZE(mixed_cdft%dest_list), 1, i)
3016 : END IF
3017 4 : targets(2, my_pos) = my_target
3018 : ! Equalize the load on the target processors
3019 4 : IF (.NOT. mixed_cdft%is_special) THEN
3020 4 : IF (send_total > NINT(work_factor*nsend_limit)) send_total = NINT(work_factor*nsend_limit) - 1
3021 4 : nsend = NINT(REAL(send_total, dp)/REAL(SIZE(mixed_cdft%dlb_control%target_list, 2), dp))
3022 8 : mixed_cdft%dlb_control%target_list(2, :) = nsend
3023 : END IF
3024 : ELSE
3025 4 : DO i = 1, no_underloaded
3026 4 : IF (work_index(i) == force_env%para_env%mepos + 1) EXIT
3027 : END DO
3028 : my_pos = i
3029 : END IF
3030 104 : CALL force_env%para_env%sum(targets)
3031 : IF (debug_this_module) THEN
3032 : CALL force_env%para_env%sum(should_warn)
3033 : IF (ANY(should_warn == 1)) THEN
3034 : CALL cp_warn(__LOCATION__, &
3035 : "MIXED_CDFT DLB: Attempted to redistribute more array"// &
3036 : " slices than actually available. Leaving a fraction of the total"// &
3037 : " slices on the overloaded processor. Perhaps you have set LOAD_SCALE too high?")
3038 : END IF
3039 : DEALLOCATE (should_warn)
3040 : END IF
3041 : ! check that there is one-to-one mapping between over- and underloaded processors
3042 8 : IF (force_env%para_env%is_source()) THEN
3043 4 : consistent = .TRUE.
3044 4 : DO i = force_env%para_env%num_pe - 1, force_env%para_env%num_pe - no_overloaded + 1, -1
3045 0 : IF (targets(1, i) > no_underloaded) consistent = .FALSE.
3046 4 : IF (targets(1, i) > targets(2, i + 1)) THEN
3047 : CYCLE
3048 : ELSE
3049 0 : consistent = .FALSE.
3050 : END IF
3051 : END DO
3052 4 : IF (.NOT. consistent) THEN
3053 : IF (debug_this_module .AND. iounit > 0) THEN
3054 : DO i = force_env%para_env%num_pe - 1, force_env%para_env%num_pe - no_overloaded + 1, -1
3055 : WRITE (iounit, '(A,I8,I8,I8,I8,I8)') &
3056 : 'load balancing info', load_imbalance(i), work_index(i), &
3057 : work_size(i), targets(1, i), targets(2, i)
3058 : END DO
3059 : END IF
3060 : CALL cp_abort(__LOCATION__, &
3061 : "Load balancing error: too much data to redistribute."// &
3062 : " Increase LOAD_SCALE or change the number of processors."// &
3063 : " If the confinement cavity occupies a large volume relative"// &
3064 0 : " to the total system volume, it might be worth disabling DLB.")
3065 : END IF
3066 : END IF
3067 : ! Tell the target processors which grid points they should compute
3068 8 : IF (my_pos <= no_underloaded) THEN
3069 4 : DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
3070 4 : IF (targets(1, i) <= my_pos .AND. targets(2, i) >= my_pos) THEN
3071 4 : mixed_cdft%dlb_control%recv_work = .TRUE.
3072 4 : mixed_cdft%dlb_control%my_source = work_index(i) - 1
3073 4 : EXIT
3074 : END IF
3075 : END DO
3076 4 : IF (mixed_cdft%dlb_control%recv_work) THEN
3077 4 : IF (.NOT. mixed_cdft%is_special) THEN
3078 4 : ALLOCATE (mixed_cdft%dlb_control%bo(12))
3079 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%bo, source=mixed_cdft%dlb_control%my_source, &
3080 4 : request=req(1))
3081 4 : CALL req(1)%wait()
3082 12 : mixed_cdft%dlb_control%my_dest_repl = [mixed_cdft%dlb_control%bo(11), mixed_cdft%dlb_control%bo(12)]
3083 12 : mixed_cdft%dlb_control%dest_tags_repl = [mixed_cdft%dlb_control%bo(9), mixed_cdft%dlb_control%bo(10)]
3084 : ALLOCATE (mixed_cdft%dlb_control%cavity(mixed_cdft%dlb_control%bo(1):mixed_cdft%dlb_control%bo(2), &
3085 : mixed_cdft%dlb_control%bo(3):mixed_cdft%dlb_control%bo(4), &
3086 20 : mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3087 : ALLOCATE (mixed_cdft%dlb_control%weight(mixed_cdft%dlb_control%bo(1):mixed_cdft%dlb_control%bo(2), &
3088 : mixed_cdft%dlb_control%bo(3):mixed_cdft%dlb_control%bo(4), &
3089 20 : mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3090 : ALLOCATE (mixed_cdft%dlb_control%gradients(3*natom, &
3091 : mixed_cdft%dlb_control%bo(1):mixed_cdft%dlb_control%bo(2), &
3092 : mixed_cdft%dlb_control%bo(3):mixed_cdft%dlb_control%bo(4), &
3093 24 : mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3094 22724 : mixed_cdft%dlb_control%gradients = 0.0_dp
3095 3524 : mixed_cdft%dlb_control%weight = 0.0_dp
3096 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%cavity, source=mixed_cdft%dlb_control%my_source, &
3097 4 : request=req(1))
3098 4 : CALL req(1)%wait()
3099 4 : DEALLOCATE (mixed_cdft%dlb_control%bo)
3100 : ELSE
3101 0 : ALLOCATE (buffsize(1))
3102 : CALL force_env%para_env%irecv(msgout=buffsize, source=mixed_cdft%dlb_control%my_source, &
3103 0 : request=req(1))
3104 0 : CALL req(1)%wait()
3105 0 : ALLOCATE (mixed_cdft%dlb_control%bo(12*buffsize(1)))
3106 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%bo, source=mixed_cdft%dlb_control%my_source, &
3107 0 : request=req(1))
3108 0 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(buffsize(1)))
3109 0 : ALLOCATE (req_recv(buffsize(1)))
3110 0 : DEALLOCATE (buffsize)
3111 0 : CALL req(1)%wait()
3112 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
3113 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity(mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3114 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3115 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3116 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3117 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3118 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3119 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%sendbuff(j)%cavity, &
3120 : source=mixed_cdft%dlb_control%my_source, &
3121 0 : request=req_recv(j))
3122 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight(mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3123 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3124 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3125 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3126 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3127 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3128 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients(3*natom, &
3129 : mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3130 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3131 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3132 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3133 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3134 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3135 0 : mixed_cdft%dlb_control%sendbuff(j)%weight = 0.0_dp
3136 0 : mixed_cdft%dlb_control%sendbuff(j)%gradients = 0.0_dp
3137 : mixed_cdft%dlb_control%sendbuff(j)%tag = [mixed_cdft%dlb_control%bo(12*(j - 1) + 9), &
3138 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 10)]
3139 : mixed_cdft%dlb_control%sendbuff(j)%rank = [mixed_cdft%dlb_control%bo(12*(j - 1) + 11), &
3140 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 12)]
3141 : END DO
3142 0 : CALL mp_waitall(req_recv)
3143 0 : DEALLOCATE (req_recv)
3144 : END IF
3145 : END IF
3146 : ELSE
3147 4 : IF (.NOT. mixed_cdft%is_special) THEN
3148 4 : offset = 0
3149 4 : ALLOCATE (sendbuffer(12))
3150 4 : send_total = 0
3151 8 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3152 : tags = [(i - 1)*3 + 1 + force_env%para_env%mepos*6*max_targets, &
3153 12 : (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets] ! Unique communicator tags
3154 4 : mixed_cdft%dlb_control%target_list(3, i) = tags(1)
3155 4 : IF (mixed_cdft%is_pencil) THEN
3156 : sendbuffer = [bo_conf(1, 1) + offset, &
3157 : bo_conf(1, 1) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3158 : bo_conf(1, 2), bo_conf(2, 2), bo(1, 3), bo(2, 3), bo_conf(1, 3), bo_conf(2, 3), &
3159 0 : tags(1), tags(2), mixed_cdft%dest_list(1), mixed_cdft%dest_list(2)]
3160 : ELSE
3161 : sendbuffer = [bo_conf(1, 1), bo_conf(2, 1), &
3162 : bo_conf(1, 2) + offset, &
3163 : bo_conf(1, 2) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3164 : bo(1, 3), bo(2, 3), bo_conf(1, 3), bo_conf(2, 3), tags(1), tags(2), &
3165 52 : mixed_cdft%dest_list(1), mixed_cdft%dest_list(2)]
3166 : END IF
3167 4 : send_total = send_total + mixed_cdft%dlb_control%target_list(2, i) - 1
3168 : CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dlb_control%target_list(1, i), &
3169 4 : request=req(1))
3170 4 : CALL req(1)%wait()
3171 4 : IF (mixed_cdft%is_pencil) THEN
3172 : ALLOCATE (cavity(bo_conf(1, 1) + offset: &
3173 : bo_conf(1, 1) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3174 0 : bo_conf(1, 2):bo_conf(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
3175 : cavity = cdft_control%becke_control%cavity%array(bo_conf(1, 1) + offset: &
3176 : bo_conf(1, 1) + offset + &
3177 : (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3178 : bo_conf(1, 2):bo_conf(2, 2), &
3179 0 : bo_conf(1, 3):bo_conf(2, 3))
3180 : ELSE
3181 : ALLOCATE (cavity(bo_conf(1, 1):bo_conf(2, 1), &
3182 : bo_conf(1, 2) + offset: &
3183 : bo_conf(1, 2) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3184 20 : bo_conf(1, 3):bo_conf(2, 3)))
3185 : cavity = cdft_control%becke_control%cavity%array(bo_conf(1, 1):bo_conf(2, 1), &
3186 : bo_conf(1, 2) + offset: &
3187 : bo_conf(1, 2) + offset + &
3188 : (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3189 7044 : bo_conf(1, 3):bo_conf(2, 3))
3190 : END IF
3191 : CALL force_env%para_env%isend(msgin=cavity, &
3192 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3193 4 : request=req(1))
3194 4 : CALL req(1)%wait()
3195 4 : offset = offset + mixed_cdft%dlb_control%target_list(2, i)
3196 8 : DEALLOCATE (cavity)
3197 : END DO
3198 4 : IF (mixed_cdft%is_pencil) THEN
3199 0 : mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 1)
3200 0 : mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 1) + offset - 1
3201 : ELSE
3202 4 : mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 2)
3203 4 : mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 2) + offset - 1
3204 : END IF
3205 4 : DEALLOCATE (sendbuffer)
3206 : ELSE
3207 0 : ALLOCATE (buffsize(1))
3208 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3209 0 : buffsize = mixed_cdft%dlb_control%target_list(2, i)
3210 : ! Unique communicator tags (dont actually need these, should be removed)
3211 : tags = [(i - 1)*3 + 1 + force_env%para_env%mepos*6*max_targets, &
3212 0 : (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets]
3213 0 : DO j = 4, SIZE(mixed_cdft%dlb_control%target_list, 1)
3214 0 : IF (mixed_cdft%dlb_control%target_list(j, i) > uninitialized) EXIT
3215 : END DO
3216 0 : offset_special = j
3217 0 : offset_proc = j - 4 - (j - 4)/2
3218 : CALL force_env%para_env%isend(msgin=buffsize, &
3219 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3220 0 : request=req(1))
3221 0 : CALL req(1)%wait()
3222 0 : ALLOCATE (sendbuffer(12*buffsize(1)))
3223 0 : DO j = 1, buffsize(1)
3224 : sendbuffer(12*(j - 1) + 1:12*(j - 1) + 12) = [mixed_cdft%dlb_control%target_list(offset_special + 2*(j - 1), i), &
3225 : mixed_cdft%dlb_control%target_list(offset_special + 2*j - 1, i), &
3226 : bo_conf(1, 2), bo_conf(2, 2), bo(1, 3), bo(2, 3), &
3227 : bo_conf(1, 3), bo_conf(2, 3), tags(1), tags(2), &
3228 : mixed_cdft%dest_list(j + offset_proc), &
3229 0 : mixed_cdft%dest_list(j + offset_proc) + force_env%para_env%num_pe/2]
3230 : END DO
3231 : CALL force_env%para_env%isend(msgin=sendbuffer, &
3232 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3233 0 : request=req(1))
3234 0 : CALL req(1)%wait()
3235 0 : DEALLOCATE (sendbuffer)
3236 0 : DO j = 1, buffsize(1)
3237 : ALLOCATE (cavity(mixed_cdft%dlb_control%target_list(offset_special + 2*(j - 1), i): &
3238 : mixed_cdft%dlb_control%target_list(offset_special + 2*j - 1, i), &
3239 0 : bo_conf(1, 2):bo_conf(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
3240 : cavity = cdft_control%becke_control%cavity%array(LBOUND(cavity, 1):UBOUND(cavity, 1), &
3241 : bo_conf(1, 2):bo_conf(2, 2), &
3242 0 : bo_conf(1, 3):bo_conf(2, 3))
3243 : CALL force_env%para_env%isend(msgin=cavity, &
3244 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3245 0 : request=req(1))
3246 0 : CALL req(1)%wait()
3247 0 : DEALLOCATE (cavity)
3248 : END DO
3249 : END DO
3250 0 : DEALLOCATE (buffsize)
3251 : END IF
3252 : END IF
3253 8 : DEALLOCATE (expected_work, work_size, load_imbalance, work_index, targets)
3254 : ! Once calculated, data defined on the distributed grid points is sent directly to the processors that own the
3255 : ! grid points after the constraint is copied onto the two processor groups, instead of sending the data back to
3256 : ! the original owner
3257 8 : IF (mixed_cdft%is_special) THEN
3258 0 : my_special_work = 2
3259 0 : ALLOCATE (mask_send(SIZE(mixed_cdft%dest_list)), mask_recv(SIZE(mixed_cdft%source_list)))
3260 0 : ALLOCATE (nsend_proc(SIZE(mixed_cdft%dest_list)), nrecv(SIZE(mixed_cdft%source_list)))
3261 0 : nrecv = 0
3262 0 : nsend_proc = 0
3263 0 : mask_recv = .FALSE.
3264 0 : mask_send = .FALSE.
3265 : ELSE
3266 : my_special_work = 1
3267 : END IF
3268 72 : ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)), sbuff(SIZE(mixed_cdft%dest_list)))
3269 56 : ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%source_list) + (my_special_work**2)*SIZE(mixed_cdft%dest_list)))
3270 24 : ALLOCATE (mixed_cdft%dlb_control%recv_work_repl(SIZE(mixed_cdft%source_list)))
3271 24 : DO i = 1, SIZE(mixed_cdft%source_list)
3272 16 : NULLIFY (recvbuffer(i)%bv, recvbuffer(i)%iv)
3273 16 : ALLOCATE (recvbuffer(i)%bv(1), recvbuffer(i)%iv(3))
3274 : CALL force_env%para_env%irecv(msgout=recvbuffer(i)%bv, &
3275 : source=mixed_cdft%source_list(i), &
3276 16 : request=req_total(i), tag=1)
3277 24 : IF (mixed_cdft%is_special) THEN
3278 : CALL force_env%para_env%irecv(msgout=recvbuffer(i)%iv, &
3279 : source=mixed_cdft%source_list(i), &
3280 : request=req_total(i + SIZE(mixed_cdft%source_list)), &
3281 0 : tag=2)
3282 : END IF
3283 : END DO
3284 16 : DO i = 1, my_special_work
3285 32 : DO j = 1, SIZE(mixed_cdft%dest_list)
3286 16 : IF (i == 1) THEN
3287 16 : NULLIFY (sbuff(j)%iv, sbuff(j)%bv)
3288 16 : ALLOCATE (sbuff(j)%bv(1))
3289 32 : sbuff(j)%bv = mixed_cdft%dlb_control%send_work
3290 16 : IF (mixed_cdft%is_special) THEN
3291 0 : ALLOCATE (sbuff(j)%iv(3))
3292 0 : sbuff(j)%iv(1:2) = mixed_cdft%dest_list_bo(1:2, j)
3293 0 : sbuff(j)%iv(3) = 0
3294 0 : IF (sbuff(j)%iv(1) == should_deallocate) mask_send(j) = .TRUE.
3295 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3296 0 : sbuff(j)%bv = touched(j)
3297 0 : IF (touched(j)) THEN
3298 0 : nsend = 0
3299 0 : DO ispecial = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3300 0 : IF (mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), ispecial) /= uninitialized) THEN
3301 0 : nsend = nsend + 1
3302 : END IF
3303 : END DO
3304 0 : sbuff(j)%iv(3) = nsend
3305 0 : nsend_proc(j) = nsend
3306 : END IF
3307 : END IF
3308 : END IF
3309 : END IF
3310 16 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + my_special_work*SIZE(mixed_cdft%source_list)
3311 : CALL force_env%para_env%isend(msgin=sbuff(j)%bv, &
3312 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
3313 16 : request=req_total(ind), tag=1)
3314 24 : IF (mixed_cdft%is_special) THEN
3315 : CALL force_env%para_env%isend(msgin=sbuff(j)%iv, &
3316 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
3317 0 : request=req_total(ind + 2*SIZE(mixed_cdft%dest_list)), tag=2)
3318 : END IF
3319 : END DO
3320 : END DO
3321 8 : CALL mp_waitall(req_total)
3322 8 : DEALLOCATE (req_total)
3323 24 : DO i = 1, SIZE(mixed_cdft%source_list)
3324 16 : mixed_cdft%dlb_control%recv_work_repl(i) = recvbuffer(i)%bv(1)
3325 16 : IF (mixed_cdft%is_special .AND. mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3326 0 : mixed_cdft%source_list_bo(1:2, i) = recvbuffer(i)%iv(1:2)
3327 0 : nrecv(i) = recvbuffer(i)%iv(3)
3328 0 : IF (recvbuffer(i)%iv(1) == should_deallocate) mask_recv(i) = .TRUE.
3329 : END IF
3330 16 : DEALLOCATE (recvbuffer(i)%bv)
3331 24 : IF (ASSOCIATED(recvbuffer(i)%iv)) DEALLOCATE (recvbuffer(i)%iv)
3332 : END DO
3333 24 : DO j = 1, SIZE(mixed_cdft%dest_list)
3334 16 : DEALLOCATE (sbuff(j)%bv)
3335 24 : IF (ASSOCIATED(sbuff(j)%iv)) DEALLOCATE (sbuff(j)%iv)
3336 : END DO
3337 8 : DEALLOCATE (recvbuffer)
3338 : ! For some reason if debug_this_module is true and is_special is false, the deallocate statement
3339 : ! on line 3433 gets executed no matter what (gfortran 5.3.0 bug?). Printing out the variable seems to fix it...
3340 : IF (debug_this_module) THEN
3341 : WRITE (dummy, *) mixed_cdft%is_special
3342 : END IF
3343 :
3344 8 : IF (.NOT. mixed_cdft%is_special) THEN
3345 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3346 32 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl) + 2))
3347 4 : ALLOCATE (sendbuffer(6))
3348 4 : IF (mixed_cdft%is_pencil) THEN
3349 : sendbuffer = [SIZE(mixed_cdft%dlb_control%target_list, 2), bo_conf(1, 3), bo_conf(2, 3), &
3350 0 : bo_conf(1, 1), bo_conf(1, 2), bo_conf(2, 2)]
3351 : ELSE
3352 : sendbuffer = [SIZE(mixed_cdft%dlb_control%target_list, 2), bo_conf(1, 3), bo_conf(2, 3), &
3353 28 : bo_conf(1, 2), bo_conf(1, 1), bo_conf(2, 1)]
3354 : END IF
3355 8 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3356 24 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl)))
3357 : END IF
3358 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3359 24 : ALLOCATE (mixed_cdft%dlb_control%recv_info(2))
3360 8 : NULLIFY (mixed_cdft%dlb_control%recv_info(1)%target_list, mixed_cdft%dlb_control%recv_info(2)%target_list)
3361 24 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(2))
3362 8 : NULLIFY (mixed_cdft%dlb_control%recvbuff(1)%buffs, mixed_cdft%dlb_control%recvbuff(2)%buffs)
3363 : END IF
3364 : ! First communicate which grid points were distributed
3365 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3366 12 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl) + 1
3367 12 : DO i = 1, 2
3368 : CALL force_env%para_env%isend(msgin=sendbuffer, &
3369 : dest=mixed_cdft%dest_list(i), &
3370 8 : request=req_total(ind))
3371 12 : ind = ind + 1
3372 : END DO
3373 : END IF
3374 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3375 8 : ind = 1
3376 24 : DO i = 1, 2
3377 24 : IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3378 8 : ALLOCATE (mixed_cdft%dlb_control%recv_info(i)%matrix_info(6))
3379 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recv_info(i)%matrix_info, &
3380 : source=mixed_cdft%source_list(i), &
3381 8 : request=req_total(ind))
3382 8 : ind = ind + 1
3383 : END IF
3384 : END DO
3385 : END IF
3386 8 : IF (ASSOCIATED(req_total)) THEN
3387 8 : CALL mp_waitall(req_total)
3388 : END IF
3389 : ! Now communicate which processor handles which grid points
3390 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3391 12 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl) + 1
3392 12 : DO i = 1, 2
3393 8 : IF (i == 2) THEN
3394 8 : mixed_cdft%dlb_control%target_list(3, :) = mixed_cdft%dlb_control%target_list(3, :) + 3*max_targets
3395 : END IF
3396 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%target_list, &
3397 : dest=mixed_cdft%dest_list(i), &
3398 8 : request=req_total(ind))
3399 12 : ind = ind + 1
3400 : END DO
3401 : END IF
3402 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3403 8 : ind = 1
3404 24 : DO i = 1, 2
3405 24 : IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3406 : ALLOCATE (mixed_cdft%dlb_control%recv_info(i)% &
3407 24 : target_list(3, mixed_cdft%dlb_control%recv_info(i)%matrix_info(1)))
3408 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recv_info(i)%target_list, &
3409 : source=mixed_cdft%source_list(i), &
3410 8 : request=req_total(ind))
3411 8 : ind = ind + 1
3412 : END IF
3413 : END DO
3414 : END IF
3415 8 : IF (ASSOCIATED(req_total)) THEN
3416 8 : CALL mp_waitall(req_total)
3417 8 : DEALLOCATE (req_total)
3418 : END IF
3419 8 : IF (ASSOCIATED(sendbuffer)) DEALLOCATE (sendbuffer)
3420 : ELSE
3421 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3422 0 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl) + 2*COUNT(touched)))
3423 0 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3424 0 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl)))
3425 : END IF
3426 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3427 0 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl)
3428 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3429 0 : IF (touched(j)) THEN
3430 0 : ALLOCATE (sbuff(j)%iv(4 + 3*nsend_proc(j)))
3431 0 : sbuff(j)%iv(1:4) = [bo_conf(1, 2), bo_conf(2, 2), bo_conf(1, 3), bo_conf(2, 3)]
3432 0 : offset = 5
3433 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3434 0 : IF (mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), i) /= uninitialized) THEN
3435 : sbuff(j)%iv(offset:offset + 2) = [mixed_cdft%dlb_control%target_list(1, i), &
3436 : mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), i), &
3437 0 : mixed_cdft%dlb_control%target_list(4 + 2*j - 1, i)]
3438 0 : offset = offset + 3
3439 : END IF
3440 : END DO
3441 0 : DO ispecial = 1, my_special_work
3442 : CALL force_env%para_env%isend(msgin=sbuff(j)%iv, &
3443 : dest=mixed_cdft%dest_list(j) + (ispecial - 1)*force_env%para_env%num_pe/2, &
3444 0 : request=req_total(ind + ispecial))
3445 : END DO
3446 0 : ind = ind + my_special_work
3447 : END IF
3448 : END DO
3449 : END IF
3450 0 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3451 0 : ALLOCATE (mixed_cdft%dlb_control%recv_info(SIZE(mixed_cdft%source_list)))
3452 0 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(SIZE(mixed_cdft%source_list)))
3453 0 : ind = 1
3454 0 : DO j = 1, SIZE(mixed_cdft%source_list)
3455 : NULLIFY (mixed_cdft%dlb_control%recv_info(j)%target_list, &
3456 0 : mixed_cdft%dlb_control%recvbuff(j)%buffs)
3457 0 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3458 0 : ALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info(4 + 3*nrecv(j)))
3459 : CALL force_env%para_env%irecv(mixed_cdft%dlb_control%recv_info(j)%matrix_info, &
3460 : source=mixed_cdft%source_list(j), &
3461 0 : request=req_total(ind))
3462 0 : ind = ind + 1
3463 : END IF
3464 : END DO
3465 : END IF
3466 0 : IF (ASSOCIATED(req_total)) THEN
3467 0 : CALL mp_waitall(req_total)
3468 0 : DEALLOCATE (req_total)
3469 : END IF
3470 0 : IF (ANY(mask_send)) THEN
3471 : ALLOCATE (tmp(SIZE(mixed_cdft%dest_list) - COUNT(mask_send)), &
3472 0 : tmp_bo(2, SIZE(mixed_cdft%dest_list) - COUNT(mask_send)))
3473 0 : i = 1
3474 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3475 0 : IF (.NOT. mask_send(j)) THEN
3476 0 : tmp(i) = mixed_cdft%dest_list(j)
3477 0 : tmp_bo(1:2, i) = mixed_cdft%dest_list_bo(1:2, j)
3478 0 : i = i + 1
3479 : END IF
3480 : END DO
3481 0 : DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo)
3482 0 : ALLOCATE (mixed_cdft%dest_list(SIZE(tmp)), mixed_cdft%dest_list_bo(2, SIZE(tmp)))
3483 0 : mixed_cdft%dest_list = tmp
3484 0 : mixed_cdft%dest_list_bo = tmp_bo
3485 0 : DEALLOCATE (tmp, tmp_bo)
3486 : END IF
3487 0 : IF (ANY(mask_recv)) THEN
3488 : ALLOCATE (tmp(SIZE(mixed_cdft%source_list) - COUNT(mask_recv)), &
3489 0 : tmp_bo(4, SIZE(mixed_cdft%source_list) - COUNT(mask_recv)))
3490 0 : i = 1
3491 0 : DO j = 1, SIZE(mixed_cdft%source_list)
3492 0 : IF (.NOT. mask_recv(j)) THEN
3493 0 : tmp(i) = mixed_cdft%source_list(j)
3494 0 : tmp_bo(1:4, i) = mixed_cdft%source_list_bo(1:4, j)
3495 0 : i = i + 1
3496 : END IF
3497 : END DO
3498 0 : DEALLOCATE (mixed_cdft%source_list, mixed_cdft%source_list_bo)
3499 0 : ALLOCATE (mixed_cdft%source_list(SIZE(tmp)), mixed_cdft%source_list_bo(4, SIZE(tmp)))
3500 0 : mixed_cdft%source_list = tmp
3501 0 : mixed_cdft%source_list_bo = tmp_bo
3502 0 : DEALLOCATE (tmp, tmp_bo)
3503 : END IF
3504 0 : DEALLOCATE (mask_recv, mask_send)
3505 0 : DEALLOCATE (nsend_proc, nrecv)
3506 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3507 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3508 0 : IF (touched(j)) DEALLOCATE (sbuff(j)%iv)
3509 : END DO
3510 0 : IF (ASSOCIATED(touched)) DEALLOCATE (touched)
3511 : END IF
3512 : END IF
3513 8 : DEALLOCATE (sbuff)
3514 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
3515 8 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3516 8 : CALL timestop(handle)
3517 :
3518 16 : END SUBROUTINE mixed_becke_constraint_dlb
3519 :
3520 : ! **************************************************************************************************
3521 : !> \brief Low level routine to build mixed Becke constraint and gradients
3522 : !> \param force_env the force_env that holds the CDFT states
3523 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
3524 : !> \param in_memory decides whether to build the weight function gradients in parallel before solving
3525 : !> the CDFT states or later during the SCF procedure of the individual states
3526 : !> \param is_constraint a list used to determine which atoms in the system define the constraint
3527 : !> \param store_vectors should temporary arrays be stored in memory to accelerate the calculation
3528 : !> \param R12 temporary array holding the pairwise atomic distances
3529 : !> \param position_vecs temporary array holding the pbc corrected atomic position vectors
3530 : !> \param pair_dist_vecs temporary array holding the pairwise displament vectors
3531 : !> \param coefficients array that determines how atoms should be summed to form the constraint
3532 : !> \param catom temporary array to map the global index of constraint atoms to their position
3533 : !> in a list that holds only constraint atoms
3534 : !> \par History
3535 : !> 03.2016 created [Nico Holmberg]
3536 : ! **************************************************************************************************
3537 36 : SUBROUTINE mixed_becke_constraint_low(force_env, mixed_cdft, in_memory, &
3538 : is_constraint, store_vectors, R12, position_vecs, &
3539 : pair_dist_vecs, coefficients, catom)
3540 : TYPE(force_env_type), POINTER :: force_env
3541 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
3542 : LOGICAL, INTENT(IN) :: in_memory
3543 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: is_constraint
3544 : LOGICAL, INTENT(IN) :: store_vectors
3545 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3546 : INTENT(INOUT) :: R12, position_vecs
3547 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3548 : INTENT(INOUT) :: pair_dist_vecs
3549 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
3550 : INTENT(INOUT) :: coefficients
3551 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: catom
3552 :
3553 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_low'
3554 :
3555 : INTEGER :: handle, i, iatom, icomm, iforce_eval, index, iounit, ip, ispecial, iwork, j, &
3556 : jatom, jcomm, k, my_special_work, my_work, natom, nbuffs, nforce_eval, np(3), &
3557 : nsent_total, nskipped, nwork, offset, offset_repl
3558 36 : INTEGER, DIMENSION(:), POINTER :: work, work_dlb
3559 36 : INTEGER, DIMENSION(:, :), POINTER :: nsent
3560 : LOGICAL :: completed_recv, should_communicate
3561 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: skip_me
3562 36 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: completed
3563 : REAL(kind=dp) :: dist1, dist2, dmyexp, my1, my1_homo, &
3564 : myexp, sum_cell_f_all, &
3565 : sum_cell_f_constr, th, tmp_const
3566 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, distances, ds_dR_i, &
3567 36 : ds_dR_j
3568 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_const_dR, d_sum_Pm_dR, &
3569 36 : distance_vecs, dP_i_dRi
3570 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dP_i_dRj
3571 : REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dR_i, dmy_dR_j, &
3572 : dr, dr1_r2, dr_i_dR, dr_ij_dR, &
3573 : dr_j_dR, grid_p, r, r1, shift
3574 36 : REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
3575 36 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cavity, weight
3576 36 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: gradients
3577 : TYPE(cdft_control_type), POINTER :: cdft_control
3578 : TYPE(cell_type), POINTER :: cell
3579 : TYPE(cp_logger_type), POINTER :: logger
3580 : TYPE(cp_subsys_type), POINTER :: subsys_mix
3581 : TYPE(force_env_type), POINTER :: force_env_qs
3582 36 : TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
3583 36 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_send
3584 : TYPE(particle_list_type), POINTER :: particles
3585 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3586 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3587 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
3588 :
3589 72 : logger => cp_get_default_logger()
3590 36 : NULLIFY (work, req_recv, req_send, work_dlb, nsent, cutoffs, cavity, &
3591 36 : weight, gradients, cell, subsys_mix, force_env_qs, &
3592 36 : particle_set, particles, auxbas_pw_pool, force_env_section, &
3593 36 : print_section, cdft_control)
3594 36 : CALL timeset(routineN, handle)
3595 36 : nforce_eval = SIZE(force_env%sub_force_env)
3596 36 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
3597 36 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3598 36 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
3599 36 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
3600 : CALL force_env_get(force_env=force_env, &
3601 : subsys=subsys_mix, &
3602 28 : cell=cell)
3603 : CALL cp_subsys_get(subsys=subsys_mix, &
3604 : particles=particles, &
3605 28 : particle_set=particle_set)
3606 : ELSE
3607 24 : DO iforce_eval = 1, nforce_eval
3608 16 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
3609 24 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
3610 : END DO
3611 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
3612 : cp_subsys=subsys_mix, &
3613 8 : cell=cell)
3614 : CALL cp_subsys_get(subsys=subsys_mix, &
3615 : particles=particles, &
3616 8 : particle_set=particle_set)
3617 : END IF
3618 36 : natom = SIZE(particles%els)
3619 36 : cdft_control => mixed_cdft%cdft_control
3620 36 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
3621 144 : np = auxbas_pw_pool%pw_grid%npts
3622 144 : dr = auxbas_pw_pool%pw_grid%dr
3623 144 : shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp
3624 180 : ALLOCATE (cell_functions(natom), skip_me(natom))
3625 36 : IF (store_vectors) THEN
3626 72 : ALLOCATE (distances(natom))
3627 108 : ALLOCATE (distance_vecs(3, natom))
3628 : END IF
3629 36 : IF (in_memory) THEN
3630 24 : ALLOCATE (ds_dR_j(3))
3631 24 : ALLOCATE (ds_dR_i(3))
3632 72 : ALLOCATE (d_sum_Pm_dR(3, natom))
3633 48 : ALLOCATE (d_sum_const_dR(3, natom))
3634 96 : ALLOCATE (dP_i_dRj(3, natom, natom))
3635 48 : ALLOCATE (dP_i_dRi(3, natom))
3636 24 : th = 1.0e-8_dp
3637 : END IF
3638 36 : IF (mixed_cdft%dlb) THEN
3639 32 : ALLOCATE (work(force_env%para_env%num_pe), work_dlb(force_env%para_env%num_pe))
3640 24 : work = 0
3641 24 : work_dlb = 0
3642 : END IF
3643 36 : my_work = 1
3644 36 : my_special_work = 1
3645 : ! Load balancing: allocate storage for receiving buffers and post recv requests
3646 36 : IF (mixed_cdft%dlb) THEN
3647 8 : IF (mixed_cdft%dlb_control%recv_work) THEN
3648 4 : my_work = 2
3649 4 : IF (.NOT. mixed_cdft%is_special) THEN
3650 40 : ALLOCATE (req_send(2, 3))
3651 : ELSE
3652 0 : ALLOCATE (req_send(2, 3*SIZE(mixed_cdft%dlb_control%sendbuff)))
3653 : END IF
3654 : END IF
3655 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3656 8 : IF (.NOT. mixed_cdft%is_special) THEN
3657 8 : offset_repl = 0
3658 8 : IF (mixed_cdft%dlb_control%recv_work_repl(1) .AND. mixed_cdft%dlb_control%recv_work_repl(2)) THEN
3659 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2) + &
3660 0 : SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3661 0 : offset_repl = 3*SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2)
3662 8 : ELSE IF (mixed_cdft%dlb_control%recv_work_repl(1)) THEN
3663 0 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2))))
3664 : ELSE
3665 48 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3666 : END IF
3667 : ELSE
3668 0 : nbuffs = 0
3669 0 : offset_repl = 1
3670 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3671 0 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3672 0 : nbuffs = nbuffs + (SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3
3673 : END IF
3674 : END DO
3675 0 : ALLOCATE (req_recv(3*nbuffs))
3676 : END IF
3677 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3678 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3679 8 : IF (.NOT. mixed_cdft%is_special) THEN
3680 8 : offset = 0
3681 8 : index = j + (j/2)
3682 64 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)))
3683 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)
3684 8 : IF (mixed_cdft%is_pencil) THEN
3685 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3686 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3687 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3688 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3689 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3690 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3691 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3692 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3693 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3694 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3695 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3696 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3697 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3698 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3699 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3700 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3701 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3702 : gradients(3*natom, &
3703 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3704 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3705 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3706 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3707 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3708 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3709 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3710 : ELSE
3711 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3712 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3713 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3714 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3715 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3716 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3717 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3718 40 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3719 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3720 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3721 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3722 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3723 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3724 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3725 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3726 40 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3727 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3728 : gradients(3*natom, &
3729 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3730 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3731 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3732 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3733 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3734 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3735 48 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3736 : END IF
3737 :
3738 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, &
3739 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3740 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 1), &
3741 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i))
3742 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, &
3743 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3744 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 2), &
3745 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i) + 1)
3746 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, &
3747 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3748 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 3), &
3749 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i) + 2)
3750 16 : offset = offset + mixed_cdft%dlb_control%recv_info(j)%target_list(2, i)
3751 : END DO
3752 8 : DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info)
3753 : ELSE
3754 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)% &
3755 0 : buffs((SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3))
3756 0 : index = 6
3757 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
3758 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3759 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3760 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3761 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3762 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3763 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3764 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3765 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3766 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3767 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3768 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3769 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3770 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3771 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3772 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3773 : gradients(3*natom, mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3774 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3775 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3776 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3777 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3778 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3779 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, &
3780 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3781 0 : request=req_recv(offset_repl), tag=1)
3782 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, &
3783 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3784 0 : request=req_recv(offset_repl + 1), tag=2)
3785 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, &
3786 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3787 0 : request=req_recv(offset_repl + 2), tag=3)
3788 0 : index = index + 3
3789 0 : offset_repl = offset_repl + 3
3790 : END DO
3791 0 : DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info)
3792 : END IF
3793 : END IF
3794 : END DO
3795 : END IF
3796 : END IF
3797 36 : cutoffs => cdft_control%becke_control%cutoffs
3798 36 : should_communicate = .FALSE.
3799 144 : DO i = 1, 3
3800 144 : cell_v(i) = cell%hmat(i, i)
3801 : END DO
3802 76 : DO iwork = my_work, 1, -1
3803 40 : IF (iwork == 2) THEN
3804 4 : IF (.NOT. mixed_cdft%is_special) THEN
3805 4 : cavity => mixed_cdft%dlb_control%cavity
3806 4 : weight => mixed_cdft%dlb_control%weight
3807 4 : gradients => mixed_cdft%dlb_control%gradients
3808 4 : ALLOCATE (completed(2, 3), nsent(2, 3))
3809 : ELSE
3810 0 : my_special_work = SIZE(mixed_cdft%dlb_control%sendbuff)
3811 0 : ALLOCATE (completed(2, 3*my_special_work), nsent(2, 3*my_special_work))
3812 : END IF
3813 4 : completed = .FALSE.
3814 40 : nsent = 0
3815 : ELSE
3816 36 : IF (.NOT. mixed_cdft%is_special) THEN
3817 36 : weight => mixed_cdft%weight
3818 36 : cavity => mixed_cdft%cavity
3819 36 : gradients => cdft_control%group(1)%gradients
3820 : ELSE
3821 0 : my_special_work = SIZE(mixed_cdft%dest_list)
3822 : END IF
3823 : END IF
3824 116 : DO ispecial = 1, my_special_work
3825 40 : nwork = 0
3826 40 : IF (mixed_cdft%is_special) THEN
3827 0 : IF (iwork == 1) THEN
3828 0 : weight => mixed_cdft%sendbuff(ispecial)%weight
3829 0 : cavity => mixed_cdft%sendbuff(ispecial)%cavity
3830 0 : gradients => mixed_cdft%sendbuff(ispecial)%gradients
3831 : ELSE
3832 0 : weight => mixed_cdft%dlb_control%sendbuff(ispecial)%weight
3833 0 : cavity => mixed_cdft%dlb_control%sendbuff(ispecial)%cavity
3834 0 : gradients => mixed_cdft%dlb_control%sendbuff(ispecial)%gradients
3835 : END IF
3836 : END IF
3837 1021 : DO k = LBOUND(weight, 1), UBOUND(weight, 1)
3838 901 : IF (mixed_cdft%dlb .AND. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3839 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3840 0 : IF (k >= mixed_cdft%dlb_control%distributed(1) .AND. &
3841 : k <= mixed_cdft%dlb_control%distributed(2)) THEN
3842 : CYCLE
3843 : END IF
3844 : END IF
3845 : END IF
3846 41952 : DO j = LBOUND(weight, 2), UBOUND(weight, 2)
3847 39209 : IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3848 6400 : IF (mixed_cdft%dlb_control%send_work) THEN
3849 3120 : IF (j >= mixed_cdft%dlb_control%distributed(1) .AND. &
3850 : j <= mixed_cdft%dlb_control%distributed(2)) THEN
3851 : CYCLE
3852 : END IF
3853 : END IF
3854 : END IF
3855 : ! Check if any of the buffers have become available for deallocation
3856 39209 : IF (should_communicate) THEN
3857 108 : DO icomm = 1, SIZE(nsent, 2)
3858 270 : DO jcomm = 1, SIZE(nsent, 1)
3859 162 : IF (nsent(jcomm, icomm) == 1) CYCLE
3860 55 : completed(jcomm, icomm) = req_send(jcomm, icomm)%test()
3861 55 : IF (completed(jcomm, icomm)) THEN
3862 24 : nsent(jcomm, icomm) = nsent(jcomm, icomm) + 1
3863 24 : nsent_total = nsent_total + 1
3864 24 : IF (nsent_total == SIZE(nsent, 1)*SIZE(nsent, 2)) should_communicate = .FALSE.
3865 : END IF
3866 187 : IF (ALL(completed(:, icomm))) THEN
3867 12 : IF (MODULO(icomm, 3) == 1) THEN
3868 4 : IF (.NOT. mixed_cdft%is_special) THEN
3869 4 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
3870 : ELSE
3871 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%cavity)
3872 : END IF
3873 8 : ELSE IF (MODULO(icomm, 3) == 2) THEN
3874 4 : IF (.NOT. mixed_cdft%is_special) THEN
3875 4 : DEALLOCATE (mixed_cdft%dlb_control%weight)
3876 : ELSE
3877 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%weight)
3878 : END IF
3879 : ELSE
3880 4 : IF (.NOT. mixed_cdft%is_special) THEN
3881 4 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
3882 : ELSE
3883 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%gradients)
3884 : END IF
3885 : END IF
3886 : END IF
3887 : END DO
3888 : END DO
3889 : END IF
3890 : ! Poll to prevent starvation
3891 39209 : IF (ASSOCIATED(req_recv)) THEN
3892 6400 : completed_recv = mp_testall(req_recv)
3893 : END IF
3894 : !
3895 1926389 : DO i = LBOUND(weight, 3), UBOUND(weight, 3)
3896 1807861 : IF (cdft_control%becke_control%cavity_confine) THEN
3897 1716736 : IF (cavity(k, j, i) < cdft_control%becke_control%eps_cavity) CYCLE
3898 : END IF
3899 988401 : grid_p(1) = k*dr(1) + shift(1)
3900 988401 : grid_p(2) = j*dr(2) + shift(2)
3901 988401 : grid_p(3) = i*dr(3) + shift(3)
3902 988401 : nskipped = 0
3903 2965203 : cell_functions = 1.0_dp
3904 988401 : skip_me = .FALSE.
3905 988401 : IF (store_vectors) distances = 0.0_dp
3906 988401 : IF (in_memory) THEN
3907 675116 : d_sum_Pm_dR = 0.0_dp
3908 675116 : d_sum_const_dR = 0.0_dp
3909 675116 : dP_i_dRi = 0.0_dp
3910 : END IF
3911 2554568 : DO iatom = 1, natom
3912 1976802 : IF (skip_me(iatom)) THEN
3913 65541 : cell_functions(iatom) = 0.0_dp
3914 65541 : IF (cdft_control%becke_control%should_skip) THEN
3915 37642 : IF (is_constraint(iatom)) nskipped = nskipped + 1
3916 37642 : IF (nskipped == cdft_control%natoms) THEN
3917 0 : IF (in_memory) THEN
3918 0 : IF (cdft_control%becke_control%cavity_confine) THEN
3919 0 : cavity(k, j, i) = 0.0_dp
3920 : END IF
3921 : END IF
3922 : EXIT
3923 : END IF
3924 : END IF
3925 : CYCLE
3926 : END IF
3927 1911261 : IF (store_vectors) THEN
3928 1911261 : IF (distances(iatom) == 0.0_dp) THEN
3929 6889176 : r = position_vecs(:, iatom)
3930 6889176 : dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
3931 6889176 : dist1 = NORM2(dist_vec)
3932 6889176 : distance_vecs(:, iatom) = dist_vec
3933 1722294 : distances(iatom) = dist1
3934 : ELSE
3935 755868 : dist_vec = distance_vecs(:, iatom)
3936 : dist1 = distances(iatom)
3937 : END IF
3938 : ELSE
3939 0 : r = particle_set(iatom)%r
3940 0 : DO ip = 1, 3
3941 0 : r(ip) = MODULO(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3942 : END DO
3943 0 : dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
3944 0 : dist1 = NORM2(dist_vec)
3945 : END IF
3946 2489027 : IF (dist1 <= cutoffs(iatom)) THEN
3947 509029 : IF (in_memory) THEN
3948 : IF (dist1 <= th) dist1 = th
3949 1184700 : dr_i_dR(:) = dist_vec(:)/dist1
3950 : END IF
3951 1527087 : DO jatom = 1, natom
3952 1527087 : IF (jatom /= iatom) THEN
3953 509029 : IF (jatom < iatom) THEN
3954 254521 : IF (.NOT. skip_me(jatom)) CYCLE
3955 : END IF
3956 320062 : IF (store_vectors) THEN
3957 320062 : IF (distances(jatom) == 0.0_dp) THEN
3958 1018032 : r1 = position_vecs(:, jatom)
3959 1018032 : dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
3960 1018032 : dist2 = NORM2(dist_vec)
3961 1018032 : distance_vecs(:, jatom) = dist_vec
3962 254508 : distances(jatom) = dist2
3963 : ELSE
3964 262216 : dist_vec = distance_vecs(:, jatom)
3965 : dist2 = distances(jatom)
3966 : END IF
3967 : ELSE
3968 0 : r1 = particle_set(jatom)%r
3969 0 : DO ip = 1, 3
3970 0 : r1(ip) = MODULO(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3971 : END DO
3972 0 : dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
3973 0 : dist2 = NORM2(dist_vec)
3974 : END IF
3975 320062 : IF (in_memory) THEN
3976 191129 : IF (store_vectors) THEN
3977 764516 : dr1_r2 = pair_dist_vecs(:, iatom, jatom)
3978 : ELSE
3979 0 : dr1_r2 = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
3980 : END IF
3981 : IF (dist2 <= th) dist2 = th
3982 191129 : tmp_const = (R12(iatom, jatom)**3)
3983 764516 : dr_ij_dR(:) = dr1_r2(:)/tmp_const
3984 : !derivativ w.r.t. Rj
3985 764516 : dr_j_dR = dist_vec(:)/dist2
3986 764516 : dmy_dR_j(:) = -(dr_j_dR(:)/R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:))
3987 : !derivativ w.r.t. Ri
3988 764516 : dmy_dR_i(:) = dr_i_dR(:)/R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:)
3989 : END IF
3990 320062 : my1 = (dist1 - dist2)/R12(iatom, jatom)
3991 320062 : IF (cdft_control%becke_control%adjust) THEN
3992 166839 : my1_homo = my1
3993 : my1 = my1 + &
3994 166839 : cdft_control%becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
3995 : END IF
3996 320062 : myexp = 1.5_dp*my1 - 0.5_dp*my1**3
3997 320062 : IF (in_memory) THEN
3998 191129 : dmyexp = 1.5_dp - 1.5_dp*my1**2
3999 : tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
4000 191129 : (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
4001 :
4002 764516 : ds_dR_i(:) = -0.5_dp*tmp_const*dmy_dR_i(:)
4003 764516 : ds_dR_j(:) = -0.5_dp*tmp_const*dmy_dR_j(:)
4004 191129 : IF (cdft_control%becke_control%adjust) THEN
4005 102514 : tmp_const = 1.0_dp - 2.0_dp*my1_homo*cdft_control%becke_control%aij(iatom, jatom)
4006 410056 : ds_dR_i(:) = ds_dR_i(:)*tmp_const
4007 410056 : ds_dR_j(:) = ds_dR_j(:)*tmp_const
4008 : END IF
4009 : END IF
4010 320062 : myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4011 320062 : myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4012 320062 : tmp_const = 0.5_dp*(1.0_dp - myexp)
4013 320062 : cell_functions(iatom) = cell_functions(iatom)*tmp_const
4014 320062 : IF (in_memory) THEN
4015 191129 : IF (ABS(tmp_const) <= th) tmp_const = tmp_const + th
4016 764516 : dP_i_dRi(:, iatom) = dP_i_dRi(:, iatom) + ds_dR_i(:)/tmp_const
4017 764516 : dP_i_dRj(:, iatom, jatom) = ds_dR_j(:)/tmp_const
4018 : END IF
4019 :
4020 320062 : IF (dist2 <= cutoffs(jatom)) THEN
4021 188967 : tmp_const = 0.5_dp*(1.0_dp + myexp)
4022 188967 : cell_functions(jatom) = cell_functions(jatom)*tmp_const
4023 188967 : IF (in_memory) THEN
4024 105046 : IF (ABS(tmp_const) <= th) tmp_const = tmp_const + th
4025 420184 : dP_i_dRj(:, jatom, iatom) = -ds_dR_i(:)/tmp_const
4026 420184 : dP_i_dRi(:, jatom) = dP_i_dRi(:, jatom) - ds_dR_j(:)/tmp_const
4027 : END IF
4028 : ELSE
4029 131095 : skip_me(jatom) = .TRUE.
4030 : END IF
4031 : END IF
4032 : END DO
4033 509029 : IF (in_memory) THEN
4034 1184700 : dP_i_dRi(:, iatom) = cell_functions(iatom)*dP_i_dRi(:, iatom)
4035 1184700 : d_sum_Pm_dR(:, iatom) = d_sum_Pm_dR(:, iatom) + dP_i_dRi(:, iatom)
4036 296175 : IF (is_constraint(iatom)) THEN
4037 : d_sum_const_dR(:, iatom) = d_sum_const_dR(:, iatom) + dP_i_dRi(:, iatom)* &
4038 1184700 : coefficients(iatom)
4039 : END IF
4040 888525 : DO jatom = 1, natom
4041 888525 : IF (jatom /= iatom) THEN
4042 296175 : IF (jatom < iatom) THEN
4043 148094 : IF (.NOT. skip_me(jatom)) THEN
4044 420184 : dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom)
4045 420184 : d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)
4046 105046 : IF (is_constraint(iatom)) THEN
4047 : d_sum_const_dR(:, jatom) = d_sum_const_dR(:, jatom) + &
4048 : dP_i_dRj(:, iatom, jatom)* &
4049 420184 : coefficients(iatom)
4050 : END IF
4051 : CYCLE
4052 : END IF
4053 : END IF
4054 764516 : dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom)
4055 764516 : d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)
4056 191129 : IF (is_constraint(iatom)) THEN
4057 : d_sum_const_dR(:, jatom) = d_sum_const_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)* &
4058 764516 : coefficients(iatom)
4059 : END IF
4060 : END IF
4061 : END DO
4062 : END IF
4063 : ELSE
4064 1402232 : cell_functions(iatom) = 0.0_dp
4065 1402232 : skip_me(iatom) = .TRUE.
4066 1402232 : IF (cdft_control%becke_control%should_skip) THEN
4067 858928 : IF (is_constraint(iatom)) nskipped = nskipped + 1
4068 858928 : IF (nskipped == cdft_control%natoms) THEN
4069 410635 : IF (in_memory) THEN
4070 252800 : IF (cdft_control%becke_control%cavity_confine) THEN
4071 252800 : cavity(k, j, i) = 0.0_dp
4072 : END IF
4073 : END IF
4074 : EXIT
4075 : END IF
4076 : END IF
4077 : END IF
4078 : END DO
4079 988401 : IF (nskipped == cdft_control%natoms) CYCLE
4080 : sum_cell_f_constr = 0.0_dp
4081 1733298 : DO ip = 1, cdft_control%natoms
4082 : sum_cell_f_constr = sum_cell_f_constr + cell_functions(catom(ip))* &
4083 1733298 : cdft_control%group(1)%coeff(ip)
4084 : END DO
4085 577766 : sum_cell_f_all = 0.0_dp
4086 577766 : nwork = nwork + 1
4087 1733298 : DO ip = 1, natom
4088 1733298 : sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
4089 : END DO
4090 577766 : IF (in_memory) THEN
4091 1266948 : DO iatom = 1, natom
4092 1266948 : IF (ABS(sum_cell_f_all) > 0.0_dp) THEN
4093 : gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
4094 : d_sum_const_dR(:, iatom)/sum_cell_f_all - sum_cell_f_constr* &
4095 1529032 : d_sum_Pm_dR(:, iatom)/(sum_cell_f_all**2)
4096 : END IF
4097 : END DO
4098 : END IF
4099 616975 : IF (ABS(sum_cell_f_all) > 0.000001) THEN
4100 320062 : weight(k, j, i) = sum_cell_f_constr/sum_cell_f_all
4101 : END IF
4102 : END DO ! i
4103 : END DO ! j
4104 : END DO ! k
4105 : ! Load balancing: post send requests
4106 80 : IF (iwork == 2) THEN
4107 4 : IF (.NOT. mixed_cdft%is_special) THEN
4108 12 : DO i = 1, SIZE(req_send, 1)
4109 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%cavity, &
4110 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4111 : request=req_send(i, 1), &
4112 8 : tag=mixed_cdft%dlb_control%dest_tags_repl(i))
4113 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%weight, &
4114 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4115 : request=req_send(i, 2), &
4116 8 : tag=mixed_cdft%dlb_control%dest_tags_repl(i) + 1)
4117 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%gradients, &
4118 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4119 : request=req_send(i, 3), &
4120 12 : tag=mixed_cdft%dlb_control%dest_tags_repl(i) + 2)
4121 : END DO
4122 : should_communicate = .TRUE.
4123 : nsent_total = 0
4124 : ELSE
4125 0 : DO i = 1, SIZE(req_send, 1)
4126 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%cavity, &
4127 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4128 0 : request=req_send(i, 3*(ispecial - 1) + 1), tag=1)
4129 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%weight, &
4130 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4131 0 : request=req_send(i, 3*(ispecial - 1) + 2), tag=2)
4132 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%gradients, &
4133 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4134 0 : request=req_send(i, 3*(ispecial - 1) + 3), tag=3)
4135 : END DO
4136 0 : IF (ispecial == my_special_work) THEN
4137 0 : should_communicate = .TRUE.
4138 0 : nsent_total = 0
4139 : END IF
4140 : END IF
4141 4 : work(mixed_cdft%dlb_control%my_source + 1) = work(mixed_cdft%dlb_control%my_source + 1) + nwork
4142 4 : work_dlb(force_env%para_env%mepos + 1) = work_dlb(force_env%para_env%mepos + 1) + nwork
4143 : ELSE
4144 36 : IF (mixed_cdft%dlb) work(force_env%para_env%mepos + 1) = work(force_env%para_env%mepos + 1) + nwork
4145 36 : IF (mixed_cdft%dlb) work_dlb(force_env%para_env%mepos + 1) = work_dlb(force_env%para_env%mepos + 1) + nwork
4146 : END IF
4147 : END DO ! ispecial
4148 : END DO ! iwork
4149 : ! Load balancing: wait for communication and deallocate sending buffers
4150 36 : IF (mixed_cdft%dlb) THEN
4151 16 : IF (mixed_cdft%dlb_control%recv_work .AND. &
4152 : ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
4153 48 : ALLOCATE (req_total(SIZE(req_recv) + SIZE(req_send, 1)*SIZE(req_send, 2)))
4154 4 : index = SIZE(req_recv)
4155 28 : req_total(1:index) = req_recv
4156 16 : DO i = 1, SIZE(req_send, 2)
4157 40 : DO j = 1, SIZE(req_send, 1)
4158 24 : index = index + 1
4159 36 : req_total(index) = req_send(j, i)
4160 : END DO
4161 : END DO
4162 4 : CALL mp_waitall(req_total)
4163 4 : DEALLOCATE (req_total)
4164 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4165 0 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
4166 : END IF
4167 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4168 0 : DEALLOCATE (mixed_cdft%dlb_control%weight)
4169 : END IF
4170 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4171 0 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
4172 : END IF
4173 4 : IF (mixed_cdft%is_special) THEN
4174 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4175 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4176 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4177 : END IF
4178 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4179 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4180 : END IF
4181 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4182 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4183 : END IF
4184 : END DO
4185 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4186 : END IF
4187 4 : DEALLOCATE (req_send, req_recv)
4188 4 : ELSE IF (mixed_cdft%dlb_control%recv_work) THEN
4189 0 : IF (should_communicate) THEN
4190 0 : CALL mp_waitall(req_send)
4191 : END IF
4192 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4193 0 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
4194 : END IF
4195 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4196 0 : DEALLOCATE (mixed_cdft%dlb_control%weight)
4197 : END IF
4198 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4199 0 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
4200 : END IF
4201 0 : IF (mixed_cdft%is_special) THEN
4202 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4203 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4204 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4205 : END IF
4206 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4207 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4208 : END IF
4209 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4210 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4211 : END IF
4212 : END DO
4213 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4214 : END IF
4215 0 : DEALLOCATE (req_send)
4216 8 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
4217 4 : CALL mp_waitall(req_recv)
4218 4 : DEALLOCATE (req_recv)
4219 : END IF
4220 : END IF
4221 36 : IF (mixed_cdft%dlb) THEN
4222 40 : CALL force_env%para_env%sum(work)
4223 40 : CALL force_env%para_env%sum(work_dlb)
4224 8 : IF (.NOT. ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
4225 12 : ALLOCATE (mixed_cdft%dlb_control%prediction_error(force_env%para_env%num_pe))
4226 : END IF
4227 40 : mixed_cdft%dlb_control%prediction_error = mixed_cdft%dlb_control%expected_work - work
4228 : IF (debug_this_module .AND. iounit > 0) THEN
4229 : DO i = 1, SIZE(work, 1)
4230 : WRITE (iounit, '(A,I10,I10,I10)') &
4231 : 'Work', work(i), work_dlb(i), mixed_cdft%dlb_control%expected_work(i)
4232 : END DO
4233 : END IF
4234 8 : DEALLOCATE (work, work_dlb, mixed_cdft%dlb_control%expected_work)
4235 : END IF
4236 36 : NULLIFY (gradients, weight, cavity)
4237 36 : IF (ALLOCATED(coefficients)) THEN
4238 36 : DEALLOCATE (coefficients)
4239 : END IF
4240 36 : IF (in_memory) THEN
4241 24 : DEALLOCATE (ds_dR_j)
4242 24 : DEALLOCATE (ds_dR_i)
4243 24 : DEALLOCATE (d_sum_Pm_dR)
4244 24 : DEALLOCATE (d_sum_const_dR)
4245 24 : DEALLOCATE (dP_i_dRj)
4246 24 : DEALLOCATE (dP_i_dRi)
4247 24 : NULLIFY (gradients)
4248 24 : IF (store_vectors) THEN
4249 24 : DEALLOCATE (pair_dist_vecs)
4250 : END IF
4251 : END IF
4252 36 : NULLIFY (cutoffs)
4253 36 : IF (ALLOCATED(is_constraint)) THEN
4254 36 : DEALLOCATE (is_constraint)
4255 : END IF
4256 36 : DEALLOCATE (catom)
4257 36 : DEALLOCATE (R12)
4258 36 : DEALLOCATE (cell_functions)
4259 36 : DEALLOCATE (skip_me)
4260 36 : IF (ALLOCATED(completed)) THEN
4261 4 : DEALLOCATE (completed)
4262 : END IF
4263 36 : IF (ASSOCIATED(nsent)) THEN
4264 4 : DEALLOCATE (nsent)
4265 : END IF
4266 36 : IF (store_vectors) THEN
4267 36 : DEALLOCATE (distances)
4268 36 : DEALLOCATE (distance_vecs)
4269 36 : DEALLOCATE (position_vecs)
4270 : END IF
4271 36 : IF (ASSOCIATED(req_send)) THEN
4272 0 : DEALLOCATE (req_send)
4273 : END IF
4274 36 : IF (ASSOCIATED(req_recv)) THEN
4275 0 : DEALLOCATE (req_recv)
4276 : END IF
4277 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
4278 36 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
4279 36 : CALL timestop(handle)
4280 :
4281 72 : END SUBROUTINE mixed_becke_constraint_low
4282 :
4283 0 : END MODULE mixed_cdft_methods
|
