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 532 : 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 532 : NULLIFY (subsys_mix, force_env_qs, force_env_section, print_section, &
174 532 : root_section, mixed_section, md_section, mixed_env, mixed_cdft, &
175 532 : 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 532 : is_qmmm = .FALSE.
183 1064 : logger => cp_get_default_logger()
184 532 : CPASSERT(ASSOCIATED(force_env))
185 532 : nforce_eval = SIZE(force_env%sub_force_env)
186 532 : CALL timeset(routineN, handle)
187 532 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
188 532 : mixed_env => force_env%mixed_env
189 532 : mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
190 532 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
191 532 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
192 : ! Check if a mixed CDFT calculation is requested
193 532 : CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
194 532 : IF (mixing_type == mix_cdft .AND. .NOT. ASSOCIATED(mixed_env%cdft_control)) THEN
195 80 : mixed_env%do_mixed_cdft = .TRUE.
196 160 : IF (mixed_env%do_mixed_cdft) THEN
197 : ! Sanity check
198 80 : 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 80 : mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
203 80 : CALL section_vals_get(mapping_section, explicit=explicit)
204 : ! The sub_force_envs must share the same geometrical structure
205 80 : IF (explicit) THEN
206 0 : CPABORT("Please disable section &MAPPING for mixed CDFT calculations")
207 : END IF
208 80 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%COUPLING", i_val=et_freq)
209 80 : IF (et_freq < 0) THEN
210 0 : mixed_env%do_mixed_et = .FALSE.
211 : ELSE
212 80 : mixed_env%do_mixed_et = .TRUE.
213 80 : IF (et_freq == 0) THEN
214 0 : mixed_env%et_freq = 1
215 : ELSE
216 80 : 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 264 : DO iforce_eval = 1, nforce_eval
222 184 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
223 304 : SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
224 : CASE (use_qs_force)
225 146 : 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 158 : "supported for mixed_cdft_init")
236 : END SELECT
237 238 : CPASSERT(ASSOCIATED(force_env_qs))
238 : END DO
239 : ! Get infos about the mixed subsys
240 80 : IF (.NOT. is_qmmm) THEN
241 : CALL force_env_get(force_env=force_env, &
242 72 : subsys=subsys_mix)
243 : CALL cp_subsys_get(subsys=subsys_mix, &
244 72 : 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 80 : ALLOCATE (mixed_cdft)
253 80 : CALL mixed_cdft_type_create(mixed_cdft)
254 80 : mixed_cdft%first_iteration = .TRUE.
255 : ! Determine what run type to use
256 80 : IF (mixed_env%ngroups == 1) THEN
257 : ! States treated in serial, possibly copying CDFT weight function and gradients from state to state
258 54 : 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 80 : mixed_env%do_mixed_qmmm_cdft = is_qmmm
277 : ! Setup dynamic load balancing
278 80 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%DLB", l_val=mixed_cdft%dlb)
279 80 : mixed_cdft%dlb = mixed_cdft%dlb .AND. calculate_forces ! disable if forces are not needed
280 80 : mixed_cdft%dlb = mixed_cdft%dlb .AND. (mixed_cdft%run_type == mixed_cdft_parallel) ! disable if not parallel
281 80 : 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 80 : mixed_cdft%calculate_metric = .FALSE.
298 80 : mixed_cdft%wfn_overlap_method = .FALSE.
299 80 : mixed_cdft%use_lowdin = .FALSE.
300 80 : mixed_cdft%do_ci = .FALSE.
301 80 : mixed_cdft%nonortho_coupling = .FALSE.
302 80 : mixed_cdft%identical_constraints = .TRUE.
303 80 : mixed_cdft%block_diagonalize = .FALSE.
304 80 : IF (mixed_env%do_mixed_et) THEN
305 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%METRIC", &
306 80 : l_val=mixed_cdft%calculate_metric)
307 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%WFN_OVERLAP", &
308 80 : l_val=mixed_cdft%wfn_overlap_method)
309 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOWDIN", &
310 80 : l_val=mixed_cdft%use_lowdin)
311 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%CI", &
312 80 : l_val=mixed_cdft%do_ci)
313 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%NONORTHOGONAL_COUPLING", &
314 80 : l_val=mixed_cdft%nonortho_coupling)
315 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%BLOCK_DIAGONALIZE", &
316 80 : l_val=mixed_cdft%block_diagonalize)
317 : END IF
318 : ! Inversion method
319 80 : CALL section_vals_val_get(mixed_section, "MIXED_CDFT%EPS_SVD", r_val=mixed_cdft%eps_svd)
320 80 : 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 80 : CALL force_env_get(force_env, root_section=root_section)
325 80 : md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
326 80 : CALL section_vals_val_get(md_section, "TIMESTEP", r_val=mixed_cdft%sim_dt)
327 80 : CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=mixed_cdft%sim_step)
328 80 : mixed_cdft%sim_step = mixed_cdft%sim_step - 1 ! to get the first step correct
329 80 : 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 80 : settings, natom=SIZE(particles_mix%els))
333 : ! Transfer settings to mixed_cdft
334 80 : CALL mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
335 : ! Initilize necessary structures
336 80 : CALL mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
337 : ! Write information about the mixed CDFT calculation
338 80 : IF (iounit > 0) THEN
339 40 : WRITE (iounit, *) ""
340 : WRITE (iounit, FMT="(T2,A,T71)") &
341 40 : "MIXED_CDFT| Activating mixed CDFT calculation"
342 : WRITE (iounit, FMT="(T2,A,T71,I10)") &
343 40 : "MIXED_CDFT| Number of CDFT states: ", nforce_eval
344 52 : 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 27 : "MIXED_CDFT| CDFT states calculation mode: serial"
355 27 : IF (mixed_cdft%identical_constraints) THEN
356 : WRITE (iounit, FMT="(T2,A,T71)") &
357 26 : "MIXED_CDFT| The constraints are built before the SCF procedure of the first"
358 : WRITE (iounit, FMT="(T2,A,T71)") &
359 26 : " 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 40 : CPABORT("Unknown mixed CDFT run type.")
371 : END SELECT
372 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
373 40 : "MIXED_CDFT| Calculating electronic coupling between states: ", mixed_env%do_mixed_et
374 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
375 40 : "MIXED_CDFT| Calculating electronic coupling reliability metric: ", mixed_cdft%calculate_metric
376 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
377 40 : "MIXED_CDFT| Configuration interaction (CDFT-CI) was requested: ", mixed_cdft%do_ci
378 : WRITE (iounit, FMT="(T2,A,T71,L10)") &
379 40 : "MIXED_CDFT| Block diagonalizing the mixed CDFT Hamiltonian: ", mixed_cdft%block_diagonalize
380 40 : 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 40 : IF (mixed_env%do_mixed_et) THEN
394 40 : IF (mixed_cdft%eps_svd == 0.0_dp) THEN
395 32 : 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 80 : 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 532 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
407 532 : CALL timestop(handle)
408 :
409 5320 : 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 292 : 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 292 : NULLIFY (mixed_cdft)
427 292 : CPASSERT(ASSOCIATED(force_env))
428 292 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
429 292 : CPASSERT(ASSOCIATED(mixed_cdft))
430 292 : IF (.NOT. PRESENT(iforce_eval)) THEN
431 138 : 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 102 : CALL mixed_cdft_set_flags(force_env)
436 : CASE DEFAULT
437 : ! Do nothing
438 : END SELECT
439 : ELSE
440 342 : SELECT CASE (mixed_cdft%run_type)
441 : CASE (mixed_cdft_serial)
442 190 : 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 292 : 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 304 : 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 152 : NULLIFY (mixed_cdft, dft_control_source, dft_control_target, force_env_qs_source, &
905 152 : force_env_qs_target, pw_env_source, pw_env_target, auxbas_pw_pool_source, &
906 152 : auxbas_pw_pool_target, qs_env_source, qs_env_target, mixed_env, &
907 152 : cdft_control_source, cdft_control_target)
908 152 : mixed_env => force_env%mixed_env
909 152 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
910 152 : CALL timeset(routineN, handle)
911 152 : IF (iforce_eval == 1) THEN
912 : jforce_eval = 1
913 : ELSE
914 88 : jforce_eval = iforce_eval - 1
915 : END IF
916 152 : nforce_eval = SIZE(force_env%sub_force_env)
917 304 : SELECT CASE (force_env%sub_force_env(jforce_eval)%force_env%in_use)
918 : CASE (use_qs_force, use_qmmm)
919 152 : force_env_qs_source => force_env%sub_force_env(jforce_eval)%force_env
920 152 : 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 152 : "supported for mixed_cdft_transfer_weight")
925 : END SELECT
926 152 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
927 : CALL force_env_get(force_env=force_env_qs_source, &
928 136 : qs_env=qs_env_source)
929 : CALL force_env_get(force_env=force_env_qs_target, &
930 136 : 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 152 : 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 64 : CALL get_qs_env(qs_env_source, dft_control=dft_control_source)
939 64 : cdft_control_source => dft_control_source%qs_control%cdft_control
940 64 : cdft_control_source%external_control = .FALSE.
941 64 : cdft_control_source%need_pot = .TRUE.
942 64 : IF (mixed_cdft%identical_constraints) THEN
943 62 : cdft_control_source%transfer_pot = .TRUE.
944 : ELSE
945 2 : cdft_control_source%transfer_pot = .FALSE.
946 : END IF
947 64 : 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 88 : pw_env=pw_env_source)
952 88 : CALL pw_env_get(pw_env_source, auxbas_pw_pool=auxbas_pw_pool_source)
953 88 : cdft_control_source => dft_control_source%qs_control%cdft_control
954 : CALL get_qs_env(qs_env_target, dft_control=dft_control_target, &
955 88 : pw_env=pw_env_target)
956 88 : CALL pw_env_get(pw_env_target, auxbas_pw_pool=auxbas_pw_pool_target)
957 88 : 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 88 : IF (mixed_cdft%identical_constraints) THEN
960 : ! Weight function
961 174 : DO igroup = 1, SIZE(cdft_control_target%group)
962 88 : ALLOCATE (cdft_control_target%group(igroup)%weight)
963 88 : 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 88 : CALL pw_copy(cdft_control_source%group(igroup)%weight, cdft_control_target%group(igroup)%weight)
966 88 : CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%group(igroup)%weight)
967 174 : DEALLOCATE (cdft_control_source%group(igroup)%weight)
968 : END DO
969 : ! Cavity
970 86 : IF (cdft_control_source%type == outer_scf_becke_constraint) THEN
971 80 : 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 86 : 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 86 : 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 86 : cdft_control_target%external_control = .FALSE.
1008 86 : cdft_control_target%need_pot = .FALSE.
1009 : ! For states i+1 < nforce_eval, prevent deallocation of constraint
1010 86 : IF (iforce_eval == nforce_eval) THEN
1011 62 : cdft_control_target%transfer_pot = .FALSE.
1012 : ELSE
1013 24 : cdft_control_target%transfer_pot = .TRUE.
1014 : END IF
1015 86 : 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 152 : IF (mixed_env%do_mixed_et) THEN
1030 152 : IF (MODULO(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
1031 152 : IF (iforce_eval == 1) THEN
1032 64 : cdft_control_source%do_et = .TRUE.
1033 64 : cdft_control_source%calculate_metric = mixed_cdft%calculate_metric
1034 : ELSE
1035 88 : cdft_control_target%do_et = .TRUE.
1036 88 : 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 152 : IF (iforce_eval == nforce_eval .AND. first_call) first_call = .FALSE.
1049 152 : CALL timestop(handle)
1050 :
1051 152 : 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 204 : 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 102 : CPASSERT(ASSOCIATED(force_env))
1137 102 : CALL timeset(routineN, handle)
1138 : ! Move needed arrays from individual CDFT states to the mixed CDFT env
1139 102 : 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 102 : CALL mixed_cdft_interaction_matrices(force_env)
1144 : ! Calculate eletronic couplings between states (Lowdin/rotation)
1145 102 : CALL mixed_cdft_calculate_coupling_low(force_env)
1146 : ! Print out couplings
1147 102 : CALL mixed_cdft_print_couplings(force_env)
1148 : ! Block diagonalize the mixed CDFT Hamiltonian matrix
1149 102 : CALL mixed_cdft_block_diag(force_env)
1150 : ! CDFT Configuration Interaction
1151 102 : CALL mixed_cdft_configuration_interaction(force_env)
1152 : ! Clean up
1153 102 : CALL mixed_cdft_release_work(force_env)
1154 102 : CALL timestop(handle)
1155 :
1156 102 : 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 102 : 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 102 : 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 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_nonortho, eigenv, energy, Sda
1180 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, S_det, S_mat, strength, tmp_mat, &
1181 102 : W_diagonal, Wad, Wda
1182 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: a, b
1183 102 : 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 102 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_overlap
1187 102 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: w_matrix_mo
1188 102 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
1189 : TYPE(cp_logger_type), POINTER :: logger
1190 102 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix, density_matrix_diff, &
1191 102 : 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 102 : NULLIFY (force_env_section, print_section, mixed_cdft_section, &
1202 102 : mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
1203 102 : density_matrix_diff, mo_mo_fmstruct, &
1204 102 : mixed_mo_coeff, mixed_matrix_s, &
1205 102 : density_matrix, energy_qs, w_matrix, mo_eigval)
1206 204 : logger => cp_get_default_logger()
1207 102 : CPASSERT(ASSOCIATED(force_env))
1208 102 : CALL timeset(routineN, handle)
1209 : CALL force_env_get(force_env=force_env, &
1210 102 : force_env_section=force_env_section)
1211 102 : mixed_env => force_env%mixed_env
1212 102 : nforce_eval = SIZE(force_env%sub_force_env)
1213 102 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1214 102 : 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 102 : 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 102 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1227 : ! Get redistributed work matrices
1228 102 : CPASSERT(ASSOCIATED(mixed_cdft))
1229 102 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%mixed_mo_coeff))
1230 102 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%w_matrix))
1231 102 : CPASSERT(ASSOCIATED(mixed_cdft%matrix%mixed_matrix_s))
1232 102 : mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
1233 102 : w_matrix => mixed_cdft%matrix%w_matrix
1234 102 : mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
1235 102 : 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 102 : nvar = SIZE(w_matrix, 2)
1241 102 : nspins = SIZE(mixed_mo_coeff, 2)
1242 : ! Check that the number of MOs/AOs is equal in every CDFT state
1243 408 : ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
1244 298 : DO ispin = 1, nspins
1245 196 : CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=check_mo(1), nrow_global=check_ao(1))
1246 542 : DO iforce_eval = 2, nforce_eval
1247 244 : CALL cp_fm_get_info(mixed_mo_coeff(iforce_eval, ispin), ncol_global=check_mo(2), nrow_global=check_ao(2))
1248 244 : 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 440 : 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 102 : npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
1260 1450 : ALLOCATE (w_matrix_mo(nforce_eval, nforce_eval, nvar))
1261 796 : ALLOCATE (mo_overlap(npermutations), S_det(npermutations, nspins))
1262 816 : ALLOCATE (a(nspins, nvar, npermutations), b(nspins, nvar, npermutations))
1263 102 : a = 0.0_dp
1264 102 : b = 0.0_dp
1265 102 : 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 102 : uniform_occupation = .NOT. ALLOCATED(mixed_cdft%occupations)
1279 102 : should_scale = .FALSE.
1280 102 : 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 298 : 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 196 : NULLIFY (fm_struct_mo, mo_mo_fmstruct)
1300 196 : CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=ncol_mo(ispin), nrow_global=nrow_mo(ispin))
1301 196 : nao = nrow_mo(ispin)
1302 196 : IF (uniform_occupation) THEN
1303 168 : 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 196 : 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 196 : 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 196 : name="ET_TMP_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1357 196 : CALL cp_fm_struct_release(fm_struct_mo)
1358 : CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
1359 196 : name="INVERSE_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1360 : CALL cp_fm_create(matrix=Tinverse, matrix_struct=mo_mo_fmstruct, &
1361 196 : name="T_INVERSE_"//TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1362 556 : 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 556 : TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
1366 : END DO
1367 396 : DO ivar = 1, nvar
1368 844 : DO istate = 1, nforce_eval
1369 1824 : DO jstate = 1, nforce_eval
1370 1176 : 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 1624 : TRIM(ADJUSTL(cp_to_string(ivar)))//"_MATRIX")
1375 : END DO
1376 : END DO
1377 : END DO
1378 196 : CALL cp_fm_struct_release(mo_mo_fmstruct)
1379 : ! Remove empty MOs and (possibly) scale rest with occupation numbers
1380 196 : 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 196 : ipermutation = 0
1398 636 : DO istate = 1, nforce_eval
1399 996 : DO jstate = istate + 1, nforce_eval
1400 360 : ipermutation = ipermutation + 1
1401 : CALL cp_dbcsr_sm_fm_multiply(mixed_matrix_s, mixed_mo_coeff(istate, ispin), &
1402 360 : 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 360 : tmp2, 0.0_dp, mo_overlap(ipermutation))
1406 800 : 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 396 : DO ivar = 1, nvar
1417 844 : DO jstate = 1, nforce_eval
1418 1824 : DO istate = 1, nforce_eval
1419 1176 : 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 728 : 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 1624 : tmp2, 0.0_dp, w_matrix_mo(istate, jstate, ivar))
1427 : END DO
1428 : END DO
1429 : END DO
1430 556 : DO ipermutation = 1, npermutations
1431 : ! Invert and calculate determinant of MO overlaps
1432 360 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1433 360 : 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 332 : S_det(ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
1457 : END IF
1458 360 : 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 928 : DO j = 1, ncol_local
1461 1442 : DO k = 1, nrow_local
1462 1601 : 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 1033 : 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 360 : CALL cp_fm_transpose(inverse_mat, Tinverse)
1471 928 : DO j = 1, ncol_local
1472 1442 : DO k = 1, nrow_local
1473 1601 : 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 1033 : Tinverse%local_data(k, j)
1477 : END DO
1478 : END DO
1479 : END DO
1480 : ! Handle different constraint types
1481 724 : DO ivar = 1, nvar
1482 364 : 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 364 : CPABORT("Unknown constraint type.")
1495 : END SELECT
1496 360 : 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 364 : CPABORT("Unknown constraint type.")
1509 : END SELECT
1510 : END DO
1511 : ! Compute density matrix difference P = P_j - P_i
1512 360 : 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 1088 : CALL force_env%para_env%sum(a(ispin, :, ipermutation))
1518 2004 : CALL force_env%para_env%sum(b(ispin, :, ipermutation))
1519 : END DO
1520 : ! Release work
1521 196 : CALL cp_fm_release(tmp2)
1522 396 : DO ivar = 1, nvar
1523 844 : DO jstate = 1, nforce_eval
1524 1824 : DO istate = 1, nforce_eval
1525 1176 : IF (istate == jstate) CYCLE
1526 1624 : CALL cp_fm_release(w_matrix_mo(istate, jstate, ivar))
1527 : END DO
1528 : END DO
1529 : END DO
1530 556 : DO ipermutation = 1, npermutations
1531 556 : CALL cp_fm_release(mo_overlap(ipermutation))
1532 : END DO
1533 196 : CALL cp_fm_release(Tinverse)
1534 298 : CALL cp_fm_release(inverse_mat)
1535 : END DO
1536 102 : DEALLOCATE (mo_overlap)
1537 102 : DEALLOCATE (w_matrix_mo)
1538 102 : IF (.NOT. uniform_occupation) THEN
1539 14 : DEALLOCATE (homo)
1540 14 : DEALLOCATE (mixed_cdft%occupations)
1541 : END IF
1542 102 : 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 102 : 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 408 : ALLOCATE (Wda(nvar, npermutations))
1552 306 : ALLOCATE (Sda(npermutations))
1553 106 : IF (.NOT. mixed_cdft%identical_constraints) ALLOCATE (Wad(nvar, npermutations))
1554 286 : DO ipermutation = 1, npermutations
1555 184 : IF (nspins == 2) THEN
1556 176 : Sda(ipermutation) = ABS(S_det(ipermutation, 1)*S_det(ipermutation, 2))
1557 : ELSE
1558 8 : Sda(ipermutation) = S_det(ipermutation, 1)**2
1559 : END IF
1560 : ! Finalize <Psi_j | w_i(r) | Psi_i> by multiplication with Sda
1561 472 : DO ivar = 1, nvar
1562 370 : IF (mixed_cdft%identical_constraints) THEN
1563 : Wda(ivar, ipermutation) = (SUM(a(:, ivar, ipermutation)) + SUM(b(:, ivar, ipermutation)))* &
1564 904 : 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 102 : DEALLOCATE (a, b, S_det)
1572 : ! Transfer info about the constraint calculations
1573 816 : ALLOCATE (W_diagonal(nvar, nforce_eval), strength(nvar, nforce_eval), energy(nforce_eval))
1574 102 : W_diagonal = 0.0_dp
1575 330 : DO iforce_eval = 1, nforce_eval
1576 562 : strength(:, iforce_eval) = mixed_env%strength(iforce_eval, :)
1577 : END DO
1578 102 : energy = 0.0_dp
1579 330 : DO iforce_eval = 1, nforce_eval
1580 228 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
1581 190 : 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 166 : CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1585 : END IF
1586 190 : CALL get_qs_env(qs_env, energy=energy_qs, dft_control=dft_control)
1587 292 : IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1588 230 : W_diagonal(:, iforce_eval) = dft_control%qs_control%cdft_control%value(:)
1589 : ! The mixed CDFT Hamiltonian is built from the physical KS energies.
1590 : ! Remove the Lagrange term sum_k lambda_k*(C_k - target_k), which
1591 : ! need not vanish at finite constraint tolerances.
1592 152 : energy(iforce_eval) = energy_qs%total - energy_qs%cdft
1593 : END IF
1594 : END DO
1595 102 : CALL force_env%para_env%sum(W_diagonal)
1596 102 : CALL force_env%para_env%sum(energy)
1597 : CALL mixed_cdft_result_type_set(mixed_cdft%results, Wda=Wda, W_diagonal=W_diagonal, &
1598 102 : energy=energy, strength=strength)
1599 102 : IF (.NOT. mixed_cdft%identical_constraints) CALL mixed_cdft_result_type_set(mixed_cdft%results, Wad=Wad)
1600 : ! Construct S
1601 408 : ALLOCATE (S_mat(nforce_eval, nforce_eval))
1602 330 : DO istate = 1, nforce_eval
1603 330 : S_mat(istate, istate) = 1.0_dp
1604 : END DO
1605 286 : DO ipermutation = 1, npermutations
1606 184 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1607 184 : S_mat(istate, jstate) = Sda(ipermutation)
1608 286 : S_mat(jstate, istate) = Sda(ipermutation)
1609 : END DO
1610 102 : CALL mixed_cdft_result_type_set(mixed_cdft%results, S=S_mat)
1611 : ! Invert S via eigendecomposition and compute S^-(1/2)
1612 408 : ALLOCATE (eigenv(nforce_eval), tmp_mat(nforce_eval, nforce_eval))
1613 102 : CALL diamat_all(S_mat, eigenv, .TRUE.)
1614 102 : tmp_mat = 0.0_dp
1615 330 : DO istate = 1, nforce_eval
1616 228 : IF (eigenv(istate) < 1.0e-14_dp) THEN
1617 : ! Safeguard against division with 0 and negative numbers
1618 10 : eigenv(istate) = 1.0e-14_dp
1619 : CALL cp_warn(__LOCATION__, &
1620 : "The overlap matrix is numerically nearly singular. "// &
1621 10 : "Calculation proceeds but the results might be meaningless.")
1622 : END IF
1623 330 : tmp_mat(istate, istate) = 1.0_dp/SQRT(eigenv(istate))
1624 : END DO
1625 5998 : tmp_mat(:, :) = MATMUL(tmp_mat, TRANSPOSE(S_mat))
1626 10450 : S_mat(:, :) = MATMUL(S_mat, tmp_mat) ! S^(-1/2)
1627 102 : CALL mixed_cdft_result_type_set(mixed_cdft%results, S_minushalf=S_mat)
1628 102 : DEALLOCATE (eigenv, tmp_mat, S_mat)
1629 : ! Construct nonorthogonal diabatic Hamiltonian matrix H''
1630 306 : ALLOCATE (H_mat(nforce_eval, nforce_eval))
1631 120 : IF (mixed_cdft%nonortho_coupling) ALLOCATE (coupling_nonortho(npermutations))
1632 330 : DO istate = 1, nforce_eval
1633 330 : H_mat(istate, istate) = energy(istate)
1634 : END DO
1635 286 : DO ipermutation = 1, npermutations
1636 184 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1637 184 : sum_a = 0.0_dp
1638 184 : sum_b = 0.0_dp
1639 370 : DO ivar = 1, nvar
1640 : ! V_J * <Psi_J | w_J(r) | Psi_J>
1641 186 : sum_b(1) = sum_b(1) + strength(ivar, jstate)*W_diagonal(ivar, jstate)
1642 : ! V_I * <Psi_I | w_I(r) | Psi_I>
1643 186 : sum_a(1) = sum_a(1) + strength(ivar, istate)*W_diagonal(ivar, istate)
1644 370 : IF (mixed_cdft%identical_constraints) THEN
1645 : ! V_J * W_IJ
1646 184 : sum_b(2) = sum_b(2) + strength(ivar, jstate)*Wda(ivar, ipermutation)
1647 : ! V_I * W_JI
1648 184 : sum_a(2) = sum_a(2) + strength(ivar, istate)*Wda(ivar, ipermutation)
1649 : ELSE
1650 : ! V_J * W_IJ
1651 2 : sum_b(2) = sum_b(2) + strength(ivar, jstate)*Wad(ivar, ipermutation)
1652 : ! V_I * W_JI
1653 2 : sum_a(2) = sum_a(2) + strength(ivar, istate)*Wda(ivar, ipermutation)
1654 : END IF
1655 : END DO
1656 : ! 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>
1657 : ! H_IJ = F_J*S_IJ - V_J * W_IJ
1658 184 : c(1) = (energy(jstate) + sum_b(1))*Sda(ipermutation) - sum_b(2)
1659 : ! H_JI = F_I*S_JI - V_I * W_JI
1660 184 : c(2) = (energy(istate) + sum_a(1))*Sda(ipermutation) - sum_a(2)
1661 : ! H''(I,J) = 0.5*(H_IJ+H_JI) = H''(J,I)
1662 184 : H_mat(istate, jstate) = (c(1) + c(2))*0.5_dp
1663 184 : H_mat(jstate, istate) = H_mat(istate, jstate)
1664 470 : IF (mixed_cdft%nonortho_coupling) coupling_nonortho(ipermutation) = H_mat(istate, jstate)
1665 : END DO
1666 102 : CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat)
1667 102 : DEALLOCATE (H_mat, W_diagonal, Wda, strength, energy, Sda)
1668 102 : IF (ALLOCATED(Wad)) DEALLOCATE (Wad)
1669 102 : IF (mixed_cdft%nonortho_coupling) THEN
1670 18 : CALL mixed_cdft_result_type_set(mixed_cdft%results, nonortho=coupling_nonortho)
1671 18 : DEALLOCATE (coupling_nonortho)
1672 : END IF
1673 : ! Compute metric to assess reliability of coupling
1674 102 : IF (mixed_cdft%calculate_metric) CALL mixed_cdft_calculate_metric(force_env, mixed_cdft, density_matrix_diff, ncol_mo)
1675 : ! Compute coupling also with the wavefunction overlap method, see Migliore2009
1676 : ! Requires the unconstrained KS ground state wavefunction as input
1677 102 : IF (mixed_cdft%wfn_overlap_method) THEN
1678 8 : IF (.NOT. uniform_occupation) THEN
1679 : CALL cp_abort(__LOCATION__, &
1680 0 : "Wavefunction overlap method supports only uniformly occupied MOs.")
1681 : END IF
1682 8 : CALL mixed_cdft_wfn_overlap_method(force_env, mixed_cdft, ncol_mo, nrow_mo)
1683 : END IF
1684 : ! Release remaining work
1685 102 : DEALLOCATE (nrow_mo, ncol_mo)
1686 102 : CALL mixed_cdft_work_type_release(mixed_cdft%matrix)
1687 102 : CALL timestop(handle)
1688 :
1689 306 : END SUBROUTINE mixed_cdft_interaction_matrices
1690 :
1691 : ! **************************************************************************************************
1692 : !> \brief Routine to calculate the CDFT electronic couplings.
1693 : !> \param force_env the force_env that holds the CDFT states
1694 : !> \par History
1695 : !> 11.17 created [Nico Holmberg]
1696 : ! **************************************************************************************************
1697 102 : SUBROUTINE mixed_cdft_calculate_coupling_low(force_env)
1698 : TYPE(force_env_type), POINTER :: force_env
1699 :
1700 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_calculate_coupling_low'
1701 :
1702 : INTEGER :: handle, ipermutation, istate, jstate, &
1703 : nforce_eval, npermutations, nvar
1704 : LOGICAL :: use_lowdin, use_rotation
1705 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_lowdin, coupling_rotation, &
1706 102 : eigenv
1707 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, W_mat
1708 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1709 :
1710 102 : NULLIFY (mixed_cdft)
1711 102 : CPASSERT(ASSOCIATED(force_env))
1712 102 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1713 102 : CALL timeset(routineN, handle)
1714 102 : CPASSERT(ASSOCIATED(mixed_cdft))
1715 102 : CPASSERT(ALLOCATED(mixed_cdft%results%W_diagonal))
1716 102 : CPASSERT(ALLOCATED(mixed_cdft%results%Wda))
1717 102 : CPASSERT(ALLOCATED(mixed_cdft%results%S_minushalf))
1718 102 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1719 : ! Decide which methods to use for computing the coupling
1720 : ! Default behavior is to use rotation when a single constraint is active, otherwise uses Lowdin orthogonalization
1721 : ! The latter can also be explicitly requested when a single constraint is active
1722 : ! Possibly computes the coupling additionally with the wavefunction overlap method
1723 102 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
1724 102 : nvar = SIZE(mixed_cdft%results%Wda, 1)
1725 102 : npermutations = nforce_eval*(nforce_eval - 1)/2
1726 408 : ALLOCATE (tmp_mat(nforce_eval, nforce_eval))
1727 102 : IF (nvar == 1 .AND. mixed_cdft%identical_constraints) THEN
1728 98 : use_rotation = .TRUE.
1729 98 : use_lowdin = mixed_cdft%use_lowdin
1730 : ELSE
1731 : use_rotation = .FALSE.
1732 : use_lowdin = .TRUE.
1733 : END IF
1734 : ! Calculate coupling by rotating the CDFT states to eigenstates of the weight matrix W (single constraint only)
1735 : IF (use_rotation) THEN
1736 : ! Construct W
1737 490 : ALLOCATE (W_mat(nforce_eval, nforce_eval), coupling_rotation(npermutations))
1738 294 : ALLOCATE (eigenv(nforce_eval))
1739 : ! W_mat(i, i) = N_i where N_i is the value of the constraint in state i
1740 318 : DO istate = 1, nforce_eval
1741 538 : W_mat(istate, istate) = SUM(mixed_cdft%results%W_diagonal(:, istate))
1742 : END DO
1743 : ! W_mat(i, j) = <Psi_i|w(r)|Psi_j>
1744 278 : DO ipermutation = 1, npermutations
1745 180 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1746 360 : W_mat(istate, jstate) = SUM(mixed_cdft%results%Wda(:, ipermutation))
1747 458 : W_mat(jstate, istate) = W_mat(istate, jstate)
1748 : END DO
1749 : ! Solve generalized eigenvalue equation WV = SVL
1750 : ! Convert to standard eigenvalue problem via symmetric orthogonalisation
1751 5082 : tmp_mat(:, :) = MATMUL(W_mat, mixed_cdft%results%S_minushalf) ! W * S^(-1/2)
1752 5082 : W_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, tmp_mat) ! W' = S^(-1/2) * W * S^(-1/2)
1753 98 : CALL diamat_all(W_mat, eigenv, .TRUE.) ! Solve W'V' = AV'
1754 9266 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, W_mat) ! Reverse transformation V = S^(-1/2) V'
1755 : ! Construct final, orthogonal diabatic Hamiltonian matrix H
1756 9266 : W_mat(:, :) = MATMUL(mixed_cdft%results%H, tmp_mat) ! H'' * V
1757 10262 : W_mat(:, :) = MATMUL(TRANSPOSE(tmp_mat), W_mat) ! H = V^T * H'' * V
1758 278 : DO ipermutation = 1, npermutations
1759 180 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1760 278 : coupling_rotation(ipermutation) = W_mat(istate, jstate)
1761 : END DO
1762 98 : CALL mixed_cdft_result_type_set(mixed_cdft%results, rotation=coupling_rotation)
1763 98 : DEALLOCATE (W_mat, coupling_rotation, eigenv)
1764 : END IF
1765 : ! Calculate coupling by Lowdin orthogonalization
1766 102 : IF (use_lowdin) THEN
1767 60 : ALLOCATE (coupling_lowdin(npermutations))
1768 780 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%H, mixed_cdft%results%S_minushalf) ! H'' * S^(-1/2)
1769 : ! Final orthogonal diabatic Hamiltonian matrix H
1770 580 : tmp_mat(:, :) = MATMUL(mixed_cdft%results%S_minushalf, tmp_mat) ! H = S^(-1/2) * H'' * S^(-1/2)
1771 40 : DO ipermutation = 1, npermutations
1772 20 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1773 40 : coupling_lowdin(ipermutation) = tmp_mat(istate, jstate)
1774 : END DO
1775 20 : CALL mixed_cdft_result_type_set(mixed_cdft%results, lowdin=coupling_lowdin)
1776 20 : DEALLOCATE (coupling_lowdin)
1777 : END IF
1778 102 : DEALLOCATE (tmp_mat)
1779 102 : CALL timestop(handle)
1780 :
1781 204 : END SUBROUTINE mixed_cdft_calculate_coupling_low
1782 :
1783 : ! **************************************************************************************************
1784 : !> \brief Performs a configuration interaction calculation in the basis spanned by the CDFT states.
1785 : !> \param force_env the force_env that holds the CDFT states
1786 : !> \par History
1787 : !> 11.17 created [Nico Holmberg]
1788 : ! **************************************************************************************************
1789 102 : SUBROUTINE mixed_cdft_configuration_interaction(force_env)
1790 : TYPE(force_env_type), POINTER :: force_env
1791 :
1792 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_configuration_interaction'
1793 :
1794 : INTEGER :: handle, info, iounit, istate, ivar, &
1795 : nforce_eval, work_array_size
1796 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenv, work
1797 102 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: H_mat, H_mat_copy, S_mat, S_mat_copy
1798 : REAL(KIND=dp), EXTERNAL :: dnrm2
1799 : TYPE(cp_logger_type), POINTER :: logger
1800 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1801 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
1802 :
1803 : EXTERNAL :: dsygv
1804 :
1805 102 : NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1806 :
1807 102 : CPASSERT(ASSOCIATED(force_env))
1808 102 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1809 102 : CPASSERT(ASSOCIATED(mixed_cdft))
1810 :
1811 102 : IF (.NOT. mixed_cdft%do_ci) RETURN
1812 :
1813 22 : logger => cp_get_default_logger()
1814 22 : CALL timeset(routineN, handle)
1815 : CALL force_env_get(force_env=force_env, &
1816 22 : force_env_section=force_env_section)
1817 22 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1818 22 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1819 :
1820 22 : CPASSERT(ALLOCATED(mixed_cdft%results%S))
1821 22 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1822 22 : nforce_eval = SIZE(mixed_cdft%results%S, 1)
1823 132 : ALLOCATE (S_mat(nforce_eval, nforce_eval), H_mat(nforce_eval, nforce_eval))
1824 66 : ALLOCATE (eigenv(nforce_eval))
1825 182 : S_mat(:, :) = mixed_cdft%results%S(:, :)
1826 182 : H_mat(:, :) = mixed_cdft%results%H(:, :)
1827 : ! Workspace query
1828 22 : ALLOCATE (work(1))
1829 22 : info = 0
1830 110 : ALLOCATE (H_mat_copy(nforce_eval, nforce_eval), S_mat_copy(nforce_eval, nforce_eval))
1831 182 : H_mat_copy(:, :) = H_mat(:, :) ! Need explicit copies because dsygv destroys original values
1832 182 : S_mat_copy(:, :) = S_mat(:, :)
1833 22 : CALL dsygv(1, 'V', 'U', nforce_eval, H_mat_copy, nforce_eval, S_mat_copy, nforce_eval, eigenv, work, -1, info)
1834 22 : work_array_size = NINT(work(1))
1835 22 : DEALLOCATE (H_mat_copy, S_mat_copy)
1836 : ! Allocate work array
1837 22 : DEALLOCATE (work)
1838 66 : ALLOCATE (work(work_array_size))
1839 22 : work = 0.0_dp
1840 : ! Solve Hc = eSc
1841 22 : info = 0
1842 22 : CALL dsygv(1, 'V', 'U', nforce_eval, H_mat, nforce_eval, S_mat, nforce_eval, eigenv, work, work_array_size, info)
1843 22 : IF (info /= 0) THEN
1844 0 : IF (info > nforce_eval) THEN
1845 0 : CPABORT("Matrix S is not positive definite")
1846 : ELSE
1847 0 : CPABORT("Diagonalization of H matrix failed.")
1848 : END IF
1849 : END IF
1850 : ! dsygv returns eigenvectors (stored in columns of H_mat) that are normalized to H^T * S * H = I
1851 : ! Renormalize eigenvectors to H^T * H = I
1852 70 : DO ivar = 1, nforce_eval
1853 182 : H_mat(:, ivar) = H_mat(:, ivar)/dnrm2(nforce_eval, H_mat(:, ivar), 1)
1854 : END DO
1855 22 : DEALLOCATE (work)
1856 22 : IF (iounit > 0) THEN
1857 11 : WRITE (iounit, '(/,T3,A)') '------------------ CDFT Configuration Interaction (CDFT-CI) ------------------'
1858 35 : DO ivar = 1, nforce_eval
1859 24 : IF (ivar == 1) THEN
1860 11 : WRITE (iounit, '(T3,A,T58,(3X,F20.14))') 'Ground state energy:', eigenv(ivar)
1861 : ELSE
1862 13 : WRITE (iounit, '(/,T3,A,I2,A,T58,(3X,F20.14))') 'Excited state (', ivar - 1, ' ) energy:', eigenv(ivar)
1863 : END IF
1864 63 : DO istate = 1, nforce_eval, 2
1865 52 : IF (istate == 1) THEN
1866 : WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
1867 24 : 'Expansion coefficients:', H_mat(istate, ivar), H_mat(istate + 1, ivar)
1868 4 : ELSE IF (istate < nforce_eval) THEN
1869 4 : WRITE (iounit, '(T54,(3X,2F12.6))') H_mat(istate, ivar), H_mat(istate + 1, ivar)
1870 : ELSE
1871 0 : WRITE (iounit, '(T54,(3X,F12.6))') H_mat(istate, ivar)
1872 : END IF
1873 : END DO
1874 : END DO
1875 : WRITE (iounit, '(T3,A)') &
1876 11 : '------------------------------------------------------------------------------'
1877 : END IF
1878 22 : DEALLOCATE (S_mat, H_mat, eigenv)
1879 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1880 22 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1881 22 : CALL timestop(handle)
1882 :
1883 22 : END SUBROUTINE mixed_cdft_configuration_interaction
1884 : ! **************************************************************************************************
1885 : !> \brief Block diagonalizes the mixed CDFT Hamiltonian matrix.
1886 : !> \param force_env the force_env that holds the CDFT states
1887 : !> \par History
1888 : !> 11.17 created [Nico Holmberg]
1889 : !> 01.18 added recursive diagonalization
1890 : !> split to subroutines [Nico Holmberg]
1891 : ! **************************************************************************************************
1892 102 : SUBROUTINE mixed_cdft_block_diag(force_env)
1893 : TYPE(force_env_type), POINTER :: force_env
1894 :
1895 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_block_diag'
1896 :
1897 : INTEGER :: handle, i, iounit, irecursion, j, n, &
1898 : nblk, nforce_eval, nrecursion
1899 : LOGICAL :: ignore_excited
1900 102 : TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1901 102 : TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1902 102 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: H_block, S_block
1903 : TYPE(cp_logger_type), POINTER :: logger
1904 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1905 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
1906 :
1907 : EXTERNAL :: dsygv
1908 :
1909 102 : NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1910 :
1911 102 : CPASSERT(ASSOCIATED(force_env))
1912 102 : CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1913 102 : CPASSERT(ASSOCIATED(mixed_cdft))
1914 :
1915 102 : IF (.NOT. mixed_cdft%block_diagonalize) RETURN
1916 :
1917 8 : logger => cp_get_default_logger()
1918 8 : CALL timeset(routineN, handle)
1919 :
1920 8 : CPASSERT(ALLOCATED(mixed_cdft%results%S))
1921 8 : CPASSERT(ALLOCATED(mixed_cdft%results%H))
1922 8 : nforce_eval = SIZE(mixed_cdft%results%S, 1)
1923 :
1924 : CALL force_env_get(force_env=force_env, &
1925 8 : force_env_section=force_env_section)
1926 8 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1927 8 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1928 : ! Read block definitions from input
1929 8 : CALL mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
1930 8 : nblk = SIZE(blocks)
1931 : ! Start block diagonalization
1932 18 : DO irecursion = 1, nrecursion
1933 : ! Print block definitions
1934 10 : IF (iounit > 0 .AND. irecursion == 1) THEN
1935 4 : WRITE (iounit, '(/,T3,A)') '-------------------------- CDFT BLOCK DIAGONALIZATION ------------------------'
1936 4 : WRITE (iounit, '(T3,A)') 'Block diagonalizing the mixed CDFT Hamiltonian'
1937 4 : WRITE (iounit, '(T3,A,I3)') 'Number of blocks:', nblk
1938 4 : WRITE (iounit, '(T3,A,L3)') 'Ignoring excited states within blocks:', ignore_excited
1939 4 : WRITE (iounit, '(/,T3,A)') 'List of CDFT states for each block'
1940 14 : DO i = 1, nblk
1941 33 : WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array))
1942 : END DO
1943 : END IF
1944 : ! Recursive diagonalization: update counters and references
1945 10 : IF (irecursion > 1) THEN
1946 2 : nblk = nblk/2
1947 10 : ALLOCATE (blocks(nblk))
1948 2 : j = 1
1949 6 : DO i = 1, nblk
1950 4 : NULLIFY (blocks(i)%array)
1951 4 : ALLOCATE (blocks(i)%array(2))
1952 12 : blocks(i)%array = [j, j + 1]
1953 6 : j = j + 2
1954 : END DO
1955 : ! Print info
1956 2 : IF (iounit > 0) THEN
1957 1 : WRITE (iounit, '(/, T3,A)') 'Recursive block diagonalization of the mixed CDFT Hamiltonian'
1958 1 : WRITE (iounit, '(T6,A)') 'Block diagonalization is continued until only two matrix blocks remain.'
1959 1 : WRITE (iounit, '(T6,A)') 'The new blocks are formed by collecting pairs of blocks from the previous'
1960 1 : WRITE (iounit, '(T6,A)') 'block diagonalized matrix in ascending order.'
1961 1 : WRITE (iounit, '(/,T3,A,I3,A,I3)') 'Recursion step:', irecursion - 1, ' of ', nrecursion - 1
1962 1 : WRITE (iounit, '(/,T3,A)') 'List of old block indices for each new block'
1963 3 : DO i = 1, nblk
1964 7 : WRITE (iounit, '(T6,A,I3,A,6I3)') 'Block', i, ':', (blocks(i)%array(j), j=1, SIZE(blocks(i)%array))
1965 : END DO
1966 : END IF
1967 : END IF
1968 : ! Get the Hamiltonian and overlap matrices of each block
1969 10 : CALL mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block)
1970 : ! Diagonalize blocks
1971 10 : CALL mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues)
1972 : ! Assemble the block diagonalized matrices
1973 10 : IF (ignore_excited) THEN
1974 8 : n = nblk
1975 : ELSE
1976 2 : n = nforce_eval
1977 : END IF
1978 10 : CALL mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, n, iounit)
1979 : ! Deallocate work
1980 34 : DO i = 1, nblk
1981 24 : DEALLOCATE (H_block(i)%array)
1982 24 : DEALLOCATE (S_block(i)%array)
1983 24 : DEALLOCATE (eigenvalues(i)%array)
1984 34 : DEALLOCATE (blocks(i)%array)
1985 : END DO
1986 18 : DEALLOCATE (H_block, S_block, eigenvalues, blocks)
1987 : END DO ! recursion
1988 8 : IF (iounit > 0) THEN
1989 : WRITE (iounit, '(T3,A)') &
1990 4 : '------------------------------------------------------------------------------'
1991 : END IF
1992 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1993 8 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1994 8 : CALL timestop(handle)
1995 :
1996 102 : END SUBROUTINE mixed_cdft_block_diag
1997 : ! **************************************************************************************************
1998 : !> \brief Routine to calculate the CDFT electronic coupling reliability metric
1999 : !> \param force_env the force_env that holds the CDFT states
2000 : !> \param mixed_cdft the mixed_cdft env
2001 : !> \param density_matrix_diff array holding difference density matrices (P_j - P_i) for every CDFT
2002 : !> state permutation
2003 : !> \param ncol_mo the number of MOs per spin
2004 : !> \par History
2005 : !> 11.17 created [Nico Holmberg]
2006 : ! **************************************************************************************************
2007 14 : SUBROUTINE mixed_cdft_calculate_metric(force_env, mixed_cdft, density_matrix_diff, ncol_mo)
2008 : TYPE(force_env_type), POINTER :: force_env
2009 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2010 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix_diff
2011 : INTEGER, DIMENSION(:) :: ncol_mo
2012 :
2013 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_calculate_metric'
2014 :
2015 : INTEGER :: handle, ipermutation, ispin, j, &
2016 : nforce_eval, npermutations, nspins
2017 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
2018 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: metric
2019 : TYPE(dbcsr_type) :: e_vectors
2020 :
2021 14 : CALL timeset(routineN, handle)
2022 14 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
2023 14 : npermutations = nforce_eval*(nforce_eval - 1)/2
2024 14 : nspins = SIZE(density_matrix_diff, 2)
2025 56 : ALLOCATE (metric(npermutations, nspins))
2026 14 : metric = 0.0_dp
2027 14 : CALL dbcsr_create(e_vectors, template=density_matrix_diff(1, 1)%matrix)
2028 42 : DO ispin = 1, nspins
2029 84 : ALLOCATE (evals(ncol_mo(ispin)))
2030 64 : DO ipermutation = 1, npermutations
2031 : ! Take into account doubly occupied orbitals without LSD
2032 36 : IF (nspins == 1) THEN
2033 0 : CALL dbcsr_scale(density_matrix_diff(ipermutation, 1)%matrix, alpha_scalar=0.5_dp)
2034 : END IF
2035 : ! Diagonalize difference density matrix
2036 : CALL cp_dbcsr_syevd(density_matrix_diff(ipermutation, ispin)%matrix, e_vectors, evals, &
2037 36 : para_env=force_env%para_env, blacs_env=mixed_cdft%blacs_env)
2038 36 : CALL dbcsr_release_p(density_matrix_diff(ipermutation, ispin)%matrix)
2039 100 : DO j = 1, ncol_mo(ispin)
2040 72 : metric(ipermutation, ispin) = metric(ipermutation, ispin) + (evals(j)**2 - evals(j)**4)
2041 : END DO
2042 : END DO
2043 42 : DEALLOCATE (evals)
2044 : END DO
2045 14 : CALL dbcsr_release(e_vectors)
2046 14 : DEALLOCATE (density_matrix_diff)
2047 78 : metric(:, :) = metric(:, :)/4.0_dp
2048 14 : CALL mixed_cdft_result_type_set(mixed_cdft%results, metric=metric)
2049 14 : DEALLOCATE (metric)
2050 14 : CALL timestop(handle)
2051 :
2052 28 : END SUBROUTINE mixed_cdft_calculate_metric
2053 :
2054 : ! **************************************************************************************************
2055 : !> \brief Routine to calculate the electronic coupling according to the wavefunction overlap method
2056 : !> \param force_env the force_env that holds the CDFT states
2057 : !> \param mixed_cdft the mixed_cdft env
2058 : !> \param ncol_mo the number of MOs per spin
2059 : !> \param nrow_mo the number of AOs per spin
2060 : !> \par History
2061 : !> 11.17 created [Nico Holmberg]
2062 : ! **************************************************************************************************
2063 8 : SUBROUTINE mixed_cdft_wfn_overlap_method(force_env, mixed_cdft, ncol_mo, nrow_mo)
2064 : TYPE(force_env_type), POINTER :: force_env
2065 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2066 : INTEGER, DIMENSION(:) :: ncol_mo, nrow_mo
2067 :
2068 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_cdft_wfn_overlap_method'
2069 :
2070 : CHARACTER(LEN=default_path_length) :: file_name
2071 : INTEGER :: handle, ipermutation, ispin, istate, &
2072 : jstate, nao, nforce_eval, nmo, &
2073 : npermutations, nspins
2074 : LOGICAL :: exist, natom_mismatch
2075 : REAL(KIND=dp) :: energy_diff, maxocc, Sda
2076 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coupling_wfn
2077 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: overlaps
2078 8 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2079 : TYPE(cp_fm_struct_type), POINTER :: mo_mo_fmstruct
2080 : TYPE(cp_fm_type) :: inverse_mat, mo_overlap_wfn, mo_tmp
2081 8 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
2082 : TYPE(cp_logger_type), POINTER :: logger
2083 : TYPE(cp_subsys_type), POINTER :: subsys_mix
2084 : TYPE(dbcsr_type), POINTER :: mixed_matrix_s
2085 8 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_set
2086 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2087 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2088 : TYPE(section_vals_type), POINTER :: force_env_section, mixed_cdft_section
2089 :
2090 8 : NULLIFY (mixed_cdft_section, subsys_mix, particle_set, qs_kind_set, atomic_kind_set, &
2091 8 : mixed_mo_coeff, mixed_matrix_s, force_env_section)
2092 16 : logger => cp_get_default_logger()
2093 :
2094 8 : CALL timeset(routineN, handle)
2095 8 : nforce_eval = SIZE(mixed_cdft%results%H, 1)
2096 8 : npermutations = nforce_eval*(nforce_eval - 1)/2
2097 8 : nspins = SIZE(nrow_mo)
2098 8 : mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
2099 8 : mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
2100 : CALL force_env_get(force_env=force_env, &
2101 8 : force_env_section=force_env_section)
2102 : ! Create mo_set for input wfn
2103 40 : ALLOCATE (mo_set(nspins))
2104 8 : IF (nspins == 2) THEN
2105 8 : maxocc = 1.0_dp
2106 : ELSE
2107 0 : maxocc = 2.0_dp
2108 : END IF
2109 24 : DO ispin = 1, nspins
2110 16 : nao = nrow_mo(ispin)
2111 16 : nmo = ncol_mo(ispin)
2112 : ! Only OT with fully occupied orbitals is implicitly supported
2113 : CALL allocate_mo_set(mo_set(ispin), nao=nao, nmo=nmo, nelectron=INT(maxocc*nmo), &
2114 : n_el_f=REAL(maxocc*nmo, dp), maxocc=maxocc, &
2115 16 : flexible_electron_count=0.0_dp)
2116 16 : CALL set_mo_set(mo_set(ispin), uniform_occupation=.TRUE., homo=nmo)
2117 16 : ALLOCATE (mo_set(ispin)%mo_coeff)
2118 : CALL cp_fm_create(matrix=mo_set(ispin)%mo_coeff, &
2119 : matrix_struct=mixed_mo_coeff(1, ispin)%matrix_struct, &
2120 16 : name="GS_MO_COEFF"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX")
2121 48 : ALLOCATE (mo_set(ispin)%eigenvalues(nmo))
2122 40 : ALLOCATE (mo_set(ispin)%occupation_numbers(nmo))
2123 : END DO
2124 : ! Read wfn file (note we assume that the basis set is the same)
2125 8 : IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
2126 : ! This really shouldnt be a problem?
2127 : CALL cp_abort(__LOCATION__, &
2128 0 : "QMMM + wavefunction overlap method not supported.")
2129 : END IF
2130 8 : CALL force_env_get(force_env=force_env, subsys=subsys_mix)
2131 8 : mixed_cdft_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT")
2132 8 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
2133 8 : CPASSERT(ASSOCIATED(mixed_cdft%qs_kind_set))
2134 8 : IF (force_env%para_env%is_source()) THEN
2135 4 : CALL wfn_restart_file_name(file_name, exist, mixed_cdft_section, logger)
2136 : END IF
2137 8 : CALL force_env%para_env%bcast(exist)
2138 8 : CALL force_env%para_env%bcast(file_name)
2139 8 : IF (.NOT. exist) THEN
2140 : CALL cp_abort(__LOCATION__, &
2141 : "User requested to restart the wavefunction from the file named: "// &
2142 : TRIM(file_name)//". This file does not exist. Please check the existence of"// &
2143 : " the file or change properly the value of the keyword WFN_RESTART_FILE_NAME in"// &
2144 0 : " section FORCE_EVAL\MIXED\MIXED_CDFT.")
2145 : END IF
2146 : CALL read_mo_set_from_restart(mo_array=mo_set, qs_kind_set=mixed_cdft%qs_kind_set, particle_set=particle_set, &
2147 : para_env=force_env%para_env, id_nr=0, multiplicity=mixed_cdft%multiplicity, &
2148 : dft_section=mixed_cdft_section, natom_mismatch=natom_mismatch, &
2149 8 : cdft=.TRUE.)
2150 8 : IF (natom_mismatch) THEN
2151 : CALL cp_abort(__LOCATION__, &
2152 0 : "Restart wfn file has a wrong number of atoms")
2153 : END IF
2154 : ! Orthonormalize wfn
2155 24 : DO ispin = 1, nspins
2156 24 : IF (mixed_cdft%has_unit_metric) THEN
2157 0 : CALL make_basis_simple(mo_set(ispin)%mo_coeff, ncol_mo(ispin))
2158 : ELSE
2159 16 : CALL make_basis_sm(mo_set(ispin)%mo_coeff, ncol_mo(ispin), mixed_matrix_s)
2160 : END IF
2161 : END DO
2162 : ! Calculate MO overlaps between reference state (R) and CDFT state pairs I/J
2163 24 : ALLOCATE (coupling_wfn(npermutations))
2164 32 : ALLOCATE (overlaps(2, npermutations, nspins))
2165 8 : overlaps = 0.0_dp
2166 24 : DO ispin = 1, nspins
2167 : ! Allocate work
2168 16 : nao = nrow_mo(ispin)
2169 16 : nmo = ncol_mo(ispin)
2170 : CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, &
2171 16 : context=mixed_cdft%blacs_env, para_env=force_env%para_env)
2172 : CALL cp_fm_create(matrix=mo_overlap_wfn, matrix_struct=mo_mo_fmstruct, &
2173 16 : name="MO_OVERLAP_MATRIX_WFN")
2174 : CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
2175 16 : name="INVERSE_MO_OVERLAP_MATRIX_WFN")
2176 16 : CALL cp_fm_struct_release(mo_mo_fmstruct)
2177 : CALL cp_fm_create(matrix=mo_tmp, &
2178 : matrix_struct=mixed_mo_coeff(1, ispin)%matrix_struct, &
2179 16 : name="OVERLAP_MO_COEFF_WFN")
2180 40 : DO ipermutation = 1, npermutations
2181 24 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
2182 : ! S*C_r
2183 : CALL cp_dbcsr_sm_fm_multiply(mixed_matrix_s, mo_set(ispin)%mo_coeff, &
2184 24 : mo_tmp, nmo, 1.0_dp, 0.0_dp)
2185 : ! C_i^T * (S*C_r)
2186 24 : CALL cp_fm_set_all(mo_overlap_wfn, alpha=0.0_dp)
2187 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
2188 : mixed_mo_coeff(istate, ispin), &
2189 24 : mo_tmp, 0.0_dp, mo_overlap_wfn)
2190 24 : CALL cp_fm_invert(mo_overlap_wfn, inverse_mat, overlaps(1, ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
2191 : ! C_j^T * (S*C_r)
2192 24 : CALL cp_fm_set_all(mo_overlap_wfn, alpha=0.0_dp)
2193 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
2194 : mixed_mo_coeff(jstate, ispin), &
2195 24 : mo_tmp, 0.0_dp, mo_overlap_wfn)
2196 64 : CALL cp_fm_invert(mo_overlap_wfn, inverse_mat, overlaps(2, ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
2197 : END DO
2198 16 : CALL cp_fm_release(mo_overlap_wfn)
2199 16 : CALL cp_fm_release(inverse_mat)
2200 16 : CALL cp_fm_release(mo_tmp)
2201 56 : CALL deallocate_mo_set(mo_set(ispin))
2202 : END DO
2203 8 : DEALLOCATE (mo_set)
2204 20 : DO ipermutation = 1, npermutations
2205 12 : CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
2206 12 : IF (nspins == 2) THEN
2207 12 : overlaps(1, ipermutation, 1) = ABS(overlaps(1, ipermutation, 1)*overlaps(1, ipermutation, 2)) ! A in eq. 12c
2208 12 : overlaps(2, ipermutation, 1) = ABS(overlaps(2, ipermutation, 1)*overlaps(2, ipermutation, 2)) ! B in eq. 12c
2209 : ELSE
2210 0 : overlaps(1, ipermutation, 1) = overlaps(1, ipermutation, 1)**2
2211 0 : overlaps(2, ipermutation, 1) = overlaps(2, ipermutation, 1)**2
2212 : END IF
2213 : ! Calculate coupling using eq. 12c
2214 : ! The coupling is singular if A = B (i.e. states I/J are identical or charge in ground state is fully delocalized)
2215 32 : IF (ABS(overlaps(1, ipermutation, 1) - overlaps(2, ipermutation, 1)) <= 1.0e-14_dp) THEN
2216 : CALL cp_warn(__LOCATION__, &
2217 : "Coupling between states is singular and set to zero. "// &
2218 : "Potential causes: coupling is computed between identical CDFT states or the spin/charge "// &
2219 2 : "density is fully delocalized in the unconstrained ground state.")
2220 2 : coupling_wfn(ipermutation) = 0.0_dp
2221 : ELSE
2222 10 : energy_diff = mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate)
2223 10 : Sda = mixed_cdft%results%S(istate, jstate)
2224 : coupling_wfn(ipermutation) = ABS((overlaps(1, ipermutation, 1)*overlaps(2, ipermutation, 1)/ &
2225 : (overlaps(1, ipermutation, 1)**2 - overlaps(2, ipermutation, 1)**2))* &
2226 : (energy_diff)/(1.0_dp - Sda**2)* &
2227 : (1.0_dp - (overlaps(1, ipermutation, 1)**2 + overlaps(2, ipermutation, 1)**2)/ &
2228 : (2.0_dp*overlaps(1, ipermutation, 1)*overlaps(2, ipermutation, 1))* &
2229 10 : Sda))
2230 : END IF
2231 : END DO
2232 8 : DEALLOCATE (overlaps)
2233 8 : CALL mixed_cdft_result_type_set(mixed_cdft%results, wfn=coupling_wfn)
2234 8 : DEALLOCATE (coupling_wfn)
2235 8 : CALL timestop(handle)
2236 :
2237 16 : END SUBROUTINE mixed_cdft_wfn_overlap_method
2238 :
2239 : ! **************************************************************************************************
2240 : !> \brief Becke constraint adapted to mixed calculations, details in qs_cdft_methods.F
2241 : !> \param force_env the force_env that holds the CDFT states
2242 : !> \param calculate_forces determines if forces should be calculted
2243 : !> \par History
2244 : !> 02.2016 created [Nico Holmberg]
2245 : !> 03.2016 added dynamic load balancing (dlb)
2246 : !> changed pw_p_type data types to rank-3 reals to accommodate dlb
2247 : !> and to reduce overall memory footprint
2248 : !> split to subroutines [Nico Holmberg]
2249 : !> 04.2016 introduced mixed grid mapping [Nico Holmberg]
2250 : ! **************************************************************************************************
2251 36 : SUBROUTINE mixed_becke_constraint(force_env, calculate_forces)
2252 : TYPE(force_env_type), POINTER :: force_env
2253 : LOGICAL, INTENT(IN) :: calculate_forces
2254 :
2255 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint'
2256 :
2257 : INTEGER :: handle
2258 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
2259 : LOGICAL :: in_memory, store_vectors
2260 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint
2261 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients
2262 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: position_vecs, R12
2263 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pair_dist_vecs
2264 : TYPE(cp_logger_type), POINTER :: logger
2265 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2266 : TYPE(mixed_environment_type), POINTER :: mixed_env
2267 :
2268 36 : NULLIFY (mixed_env, mixed_cdft)
2269 36 : store_vectors = .TRUE.
2270 36 : logger => cp_get_default_logger()
2271 36 : CALL timeset(routineN, handle)
2272 36 : mixed_env => force_env%mixed_env
2273 36 : CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
2274 : CALL mixed_becke_constraint_init(force_env, mixed_cdft, calculate_forces, &
2275 : is_constraint, in_memory, store_vectors, &
2276 : R12, position_vecs, pair_dist_vecs, &
2277 36 : coefficients, catom)
2278 : CALL mixed_becke_constraint_low(force_env, mixed_cdft, in_memory, &
2279 : is_constraint, store_vectors, R12, &
2280 : position_vecs, pair_dist_vecs, &
2281 36 : coefficients, catom)
2282 36 : CALL timestop(handle)
2283 :
2284 36 : END SUBROUTINE mixed_becke_constraint
2285 : ! **************************************************************************************************
2286 : !> \brief Initialize the mixed Becke constraint calculation
2287 : !> \param force_env the force_env that holds the CDFT states
2288 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
2289 : !> \param calculate_forces determines if forces should be calculted
2290 : !> \param is_constraint a list used to determine which atoms in the system define the constraint
2291 : !> \param in_memory decides whether to build the weight function gradients in parallel before solving
2292 : !> the CDFT states or later during the SCF procedure of the individual states
2293 : !> \param store_vectors should temporary arrays be stored in memory to accelerate the calculation
2294 : !> \param R12 temporary array holding the pairwise atomic distances
2295 : !> \param position_vecs temporary array holding the pbc corrected atomic position vectors
2296 : !> \param pair_dist_vecs temporary array holding the pairwise displament vectors
2297 : !> \param coefficients array that determines how atoms should be summed to form the constraint
2298 : !> \param catom temporary array to map the global index of constraint atoms to their position
2299 : !> in a list that holds only constraint atoms
2300 : !> \par History
2301 : !> 03.2016 created [Nico Holmberg]
2302 : ! **************************************************************************************************
2303 36 : SUBROUTINE mixed_becke_constraint_init(force_env, mixed_cdft, calculate_forces, &
2304 : is_constraint, in_memory, store_vectors, &
2305 : R12, position_vecs, pair_dist_vecs, coefficients, &
2306 : catom)
2307 : TYPE(force_env_type), POINTER :: force_env
2308 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2309 : LOGICAL, INTENT(IN) :: calculate_forces
2310 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: is_constraint
2311 : LOGICAL, INTENT(OUT) :: in_memory
2312 : LOGICAL, INTENT(IN) :: store_vectors
2313 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
2314 : INTENT(out) :: R12, position_vecs
2315 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
2316 : INTENT(out) :: pair_dist_vecs
2317 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
2318 : INTENT(OUT) :: coefficients
2319 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(out) :: catom
2320 :
2321 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_init'
2322 :
2323 : CHARACTER(len=2) :: element_symbol
2324 : INTEGER :: atom_a, bounds(2), handle, i, iatom, iex, iforce_eval, ikind, iounit, ithread, j, &
2325 : jatom, katom, my_work, my_work_size, natom, nforce_eval, nkind, np(3), npme, nthread, &
2326 : numexp, offset_dlb, unit_nr
2327 : INTEGER, DIMENSION(2, 3) :: bo, bo_conf
2328 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores, stride
2329 : LOGICAL :: build, mpi_io
2330 : REAL(kind=dp) :: alpha, chi, coef, ircov, jrcov, ra(3), &
2331 : radius, uij
2332 : REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dr, r, r1, shift
2333 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list
2334 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
2335 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2336 : TYPE(cdft_control_type), POINTER :: cdft_control
2337 : TYPE(cell_type), POINTER :: cell
2338 : TYPE(cp_logger_type), POINTER :: logger
2339 : TYPE(cp_subsys_type), POINTER :: subsys_mix
2340 : TYPE(force_env_type), POINTER :: force_env_qs
2341 : TYPE(hirshfeld_type), POINTER :: cavity_env
2342 : TYPE(particle_list_type), POINTER :: particles
2343 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2344 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2345 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2346 : TYPE(realspace_grid_type), POINTER :: rs_cavity
2347 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
2348 :
2349 36 : NULLIFY (pab, cell, force_env_qs, particle_set, force_env_section, print_section, &
2350 36 : qs_kind_set, particles, subsys_mix, rs_cavity, cavity_env, auxbas_pw_pool, &
2351 36 : atomic_kind_set, radii_list, cdft_control)
2352 72 : logger => cp_get_default_logger()
2353 36 : nforce_eval = SIZE(force_env%sub_force_env)
2354 36 : CALL timeset(routineN, handle)
2355 36 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2356 36 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
2357 : CALL force_env_get(force_env=force_env, &
2358 : subsys=subsys_mix, &
2359 28 : cell=cell)
2360 : CALL cp_subsys_get(subsys=subsys_mix, &
2361 : particles=particles, &
2362 28 : particle_set=particle_set)
2363 : ELSE
2364 24 : DO iforce_eval = 1, nforce_eval
2365 16 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
2366 24 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
2367 : END DO
2368 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
2369 : cp_subsys=subsys_mix, &
2370 8 : cell=cell)
2371 : CALL cp_subsys_get(subsys=subsys_mix, &
2372 : particles=particles, &
2373 8 : particle_set=particle_set)
2374 : END IF
2375 36 : natom = SIZE(particles%els)
2376 36 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2377 36 : cdft_control => mixed_cdft%cdft_control
2378 36 : CPASSERT(ASSOCIATED(cdft_control))
2379 36 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cutoffs)) THEN
2380 24 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2381 72 : ALLOCATE (cdft_control%becke_control%cutoffs(natom))
2382 30 : SELECT CASE (cdft_control%becke_control%cutoff_type)
2383 : CASE (becke_cutoff_global)
2384 18 : cdft_control%becke_control%cutoffs(:) = cdft_control%becke_control%rglobal
2385 : CASE (becke_cutoff_element)
2386 18 : IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%becke_control%cutoffs_tmp)) THEN
2387 : CALL cp_abort(__LOCATION__, &
2388 : "Size of keyword BECKE_CONSTRAINT\ELEMENT_CUTOFFS does "// &
2389 0 : "not match number of atomic kinds in the input coordinate file.")
2390 : END IF
2391 54 : DO ikind = 1, SIZE(atomic_kind_set)
2392 36 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2393 90 : DO iatom = 1, katom
2394 36 : atom_a = atom_list(iatom)
2395 72 : cdft_control%becke_control%cutoffs(atom_a) = cdft_control%becke_control%cutoffs_tmp(ikind)
2396 : END DO
2397 : END DO
2398 42 : DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
2399 : END SELECT
2400 : END IF
2401 36 : build = .FALSE.
2402 36 : IF (cdft_control%becke_control%adjust .AND. &
2403 : .NOT. ASSOCIATED(cdft_control%becke_control%aij)) THEN
2404 72 : ALLOCATE (cdft_control%becke_control%aij(natom, natom))
2405 18 : build = .TRUE.
2406 : END IF
2407 108 : ALLOCATE (catom(cdft_control%natoms))
2408 : IF (cdft_control%save_pot .OR. &
2409 : cdft_control%becke_control%cavity_confine .OR. &
2410 36 : cdft_control%becke_control%should_skip .OR. &
2411 : mixed_cdft%first_iteration) THEN
2412 108 : ALLOCATE (is_constraint(natom))
2413 36 : is_constraint = .FALSE.
2414 : END IF
2415 36 : in_memory = calculate_forces .AND. cdft_control%becke_control%in_memory
2416 36 : IF (in_memory .NEQV. calculate_forces) THEN
2417 : CALL cp_abort(__LOCATION__, &
2418 : "The flag BECKE_CONSTRAINT\IN_MEMORY must be activated "// &
2419 0 : "for the calculation of mixed CDFT forces")
2420 : END IF
2421 108 : IF (in_memory .OR. mixed_cdft%first_iteration) ALLOCATE (coefficients(natom))
2422 108 : DO i = 1, cdft_control%natoms
2423 72 : catom(i) = cdft_control%atoms(i)
2424 : IF (cdft_control%save_pot .OR. &
2425 : cdft_control%becke_control%cavity_confine .OR. &
2426 72 : cdft_control%becke_control%should_skip .OR. &
2427 : mixed_cdft%first_iteration) THEN
2428 72 : is_constraint(catom(i)) = .TRUE.
2429 : END IF
2430 108 : IF (in_memory .OR. mixed_cdft%first_iteration) THEN
2431 72 : coefficients(catom(i)) = cdft_control%group(1)%coeff(i)
2432 : END IF
2433 : END DO
2434 36 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
2435 360 : bo = auxbas_pw_pool%pw_grid%bounds_local
2436 144 : np = auxbas_pw_pool%pw_grid%npts
2437 144 : dr = auxbas_pw_pool%pw_grid%dr
2438 144 : shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp
2439 36 : IF (store_vectors) THEN
2440 108 : IF (in_memory) ALLOCATE (pair_dist_vecs(3, natom, natom))
2441 108 : ALLOCATE (position_vecs(3, natom))
2442 : END IF
2443 144 : DO i = 1, 3
2444 144 : cell_v(i) = cell%hmat(i, i)
2445 : END DO
2446 144 : ALLOCATE (R12(natom, natom))
2447 72 : DO iatom = 1, natom - 1
2448 108 : DO jatom = iatom + 1, natom
2449 144 : r = particle_set(iatom)%r
2450 144 : r1 = particle_set(jatom)%r
2451 144 : DO i = 1, 3
2452 108 : r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2453 144 : r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2454 : END DO
2455 144 : dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
2456 36 : IF (store_vectors) THEN
2457 144 : position_vecs(:, iatom) = r(:)
2458 144 : IF (iatom == 1 .AND. jatom == natom) position_vecs(:, jatom) = r1(:)
2459 36 : IF (in_memory) THEN
2460 96 : pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
2461 96 : pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
2462 : END IF
2463 : END IF
2464 144 : R12(iatom, jatom) = NORM2(dist_vec)
2465 36 : R12(jatom, iatom) = R12(iatom, jatom)
2466 72 : IF (build) THEN
2467 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2468 18 : kind_number=ikind)
2469 18 : ircov = cdft_control%becke_control%radii(ikind)
2470 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
2471 18 : kind_number=ikind)
2472 18 : jrcov = cdft_control%becke_control%radii(ikind)
2473 18 : IF (ircov /= jrcov) THEN
2474 18 : chi = ircov/jrcov
2475 18 : uij = (chi - 1.0_dp)/(chi + 1.0_dp)
2476 18 : cdft_control%becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp)
2477 18 : 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 18 : ELSE IF (cdft_control%becke_control%aij(iatom, jatom) &
2481 : < -0.5_dp) THEN
2482 0 : cdft_control%becke_control%aij(iatom, jatom) = -0.5_dp
2483 : END IF
2484 : ELSE
2485 0 : cdft_control%becke_control%aij(iatom, jatom) = 0.0_dp
2486 : END IF
2487 : cdft_control%becke_control%aij(jatom, iatom) = &
2488 18 : -cdft_control%becke_control%aij(iatom, jatom)
2489 : END IF
2490 : END DO
2491 : END DO
2492 : ! Dump some additional information about the calculation
2493 36 : IF (mixed_cdft%first_iteration) THEN
2494 24 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2495 24 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2496 24 : IF (iounit > 0) THEN
2497 : WRITE (iounit, '(/,T3,A,T66)') &
2498 12 : '-------------------------- Becke atomic parameters ---------------------------'
2499 12 : IF (cdft_control%becke_control%adjust) THEN
2500 : WRITE (iounit, '(T3,A,A)') &
2501 9 : 'Atom Element Coefficient', ' Cutoff (angstrom) CDFT Radius (angstrom)'
2502 27 : DO iatom = 1, natom
2503 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2504 : element_symbol=element_symbol, &
2505 18 : kind_number=ikind)
2506 18 : ircov = cp_unit_from_cp2k(cdft_control%becke_control%radii(ikind), "angstrom")
2507 18 : IF (is_constraint(iatom)) THEN
2508 18 : coef = coefficients(iatom)
2509 : ELSE
2510 0 : coef = 0.0_dp
2511 : END IF
2512 : WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3,T73,F8.3)") &
2513 18 : iatom, ADJUSTR(element_symbol), coef, &
2514 18 : cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom"), &
2515 63 : ircov
2516 : END DO
2517 : ELSE
2518 : WRITE (iounit, '(T3,A,A)') &
2519 3 : 'Atom Element Coefficient', ' Cutoff (angstrom)'
2520 9 : DO iatom = 1, natom
2521 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2522 6 : element_symbol=element_symbol)
2523 6 : IF (is_constraint(iatom)) THEN
2524 6 : coef = coefficients(iatom)
2525 : ELSE
2526 0 : coef = 0.0_dp
2527 : END IF
2528 : WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3)") &
2529 6 : iatom, ADJUSTR(element_symbol), coef, &
2530 15 : cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom")
2531 : END DO
2532 : END IF
2533 : WRITE (iounit, '(T3,A)') &
2534 12 : '------------------------------------------------------------------------------'
2535 : END IF
2536 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
2537 24 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2538 24 : mixed_cdft%first_iteration = .FALSE.
2539 : END IF
2540 :
2541 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2542 34 : CPASSERT(ASSOCIATED(mixed_cdft%qs_kind_set))
2543 34 : cavity_env => cdft_control%becke_control%cavity_env
2544 34 : qs_kind_set => mixed_cdft%qs_kind_set
2545 34 : CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2546 34 : nkind = SIZE(qs_kind_set)
2547 34 : IF (.NOT. ASSOCIATED(cavity_env%kind_shape_fn)) THEN
2548 22 : IF (ASSOCIATED(cdft_control%becke_control%radii)) THEN
2549 54 : ALLOCATE (radii_list(SIZE(cdft_control%becke_control%radii)))
2550 54 : DO ikind = 1, SIZE(cdft_control%becke_control%radii)
2551 54 : IF (cavity_env%use_bohr) THEN
2552 0 : radii_list(ikind) = cdft_control%becke_control%radii(ikind)
2553 : ELSE
2554 36 : radii_list(ikind) = cp_unit_from_cp2k(cdft_control%becke_control%radii(ikind), "angstrom")
2555 : END IF
2556 : END DO
2557 : END IF
2558 : CALL create_shape_function(cavity_env, qs_kind_set, atomic_kind_set, &
2559 : radius=cdft_control%becke_control%rcavity, &
2560 22 : radii_list=radii_list)
2561 22 : IF (ASSOCIATED(radii_list)) THEN
2562 18 : DEALLOCATE (radii_list)
2563 : END IF
2564 : END IF
2565 34 : NULLIFY (rs_cavity)
2566 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_rs_grid=rs_cavity, &
2567 34 : auxbas_pw_pool=auxbas_pw_pool)
2568 : ! be careful in parallel nsmax is chosen with multigrid in mind!
2569 34 : CALL rs_grid_zero(rs_cavity)
2570 34 : ALLOCATE (pab(1, 1))
2571 34 : nthread = 1
2572 34 : ithread = 0
2573 102 : DO ikind = 1, SIZE(atomic_kind_set)
2574 68 : numexp = cavity_env%kind_shape_fn(ikind)%numexp
2575 68 : IF (numexp <= 0) CYCLE
2576 68 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2577 204 : ALLOCATE (cores(katom))
2578 136 : DO iex = 1, numexp
2579 68 : alpha = cavity_env%kind_shape_fn(ikind)%zet(iex)
2580 68 : coef = cavity_env%kind_shape_fn(ikind)%coef(iex)
2581 68 : npme = 0
2582 136 : cores = 0
2583 136 : DO iatom = 1, katom
2584 136 : IF (rs_cavity%desc%parallel .AND. .NOT. rs_cavity%desc%distributed) THEN
2585 : ! replicated realspace grid, split the atoms up between procs
2586 68 : IF (MODULO(iatom, rs_cavity%desc%group_size) == rs_cavity%desc%my_pos) THEN
2587 34 : npme = npme + 1
2588 34 : cores(npme) = iatom
2589 : END IF
2590 : ELSE
2591 0 : npme = npme + 1
2592 0 : cores(npme) = iatom
2593 : END IF
2594 : END DO
2595 170 : DO j = 1, npme
2596 34 : iatom = cores(j)
2597 34 : atom_a = atom_list(iatom)
2598 34 : pab(1, 1) = coef
2599 34 : IF (store_vectors) THEN
2600 136 : ra(:) = position_vecs(:, atom_a) + cell_v(:)/2._dp
2601 : ELSE
2602 0 : ra(:) = pbc(particle_set(atom_a)%r, cell)
2603 : END IF
2604 102 : IF (is_constraint(atom_a)) THEN
2605 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
2606 : ra=ra, rb=ra, rp=ra, &
2607 : zetp=alpha, eps=mixed_cdft%eps_rho_rspace, &
2608 : pab=pab, o1=0, o2=0, & ! without map_consistent
2609 34 : prefactor=1.0_dp, cutoff=0.0_dp)
2610 :
2611 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
2612 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, &
2613 : rs_cavity, &
2614 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
2615 : use_subpatch=.TRUE., &
2616 34 : subpatch_pattern=0)
2617 : END IF
2618 : END DO
2619 : END DO
2620 170 : DEALLOCATE (cores)
2621 : END DO
2622 34 : DEALLOCATE (pab)
2623 34 : CALL auxbas_pw_pool%create_pw(cdft_control%becke_control%cavity)
2624 34 : CALL transfer_rs2pw(rs_cavity, cdft_control%becke_control%cavity)
2625 : CALL hfun_zero(cdft_control%becke_control%cavity%array, &
2626 : cdft_control%becke_control%eps_cavity, &
2627 34 : just_zero=.FALSE., bounds=bounds, work=my_work)
2628 34 : IF (bounds(2) < bo(2, 3)) THEN
2629 8 : bounds(2) = bounds(2) - 1
2630 : ELSE
2631 26 : bounds(2) = bo(2, 3)
2632 : END IF
2633 34 : IF (bounds(1) > bo(1, 3)) THEN
2634 : ! In the special case bounds(1) == bounds(2) == bo(2, 3), after this check
2635 : ! bounds(1) > bounds(2) and the subsequent gradient allocation (:, :, :, bounds(1):bounds(2))
2636 : ! will correctly allocate a 0-sized array
2637 8 : bounds(1) = bounds(1) + 1
2638 : ELSE
2639 26 : bounds(1) = bo(1, 3)
2640 : END IF
2641 34 : IF (bounds(1) > bounds(2)) THEN
2642 0 : my_work_size = 0
2643 : ELSE
2644 34 : my_work_size = (bounds(2) - bounds(1) + 1)
2645 34 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2646 0 : my_work_size = my_work_size*(bo(2, 2) - bo(1, 2) + 1)
2647 : ELSE
2648 34 : my_work_size = my_work_size*(bo(2, 1) - bo(1, 1) + 1)
2649 : END IF
2650 : END IF
2651 102 : cdft_control%becke_control%confine_bounds = bounds
2652 34 : IF (cdft_control%becke_control%print_cavity) THEN
2653 : CALL hfun_zero(cdft_control%becke_control%cavity%array, &
2654 0 : cdft_control%becke_control%eps_cavity, just_zero=.TRUE.)
2655 : NULLIFY (stride)
2656 0 : ALLOCATE (stride(3))
2657 0 : stride = [2, 2, 2]
2658 0 : mpi_io = .TRUE.
2659 : unit_nr = cp_print_key_unit_nr(logger, print_section, "", &
2660 : middle_name="BECKE_CAVITY", &
2661 : extension=".cube", file_position="REWIND", &
2662 0 : log_filename=.FALSE., mpi_io=mpi_io)
2663 0 : IF (force_env%para_env%is_source() .AND. unit_nr < 1) THEN
2664 : CALL cp_abort(__LOCATION__, &
2665 0 : "Please turn on PROGRAM_RUN_INFO to print cavity")
2666 : END IF
2667 : CALL cp_pw_to_cube(cdft_control%becke_control%cavity, &
2668 : unit_nr, "CAVITY", particles=particles, &
2669 0 : stride=stride, mpi_io=mpi_io)
2670 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', mpi_io=mpi_io)
2671 0 : DEALLOCATE (stride)
2672 : END IF
2673 : END IF
2674 36 : bo_conf = bo
2675 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2676 102 : bo_conf(:, 3) = cdft_control%becke_control%confine_bounds
2677 : END IF
2678 : ! Load balance
2679 36 : IF (mixed_cdft%dlb) THEN
2680 : CALL mixed_becke_constraint_dlb(force_env, mixed_cdft, my_work, &
2681 8 : my_work_size, natom, bo, bo_conf)
2682 : END IF
2683 : ! The bounds have been finalized => time to allocate storage for working matrices
2684 36 : offset_dlb = 0
2685 36 : IF (mixed_cdft%dlb) THEN
2686 8 : IF (mixed_cdft%dlb_control%send_work .AND. .NOT. mixed_cdft%is_special) THEN
2687 8 : offset_dlb = SUM(mixed_cdft%dlb_control%target_list(2, :))
2688 : END IF
2689 : END IF
2690 36 : IF (cdft_control%becke_control%cavity_confine) THEN
2691 : ! Get rid of the zero part of the confinement cavity (cr3d -> real(:,:,:))
2692 34 : IF (mixed_cdft%is_special) THEN
2693 0 : ALLOCATE (mixed_cdft%sendbuff(SIZE(mixed_cdft%dest_list)))
2694 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2695 : ALLOCATE (mixed_cdft%sendbuff(i)%cavity(mixed_cdft%dest_list_bo(1, i):mixed_cdft%dest_list_bo(2, i), &
2696 0 : bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2697 : mixed_cdft%sendbuff(i)%cavity = cdft_control%becke_control%cavity%array(mixed_cdft%dest_list_bo(1, i): &
2698 : mixed_cdft%dest_list_bo(2, i), &
2699 : bo(1, 2):bo(2, 2), &
2700 0 : bo_conf(1, 3):bo_conf(2, 3))
2701 : END DO
2702 34 : ELSE IF (mixed_cdft%is_pencil) THEN
2703 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)))
2704 : mixed_cdft%cavity = cdft_control%becke_control%cavity%array(bo(1, 1) + offset_dlb:bo(2, 1), &
2705 : bo(1, 2):bo(2, 2), &
2706 0 : bo_conf(1, 3):bo_conf(2, 3))
2707 : ELSE
2708 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)))
2709 : mixed_cdft%cavity = cdft_control%becke_control%cavity%array(bo(1, 1):bo(2, 1), &
2710 : bo(1, 2) + offset_dlb:bo(2, 2), &
2711 3568226 : bo_conf(1, 3):bo_conf(2, 3))
2712 : END IF
2713 34 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
2714 : END IF
2715 36 : IF (mixed_cdft%is_special) THEN
2716 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2717 : ALLOCATE (mixed_cdft%sendbuff(i)%weight(mixed_cdft%dest_list_bo(1, i):mixed_cdft%dest_list_bo(2, i), &
2718 0 : bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2719 0 : mixed_cdft%sendbuff(i)%weight = 0.0_dp
2720 : END DO
2721 36 : ELSE IF (mixed_cdft%is_pencil) THEN
2722 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)))
2723 0 : mixed_cdft%weight = 0.0_dp
2724 : ELSE
2725 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)))
2726 1879397 : mixed_cdft%weight = 0.0_dp
2727 : END IF
2728 36 : IF (in_memory) THEN
2729 24 : IF (mixed_cdft%is_special) THEN
2730 0 : DO i = 1, SIZE(mixed_cdft%dest_list)
2731 : ALLOCATE (mixed_cdft%sendbuff(i)%gradients(3*natom, mixed_cdft%dest_list_bo(1, i): &
2732 : mixed_cdft%dest_list_bo(2, i), &
2733 : bo(1, 2):bo(2, 2), &
2734 0 : bo_conf(1, 3):bo_conf(2, 3)))
2735 0 : mixed_cdft%sendbuff(i)%gradients = 0.0_dp
2736 : END DO
2737 24 : ELSE IF (mixed_cdft%is_pencil) THEN
2738 : ALLOCATE (cdft_control%group(1)%gradients(3*natom, bo(1, 1) + offset_dlb:bo(2, 1), &
2739 : bo(1, 2):bo(2, 2), &
2740 0 : bo_conf(1, 3):bo_conf(2, 3)))
2741 0 : cdft_control%group(1)%gradients = 0.0_dp
2742 : ELSE
2743 : ALLOCATE (cdft_control%group(1)%gradients(3*natom, bo(1, 1):bo(2, 1), &
2744 : bo(1, 2) + offset_dlb:bo(2, 2), &
2745 144 : bo_conf(1, 3):bo_conf(2, 3)))
2746 9808936 : cdft_control%group(1)%gradients = 0.0_dp
2747 : END IF
2748 : END IF
2749 :
2750 36 : CALL timestop(handle)
2751 :
2752 36 : END SUBROUTINE mixed_becke_constraint_init
2753 :
2754 : ! **************************************************************************************************
2755 : !> \brief Setup load balancing for mixed Becke calculation
2756 : !> \param force_env the force_env that holds the CDFT states
2757 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
2758 : !> \param my_work an estimate of the work per processor
2759 : !> \param my_work_size size of the smallest array slice per processor. overloaded processors will
2760 : !> redistribute works as integer multiples of this value.
2761 : !> \param natom the total number of atoms
2762 : !> \param bo bounds of the realspace grid that holds the electron density
2763 : !> \param bo_conf same as bo, but bounds along z-direction have been compacted with confinement
2764 : !> \par History
2765 : !> 03.2016 created [Nico Holmberg]
2766 : ! **************************************************************************************************
2767 8 : SUBROUTINE mixed_becke_constraint_dlb(force_env, mixed_cdft, my_work, &
2768 : my_work_size, natom, bo, bo_conf)
2769 : TYPE(force_env_type), POINTER :: force_env
2770 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
2771 : INTEGER, INTENT(IN) :: my_work, my_work_size, natom
2772 : INTEGER, DIMENSION(2, 3) :: bo, bo_conf
2773 :
2774 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_dlb'
2775 : INTEGER, PARAMETER :: should_deallocate = 7000, &
2776 : uninitialized = -7000
2777 :
2778 : CHARACTER(len=2) :: dummy
2779 : INTEGER :: actually_sent, exhausted_work, handle, i, ind, iounit, ispecial, j, max_targets, &
2780 : more_work, my_pos, my_special_work, my_target, no_overloaded, no_underloaded, nsend, &
2781 : nsend_limit, nsend_max, offset, offset_proc, offset_special, send_total, tags(2)
2782 8 : INTEGER, DIMENSION(:), POINTER :: buffsize, cumulative_work, expected_work, load_imbalance, &
2783 16 : nrecv, nsend_proc, sendbuffer, should_warn, tmp, work_index, work_size
2784 8 : INTEGER, DIMENSION(:, :), POINTER :: targets, tmp_bo
2785 : LOGICAL :: consistent
2786 16 : LOGICAL, DIMENSION(:), POINTER :: mask_recv, mask_send, touched
2787 : REAL(kind=dp) :: average_work, load_scale, &
2788 : very_overloaded, work_factor
2789 8 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cavity
2790 16 : TYPE(buffers_bi), DIMENSION(:), POINTER :: recvbuffer, sbuff
2791 : TYPE(cdft_control_type), POINTER :: cdft_control
2792 : TYPE(cp_logger_type), POINTER :: logger
2793 40 : TYPE(mp_request_type), DIMENSION(4) :: req
2794 8 : TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
2795 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
2796 :
2797 16 : logger => cp_get_default_logger()
2798 8 : CALL timeset(routineN, handle)
2799 8 : mixed_cdft%dlb_control%recv_work = .FALSE.
2800 8 : mixed_cdft%dlb_control%send_work = .FALSE.
2801 8 : NULLIFY (expected_work, work_index, load_imbalance, work_size, &
2802 8 : cumulative_work, sendbuffer, buffsize, req_recv, req_total, &
2803 8 : tmp, nrecv, nsend_proc, targets, tmp_bo, touched, &
2804 8 : mask_recv, mask_send, cavity, recvbuffer, sbuff, force_env_section, &
2805 8 : print_section, cdft_control)
2806 8 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2807 8 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2808 8 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2809 8 : cdft_control => mixed_cdft%cdft_control
2810 : ! These numerical values control data redistribution and are system sensitive
2811 : ! Currently they are not refined during run time which may cause crashes
2812 : ! However, using too many processors or a confinement cavity that is too large relative to the
2813 : ! total system volume are more likely culprits.
2814 8 : load_scale = mixed_cdft%dlb_control%load_scale
2815 8 : very_overloaded = mixed_cdft%dlb_control%very_overloaded
2816 8 : more_work = mixed_cdft%dlb_control%more_work
2817 8 : max_targets = 40
2818 8 : work_factor = 0.8_dp
2819 : ! Reset targets/sources
2820 8 : IF (mixed_cdft%is_special) THEN
2821 0 : DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo, &
2822 0 : mixed_cdft%source_list, mixed_cdft%source_list_bo)
2823 : ALLOCATE (mixed_cdft%dest_list(SIZE(mixed_cdft%dest_list_save)), &
2824 : mixed_cdft%dest_list_bo(SIZE(mixed_cdft%dest_bo_save, 1), SIZE(mixed_cdft%dest_bo_save, 2)), &
2825 : mixed_cdft%source_list(SIZE(mixed_cdft%source_list_save)), &
2826 0 : mixed_cdft%source_list_bo(SIZE(mixed_cdft%source_bo_save, 1), SIZE(mixed_cdft%source_bo_save, 2)))
2827 0 : mixed_cdft%dest_list = mixed_cdft%dest_list_save
2828 0 : mixed_cdft%source_list = mixed_cdft%source_list_save
2829 0 : mixed_cdft%dest_list_bo = mixed_cdft%dest_bo_save
2830 0 : mixed_cdft%source_list_bo = mixed_cdft%source_bo_save
2831 : END IF
2832 : ALLOCATE (mixed_cdft%dlb_control%expected_work(force_env%para_env%num_pe), &
2833 : expected_work(force_env%para_env%num_pe), &
2834 48 : work_size(force_env%para_env%num_pe))
2835 : IF (debug_this_module) THEN
2836 : ALLOCATE (should_warn(force_env%para_env%num_pe))
2837 : should_warn = 0
2838 : END IF
2839 24 : expected_work = 0
2840 8 : expected_work(force_env%para_env%mepos + 1) = my_work
2841 24 : work_size = 0
2842 8 : work_size(force_env%para_env%mepos + 1) = my_work_size
2843 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2844 4 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2845 : work_size(force_env%para_env%mepos + 1) = work_size(force_env%para_env%mepos + 1) - &
2846 : NINT(REAL(mixed_cdft%dlb_control% &
2847 : prediction_error(force_env%para_env%mepos + 1), dp)/ &
2848 0 : REAL(bo(2, 1) - bo(1, 1) + 1, dp))
2849 : ELSE
2850 : work_size(force_env%para_env%mepos + 1) = work_size(force_env%para_env%mepos + 1) - &
2851 : NINT(REAL(mixed_cdft%dlb_control% &
2852 : prediction_error(force_env%para_env%mepos + 1), dp)/ &
2853 4 : REAL(bo(2, 2) - bo(1, 2) + 1, dp))
2854 : END IF
2855 : END IF
2856 40 : CALL force_env%para_env%sum(expected_work)
2857 40 : CALL force_env%para_env%sum(work_size)
2858 : ! We store the unsorted expected work to refine the estimate on subsequent calls to this routine
2859 40 : mixed_cdft%dlb_control%expected_work = expected_work
2860 : ! Take into account the prediction error of the last step
2861 8 : IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2862 20 : expected_work = expected_work - mixed_cdft%dlb_control%prediction_error
2863 : END IF
2864 : !
2865 24 : average_work = REAL(SUM(expected_work), dp)/REAL(force_env%para_env%num_pe, dp)
2866 : ALLOCATE (work_index(force_env%para_env%num_pe), &
2867 : load_imbalance(force_env%para_env%num_pe), &
2868 48 : targets(2, force_env%para_env%num_pe))
2869 40 : load_imbalance = expected_work - NINT(average_work)
2870 8 : no_overloaded = 0
2871 8 : no_underloaded = 0
2872 56 : targets = 0
2873 : ! Convert the load imbalance to a multiple of the actual work size
2874 24 : DO i = 1, force_env%para_env%num_pe
2875 24 : IF (load_imbalance(i) > 0) THEN
2876 8 : no_overloaded = no_overloaded + 1
2877 : ! Allow heavily overloaded processors to dump more data since most likely they have a lot of 'real' work
2878 8 : IF (expected_work(i) > NINT(very_overloaded*average_work)) THEN
2879 0 : load_imbalance(i) = (CEILING(REAL(load_imbalance(i), dp)/REAL(work_size(i), dp)) + more_work)*work_size(i)
2880 : ELSE
2881 8 : load_imbalance(i) = CEILING(REAL(load_imbalance(i), dp)/REAL(work_size(i), dp))*work_size(i)
2882 : END IF
2883 : ELSE
2884 : ! Allow the underloaded processors to take load_scale amount of additional work
2885 : ! otherwise we may be unable to exhaust all overloaded processors
2886 8 : load_imbalance(i) = NINT(load_imbalance(i)*load_scale)
2887 8 : no_underloaded = no_underloaded + 1
2888 : END IF
2889 : END DO
2890 8 : CALL sort(expected_work, force_env%para_env%num_pe, indices=work_index)
2891 : ! Redistribute work in order from the most overloaded processors to the most underloaded processors
2892 : ! Each underloaded processor is limited to one overloaded processor
2893 8 : IF (load_imbalance(force_env%para_env%mepos + 1) > 0) THEN
2894 4 : offset = 0
2895 4 : mixed_cdft%dlb_control%send_work = .TRUE.
2896 : ! Build up the total amount of work that needs redistribution
2897 12 : ALLOCATE (cumulative_work(force_env%para_env%num_pe))
2898 12 : cumulative_work = 0
2899 4 : DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
2900 4 : IF (work_index(i) == force_env%para_env%mepos + 1) THEN
2901 : EXIT
2902 : ELSE
2903 0 : offset = offset + load_imbalance(work_index(i))
2904 0 : IF (i == force_env%para_env%num_pe) THEN
2905 0 : cumulative_work(i) = load_imbalance(work_index(i))
2906 : ELSE
2907 0 : cumulative_work(i) = cumulative_work(i + 1) + load_imbalance(work_index(i))
2908 : END IF
2909 : END IF
2910 : END DO
2911 4 : my_pos = i
2912 4 : j = force_env%para_env%num_pe
2913 4 : nsend_max = load_imbalance(work_index(j))/work_size(work_index(j))
2914 4 : exhausted_work = 0
2915 : ! Determine send offset by going through all processors that are more overloaded than my_pos
2916 4 : DO i = 1, no_underloaded
2917 4 : IF (my_pos == force_env%para_env%num_pe) EXIT
2918 0 : nsend = -load_imbalance(work_index(i))/work_size(work_index(j))
2919 0 : IF (nsend < 1) nsend = 1
2920 0 : nsend_max = nsend_max - nsend
2921 0 : IF (nsend_max < 0) nsend = nsend + nsend_max
2922 0 : exhausted_work = exhausted_work + nsend*work_size(work_index(j))
2923 0 : offset = offset - nsend*work_size(work_index(j))
2924 0 : IF (offset < 0) EXIT
2925 4 : IF (exhausted_work == cumulative_work(j)) THEN
2926 0 : j = j - 1
2927 0 : nsend_max = load_imbalance(work_index(j))/work_size(work_index(j))
2928 : END IF
2929 : END DO
2930 : ! Underloaded processors were fully exhausted: rewind index
2931 : ! Load balancing will fail if this happens on multiple processors
2932 4 : IF (i > no_underloaded) THEN
2933 0 : i = no_underloaded
2934 : END IF
2935 4 : my_target = i
2936 4 : DEALLOCATE (cumulative_work)
2937 : ! Determine how much and who to send slices of my grid points
2938 4 : nsend_max = load_imbalance(force_env%para_env%mepos + 1)/work_size(force_env%para_env%mepos + 1)
2939 : ! This the actual number of available array slices
2940 4 : IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2941 0 : nsend_limit = bo(2, 1) - bo(1, 1) + 1
2942 : ELSE
2943 4 : nsend_limit = bo(2, 2) - bo(1, 2) + 1
2944 : END IF
2945 4 : IF (.NOT. mixed_cdft%is_special) THEN
2946 4 : ALLOCATE (mixed_cdft%dlb_control%target_list(3, max_targets))
2947 : ELSE
2948 0 : ALLOCATE (mixed_cdft%dlb_control%target_list(3 + 2*SIZE(mixed_cdft%dest_list), max_targets))
2949 0 : ALLOCATE (touched(SIZE(mixed_cdft%dest_list)))
2950 0 : touched = .FALSE.
2951 : END IF
2952 644 : mixed_cdft%dlb_control%target_list = uninitialized
2953 4 : i = 1
2954 4 : ispecial = 1
2955 4 : offset_special = 0
2956 4 : targets(1, my_pos) = my_target
2957 4 : send_total = 0
2958 : ! Main loop. Note, we actually allow my_pos to offload more slices than nsend_max
2959 : DO
2960 4 : nsend = -load_imbalance(work_index(my_target))/work_size(force_env%para_env%mepos + 1)
2961 4 : IF (nsend < 1) nsend = 1 ! send at least one block
2962 : ! Prevent over redistribution: leave at least (1-work_factor)*nsend_limit slices to my_pos
2963 4 : IF (nsend > NINT(work_factor*nsend_limit - send_total)) THEN
2964 : nsend = NINT(work_factor*nsend_limit - send_total)
2965 : IF (debug_this_module) THEN
2966 : should_warn(force_env%para_env%mepos + 1) = 1
2967 : END IF
2968 : END IF
2969 4 : mixed_cdft%dlb_control%target_list(1, i) = work_index(my_target) - 1 ! This is the actual processor rank
2970 4 : IF (mixed_cdft%is_special) THEN
2971 0 : mixed_cdft%dlb_control%target_list(2, i) = 0
2972 0 : actually_sent = nsend
2973 0 : DO j = ispecial, SIZE(mixed_cdft%dest_list)
2974 0 : mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + 1
2975 0 : touched(j) = .TRUE.
2976 0 : IF (nsend < mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1) THEN
2977 0 : mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2978 0 : mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(1, j) + nsend - 1
2979 0 : mixed_cdft%dest_list_bo(1, j) = mixed_cdft%dest_list_bo(1, j) + nsend
2980 0 : nsend = 0
2981 0 : EXIT
2982 : ELSE
2983 0 : mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2984 0 : mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(2, j)
2985 0 : nsend = nsend - (mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1)
2986 0 : mixed_cdft%dest_list_bo(1:2, j) = should_deallocate
2987 : END IF
2988 0 : IF (nsend <= 0) EXIT
2989 : END DO
2990 0 : IF (mixed_cdft%dest_list_bo(1, ispecial) == should_deallocate) ispecial = j + 1
2991 0 : actually_sent = actually_sent - nsend
2992 0 : nsend_max = nsend_max - actually_sent
2993 0 : send_total = send_total + actually_sent
2994 : ELSE
2995 4 : mixed_cdft%dlb_control%target_list(2, i) = nsend
2996 4 : nsend_max = nsend_max - nsend
2997 4 : send_total = send_total + nsend
2998 : END IF
2999 4 : IF (nsend_max < 0) nsend_max = 0
3000 4 : IF (nsend_max == 0) EXIT
3001 0 : IF (my_target /= no_underloaded) THEN
3002 0 : my_target = my_target + 1
3003 : ELSE
3004 : ! If multiple processors execute this block load balancing will fail
3005 0 : mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + nsend_max
3006 0 : nsend_max = 0
3007 0 : EXIT
3008 : END IF
3009 0 : i = i + 1
3010 4 : IF (i > max_targets) THEN
3011 : CALL cp_abort(__LOCATION__, &
3012 0 : "Load balancing error: increase max_targets")
3013 : END IF
3014 : END DO
3015 4 : IF (.NOT. mixed_cdft%is_special) THEN
3016 4 : CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3, 1, i)
3017 : ELSE
3018 0 : CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3 + 2*SIZE(mixed_cdft%dest_list), 1, i)
3019 : END IF
3020 4 : targets(2, my_pos) = my_target
3021 : ! Equalize the load on the target processors
3022 4 : IF (.NOT. mixed_cdft%is_special) THEN
3023 4 : IF (send_total > NINT(work_factor*nsend_limit)) send_total = NINT(work_factor*nsend_limit) - 1
3024 4 : nsend = NINT(REAL(send_total, dp)/REAL(SIZE(mixed_cdft%dlb_control%target_list, 2), dp))
3025 8 : mixed_cdft%dlb_control%target_list(2, :) = nsend
3026 : END IF
3027 : ELSE
3028 4 : DO i = 1, no_underloaded
3029 4 : IF (work_index(i) == force_env%para_env%mepos + 1) EXIT
3030 : END DO
3031 : my_pos = i
3032 : END IF
3033 104 : CALL force_env%para_env%sum(targets)
3034 : IF (debug_this_module) THEN
3035 : CALL force_env%para_env%sum(should_warn)
3036 : IF (ANY(should_warn == 1)) THEN
3037 : CALL cp_warn(__LOCATION__, &
3038 : "MIXED_CDFT DLB: Attempted to redistribute more array"// &
3039 : " slices than actually available. Leaving a fraction of the total"// &
3040 : " slices on the overloaded processor. Perhaps you have set LOAD_SCALE too high?")
3041 : END IF
3042 : DEALLOCATE (should_warn)
3043 : END IF
3044 : ! check that there is one-to-one mapping between over- and underloaded processors
3045 8 : IF (force_env%para_env%is_source()) THEN
3046 4 : consistent = .TRUE.
3047 4 : DO i = force_env%para_env%num_pe - 1, force_env%para_env%num_pe - no_overloaded + 1, -1
3048 0 : IF (targets(1, i) > no_underloaded) consistent = .FALSE.
3049 4 : IF (targets(1, i) > targets(2, i + 1)) THEN
3050 : CYCLE
3051 : ELSE
3052 0 : consistent = .FALSE.
3053 : END IF
3054 : END DO
3055 4 : IF (.NOT. consistent) THEN
3056 : IF (debug_this_module .AND. iounit > 0) THEN
3057 : DO i = force_env%para_env%num_pe - 1, force_env%para_env%num_pe - no_overloaded + 1, -1
3058 : WRITE (iounit, '(A,I8,I8,I8,I8,I8)') &
3059 : 'load balancing info', load_imbalance(i), work_index(i), &
3060 : work_size(i), targets(1, i), targets(2, i)
3061 : END DO
3062 : END IF
3063 : CALL cp_abort(__LOCATION__, &
3064 : "Load balancing error: too much data to redistribute."// &
3065 : " Increase LOAD_SCALE or change the number of processors."// &
3066 : " If the confinement cavity occupies a large volume relative"// &
3067 0 : " to the total system volume, it might be worth disabling DLB.")
3068 : END IF
3069 : END IF
3070 : ! Tell the target processors which grid points they should compute
3071 8 : IF (my_pos <= no_underloaded) THEN
3072 4 : DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
3073 4 : IF (targets(1, i) <= my_pos .AND. targets(2, i) >= my_pos) THEN
3074 4 : mixed_cdft%dlb_control%recv_work = .TRUE.
3075 4 : mixed_cdft%dlb_control%my_source = work_index(i) - 1
3076 4 : EXIT
3077 : END IF
3078 : END DO
3079 4 : IF (mixed_cdft%dlb_control%recv_work) THEN
3080 4 : IF (.NOT. mixed_cdft%is_special) THEN
3081 4 : ALLOCATE (mixed_cdft%dlb_control%bo(12))
3082 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%bo, source=mixed_cdft%dlb_control%my_source, &
3083 4 : request=req(1))
3084 4 : CALL req(1)%wait()
3085 12 : mixed_cdft%dlb_control%my_dest_repl = [mixed_cdft%dlb_control%bo(11), mixed_cdft%dlb_control%bo(12)]
3086 12 : mixed_cdft%dlb_control%dest_tags_repl = [mixed_cdft%dlb_control%bo(9), mixed_cdft%dlb_control%bo(10)]
3087 : ALLOCATE (mixed_cdft%dlb_control%cavity(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%weight(mixed_cdft%dlb_control%bo(1):mixed_cdft%dlb_control%bo(2), &
3091 : mixed_cdft%dlb_control%bo(3):mixed_cdft%dlb_control%bo(4), &
3092 20 : mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3093 : ALLOCATE (mixed_cdft%dlb_control%gradients(3*natom, &
3094 : mixed_cdft%dlb_control%bo(1):mixed_cdft%dlb_control%bo(2), &
3095 : mixed_cdft%dlb_control%bo(3):mixed_cdft%dlb_control%bo(4), &
3096 24 : mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3097 22724 : mixed_cdft%dlb_control%gradients = 0.0_dp
3098 3524 : mixed_cdft%dlb_control%weight = 0.0_dp
3099 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%cavity, source=mixed_cdft%dlb_control%my_source, &
3100 4 : request=req(1))
3101 4 : CALL req(1)%wait()
3102 4 : DEALLOCATE (mixed_cdft%dlb_control%bo)
3103 : ELSE
3104 0 : ALLOCATE (buffsize(1))
3105 : CALL force_env%para_env%irecv(msgout=buffsize, source=mixed_cdft%dlb_control%my_source, &
3106 0 : request=req(1))
3107 0 : CALL req(1)%wait()
3108 0 : ALLOCATE (mixed_cdft%dlb_control%bo(12*buffsize(1)))
3109 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%bo, source=mixed_cdft%dlb_control%my_source, &
3110 0 : request=req(1))
3111 0 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(buffsize(1)))
3112 0 : ALLOCATE (req_recv(buffsize(1)))
3113 0 : DEALLOCATE (buffsize)
3114 0 : CALL req(1)%wait()
3115 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
3116 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity(mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3117 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3118 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3119 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3120 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3121 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3122 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%sendbuff(j)%cavity, &
3123 : source=mixed_cdft%dlb_control%my_source, &
3124 0 : request=req_recv(j))
3125 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight(mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3126 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3127 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3128 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3129 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3130 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3131 : ALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients(3*natom, &
3132 : mixed_cdft%dlb_control%bo(12*(j - 1) + 1): &
3133 : mixed_cdft%dlb_control%bo(12*(j - 1) + 2), &
3134 : mixed_cdft%dlb_control%bo(12*(j - 1) + 3): &
3135 : mixed_cdft%dlb_control%bo(12*(j - 1) + 4), &
3136 : mixed_cdft%dlb_control%bo(12*(j - 1) + 7): &
3137 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3138 0 : mixed_cdft%dlb_control%sendbuff(j)%weight = 0.0_dp
3139 0 : mixed_cdft%dlb_control%sendbuff(j)%gradients = 0.0_dp
3140 : mixed_cdft%dlb_control%sendbuff(j)%tag = [mixed_cdft%dlb_control%bo(12*(j - 1) + 9), &
3141 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 10)]
3142 : mixed_cdft%dlb_control%sendbuff(j)%rank = [mixed_cdft%dlb_control%bo(12*(j - 1) + 11), &
3143 0 : mixed_cdft%dlb_control%bo(12*(j - 1) + 12)]
3144 : END DO
3145 0 : CALL mp_waitall(req_recv)
3146 0 : DEALLOCATE (req_recv)
3147 : END IF
3148 : END IF
3149 : ELSE
3150 4 : IF (.NOT. mixed_cdft%is_special) THEN
3151 4 : offset = 0
3152 4 : ALLOCATE (sendbuffer(12))
3153 4 : send_total = 0
3154 8 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3155 : tags = [(i - 1)*3 + 1 + force_env%para_env%mepos*6*max_targets, &
3156 12 : (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets] ! Unique communicator tags
3157 4 : mixed_cdft%dlb_control%target_list(3, i) = tags(1)
3158 4 : IF (mixed_cdft%is_pencil) THEN
3159 : sendbuffer = [bo_conf(1, 1) + offset, &
3160 : bo_conf(1, 1) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3161 : bo_conf(1, 2), bo_conf(2, 2), bo(1, 3), bo(2, 3), bo_conf(1, 3), bo_conf(2, 3), &
3162 0 : tags(1), tags(2), mixed_cdft%dest_list(1), mixed_cdft%dest_list(2)]
3163 : ELSE
3164 : sendbuffer = [bo_conf(1, 1), bo_conf(2, 1), &
3165 : bo_conf(1, 2) + offset, &
3166 : bo_conf(1, 2) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3167 : bo(1, 3), bo(2, 3), bo_conf(1, 3), bo_conf(2, 3), tags(1), tags(2), &
3168 52 : mixed_cdft%dest_list(1), mixed_cdft%dest_list(2)]
3169 : END IF
3170 4 : send_total = send_total + mixed_cdft%dlb_control%target_list(2, i) - 1
3171 : CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dlb_control%target_list(1, i), &
3172 4 : request=req(1))
3173 4 : CALL req(1)%wait()
3174 4 : IF (mixed_cdft%is_pencil) THEN
3175 : ALLOCATE (cavity(bo_conf(1, 1) + offset: &
3176 : bo_conf(1, 1) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3177 0 : bo_conf(1, 2):bo_conf(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
3178 : cavity = cdft_control%becke_control%cavity%array(bo_conf(1, 1) + offset: &
3179 : bo_conf(1, 1) + offset + &
3180 : (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3181 : bo_conf(1, 2):bo_conf(2, 2), &
3182 0 : bo_conf(1, 3):bo_conf(2, 3))
3183 : ELSE
3184 : ALLOCATE (cavity(bo_conf(1, 1):bo_conf(2, 1), &
3185 : bo_conf(1, 2) + offset: &
3186 : bo_conf(1, 2) + offset + (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3187 20 : bo_conf(1, 3):bo_conf(2, 3)))
3188 : cavity = cdft_control%becke_control%cavity%array(bo_conf(1, 1):bo_conf(2, 1), &
3189 : bo_conf(1, 2) + offset: &
3190 : bo_conf(1, 2) + offset + &
3191 : (mixed_cdft%dlb_control%target_list(2, i) - 1), &
3192 7044 : bo_conf(1, 3):bo_conf(2, 3))
3193 : END IF
3194 : CALL force_env%para_env%isend(msgin=cavity, &
3195 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3196 4 : request=req(1))
3197 4 : CALL req(1)%wait()
3198 4 : offset = offset + mixed_cdft%dlb_control%target_list(2, i)
3199 8 : DEALLOCATE (cavity)
3200 : END DO
3201 4 : IF (mixed_cdft%is_pencil) THEN
3202 0 : mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 1)
3203 0 : mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 1) + offset - 1
3204 : ELSE
3205 4 : mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 2)
3206 4 : mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 2) + offset - 1
3207 : END IF
3208 4 : DEALLOCATE (sendbuffer)
3209 : ELSE
3210 0 : ALLOCATE (buffsize(1))
3211 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3212 0 : buffsize = mixed_cdft%dlb_control%target_list(2, i)
3213 : ! Unique communicator tags (dont actually need these, should be removed)
3214 : tags = [(i - 1)*3 + 1 + force_env%para_env%mepos*6*max_targets, &
3215 0 : (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets]
3216 0 : DO j = 4, SIZE(mixed_cdft%dlb_control%target_list, 1)
3217 0 : IF (mixed_cdft%dlb_control%target_list(j, i) > uninitialized) EXIT
3218 : END DO
3219 0 : offset_special = j
3220 0 : offset_proc = j - 4 - (j - 4)/2
3221 : CALL force_env%para_env%isend(msgin=buffsize, &
3222 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3223 0 : request=req(1))
3224 0 : CALL req(1)%wait()
3225 0 : ALLOCATE (sendbuffer(12*buffsize(1)))
3226 0 : DO j = 1, buffsize(1)
3227 : sendbuffer(12*(j - 1) + 1:12*(j - 1) + 12) = [mixed_cdft%dlb_control%target_list(offset_special + 2*(j - 1), i), &
3228 : mixed_cdft%dlb_control%target_list(offset_special + 2*j - 1, i), &
3229 : bo_conf(1, 2), bo_conf(2, 2), bo(1, 3), bo(2, 3), &
3230 : bo_conf(1, 3), bo_conf(2, 3), tags(1), tags(2), &
3231 : mixed_cdft%dest_list(j + offset_proc), &
3232 0 : mixed_cdft%dest_list(j + offset_proc) + force_env%para_env%num_pe/2]
3233 : END DO
3234 : CALL force_env%para_env%isend(msgin=sendbuffer, &
3235 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3236 0 : request=req(1))
3237 0 : CALL req(1)%wait()
3238 0 : DEALLOCATE (sendbuffer)
3239 0 : DO j = 1, buffsize(1)
3240 : ALLOCATE (cavity(mixed_cdft%dlb_control%target_list(offset_special + 2*(j - 1), i): &
3241 : mixed_cdft%dlb_control%target_list(offset_special + 2*j - 1, i), &
3242 0 : bo_conf(1, 2):bo_conf(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
3243 : cavity = cdft_control%becke_control%cavity%array(LBOUND(cavity, 1):UBOUND(cavity, 1), &
3244 : bo_conf(1, 2):bo_conf(2, 2), &
3245 0 : bo_conf(1, 3):bo_conf(2, 3))
3246 : CALL force_env%para_env%isend(msgin=cavity, &
3247 : dest=mixed_cdft%dlb_control%target_list(1, i), &
3248 0 : request=req(1))
3249 0 : CALL req(1)%wait()
3250 0 : DEALLOCATE (cavity)
3251 : END DO
3252 : END DO
3253 0 : DEALLOCATE (buffsize)
3254 : END IF
3255 : END IF
3256 8 : DEALLOCATE (expected_work, work_size, load_imbalance, work_index, targets)
3257 : ! Once calculated, data defined on the distributed grid points is sent directly to the processors that own the
3258 : ! grid points after the constraint is copied onto the two processor groups, instead of sending the data back to
3259 : ! the original owner
3260 8 : IF (mixed_cdft%is_special) THEN
3261 0 : my_special_work = 2
3262 0 : ALLOCATE (mask_send(SIZE(mixed_cdft%dest_list)), mask_recv(SIZE(mixed_cdft%source_list)))
3263 0 : ALLOCATE (nsend_proc(SIZE(mixed_cdft%dest_list)), nrecv(SIZE(mixed_cdft%source_list)))
3264 0 : nrecv = 0
3265 0 : nsend_proc = 0
3266 0 : mask_recv = .FALSE.
3267 0 : mask_send = .FALSE.
3268 : ELSE
3269 : my_special_work = 1
3270 : END IF
3271 72 : ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)), sbuff(SIZE(mixed_cdft%dest_list)))
3272 56 : ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%source_list) + (my_special_work**2)*SIZE(mixed_cdft%dest_list)))
3273 24 : ALLOCATE (mixed_cdft%dlb_control%recv_work_repl(SIZE(mixed_cdft%source_list)))
3274 24 : DO i = 1, SIZE(mixed_cdft%source_list)
3275 16 : NULLIFY (recvbuffer(i)%bv, recvbuffer(i)%iv)
3276 16 : ALLOCATE (recvbuffer(i)%bv(1), recvbuffer(i)%iv(3))
3277 : CALL force_env%para_env%irecv(msgout=recvbuffer(i)%bv, &
3278 : source=mixed_cdft%source_list(i), &
3279 16 : request=req_total(i), tag=1)
3280 24 : IF (mixed_cdft%is_special) THEN
3281 : CALL force_env%para_env%irecv(msgout=recvbuffer(i)%iv, &
3282 : source=mixed_cdft%source_list(i), &
3283 : request=req_total(i + SIZE(mixed_cdft%source_list)), &
3284 0 : tag=2)
3285 : END IF
3286 : END DO
3287 16 : DO i = 1, my_special_work
3288 32 : DO j = 1, SIZE(mixed_cdft%dest_list)
3289 16 : IF (i == 1) THEN
3290 16 : NULLIFY (sbuff(j)%iv, sbuff(j)%bv)
3291 16 : ALLOCATE (sbuff(j)%bv(1))
3292 32 : sbuff(j)%bv = mixed_cdft%dlb_control%send_work
3293 16 : IF (mixed_cdft%is_special) THEN
3294 0 : ALLOCATE (sbuff(j)%iv(3))
3295 0 : sbuff(j)%iv(1:2) = mixed_cdft%dest_list_bo(1:2, j)
3296 0 : sbuff(j)%iv(3) = 0
3297 0 : IF (sbuff(j)%iv(1) == should_deallocate) mask_send(j) = .TRUE.
3298 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3299 0 : sbuff(j)%bv = touched(j)
3300 0 : IF (touched(j)) THEN
3301 0 : nsend = 0
3302 0 : DO ispecial = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3303 0 : IF (mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), ispecial) /= uninitialized) THEN
3304 0 : nsend = nsend + 1
3305 : END IF
3306 : END DO
3307 0 : sbuff(j)%iv(3) = nsend
3308 0 : nsend_proc(j) = nsend
3309 : END IF
3310 : END IF
3311 : END IF
3312 : END IF
3313 16 : ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + my_special_work*SIZE(mixed_cdft%source_list)
3314 : CALL force_env%para_env%isend(msgin=sbuff(j)%bv, &
3315 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
3316 16 : request=req_total(ind), tag=1)
3317 24 : IF (mixed_cdft%is_special) THEN
3318 : CALL force_env%para_env%isend(msgin=sbuff(j)%iv, &
3319 : dest=mixed_cdft%dest_list(j) + (i - 1)*force_env%para_env%num_pe/2, &
3320 0 : request=req_total(ind + 2*SIZE(mixed_cdft%dest_list)), tag=2)
3321 : END IF
3322 : END DO
3323 : END DO
3324 8 : CALL mp_waitall(req_total)
3325 8 : DEALLOCATE (req_total)
3326 24 : DO i = 1, SIZE(mixed_cdft%source_list)
3327 16 : mixed_cdft%dlb_control%recv_work_repl(i) = recvbuffer(i)%bv(1)
3328 16 : IF (mixed_cdft%is_special .AND. mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3329 0 : mixed_cdft%source_list_bo(1:2, i) = recvbuffer(i)%iv(1:2)
3330 0 : nrecv(i) = recvbuffer(i)%iv(3)
3331 0 : IF (recvbuffer(i)%iv(1) == should_deallocate) mask_recv(i) = .TRUE.
3332 : END IF
3333 16 : DEALLOCATE (recvbuffer(i)%bv)
3334 24 : IF (ASSOCIATED(recvbuffer(i)%iv)) DEALLOCATE (recvbuffer(i)%iv)
3335 : END DO
3336 24 : DO j = 1, SIZE(mixed_cdft%dest_list)
3337 16 : DEALLOCATE (sbuff(j)%bv)
3338 24 : IF (ASSOCIATED(sbuff(j)%iv)) DEALLOCATE (sbuff(j)%iv)
3339 : END DO
3340 8 : DEALLOCATE (recvbuffer)
3341 : ! For some reason if debug_this_module is true and is_special is false, the deallocate statement
3342 : ! on line 3433 gets executed no matter what (gfortran 5.3.0 bug?). Printing out the variable seems to fix it...
3343 : IF (debug_this_module) THEN
3344 : WRITE (dummy, *) mixed_cdft%is_special
3345 : END IF
3346 :
3347 8 : IF (.NOT. mixed_cdft%is_special) THEN
3348 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3349 32 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl) + 2))
3350 4 : ALLOCATE (sendbuffer(6))
3351 4 : IF (mixed_cdft%is_pencil) THEN
3352 : sendbuffer = [SIZE(mixed_cdft%dlb_control%target_list, 2), bo_conf(1, 3), bo_conf(2, 3), &
3353 0 : bo_conf(1, 1), bo_conf(1, 2), bo_conf(2, 2)]
3354 : ELSE
3355 : sendbuffer = [SIZE(mixed_cdft%dlb_control%target_list, 2), bo_conf(1, 3), bo_conf(2, 3), &
3356 28 : bo_conf(1, 2), bo_conf(1, 1), bo_conf(2, 1)]
3357 : END IF
3358 8 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3359 24 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl)))
3360 : END IF
3361 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3362 24 : ALLOCATE (mixed_cdft%dlb_control%recv_info(2))
3363 8 : NULLIFY (mixed_cdft%dlb_control%recv_info(1)%target_list, mixed_cdft%dlb_control%recv_info(2)%target_list)
3364 24 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(2))
3365 8 : NULLIFY (mixed_cdft%dlb_control%recvbuff(1)%buffs, mixed_cdft%dlb_control%recvbuff(2)%buffs)
3366 : END IF
3367 : ! First communicate which grid points were distributed
3368 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3369 12 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl) + 1
3370 12 : DO i = 1, 2
3371 : CALL force_env%para_env%isend(msgin=sendbuffer, &
3372 : dest=mixed_cdft%dest_list(i), &
3373 8 : request=req_total(ind))
3374 12 : ind = ind + 1
3375 : END DO
3376 : END IF
3377 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3378 8 : ind = 1
3379 24 : DO i = 1, 2
3380 24 : IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3381 8 : ALLOCATE (mixed_cdft%dlb_control%recv_info(i)%matrix_info(6))
3382 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recv_info(i)%matrix_info, &
3383 : source=mixed_cdft%source_list(i), &
3384 8 : request=req_total(ind))
3385 8 : ind = ind + 1
3386 : END IF
3387 : END DO
3388 : END IF
3389 8 : IF (ASSOCIATED(req_total)) THEN
3390 8 : CALL mp_waitall(req_total)
3391 : END IF
3392 : ! Now communicate which processor handles which grid points
3393 8 : IF (mixed_cdft%dlb_control%send_work) THEN
3394 12 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl) + 1
3395 12 : DO i = 1, 2
3396 8 : IF (i == 2) THEN
3397 8 : mixed_cdft%dlb_control%target_list(3, :) = mixed_cdft%dlb_control%target_list(3, :) + 3*max_targets
3398 : END IF
3399 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%target_list, &
3400 : dest=mixed_cdft%dest_list(i), &
3401 8 : request=req_total(ind))
3402 12 : ind = ind + 1
3403 : END DO
3404 : END IF
3405 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3406 8 : ind = 1
3407 24 : DO i = 1, 2
3408 24 : IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3409 : ALLOCATE (mixed_cdft%dlb_control%recv_info(i)% &
3410 24 : target_list(3, mixed_cdft%dlb_control%recv_info(i)%matrix_info(1)))
3411 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recv_info(i)%target_list, &
3412 : source=mixed_cdft%source_list(i), &
3413 8 : request=req_total(ind))
3414 8 : ind = ind + 1
3415 : END IF
3416 : END DO
3417 : END IF
3418 8 : IF (ASSOCIATED(req_total)) THEN
3419 8 : CALL mp_waitall(req_total)
3420 8 : DEALLOCATE (req_total)
3421 : END IF
3422 8 : IF (ASSOCIATED(sendbuffer)) DEALLOCATE (sendbuffer)
3423 : ELSE
3424 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3425 0 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl) + 2*COUNT(touched)))
3426 0 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3427 0 : ALLOCATE (req_total(COUNT(mixed_cdft%dlb_control%recv_work_repl)))
3428 : END IF
3429 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3430 0 : ind = COUNT(mixed_cdft%dlb_control%recv_work_repl)
3431 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3432 0 : IF (touched(j)) THEN
3433 0 : ALLOCATE (sbuff(j)%iv(4 + 3*nsend_proc(j)))
3434 0 : sbuff(j)%iv(1:4) = [bo_conf(1, 2), bo_conf(2, 2), bo_conf(1, 3), bo_conf(2, 3)]
3435 0 : offset = 5
3436 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3437 0 : IF (mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), i) /= uninitialized) THEN
3438 : sbuff(j)%iv(offset:offset + 2) = [mixed_cdft%dlb_control%target_list(1, i), &
3439 : mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), i), &
3440 0 : mixed_cdft%dlb_control%target_list(4 + 2*j - 1, i)]
3441 0 : offset = offset + 3
3442 : END IF
3443 : END DO
3444 0 : DO ispecial = 1, my_special_work
3445 : CALL force_env%para_env%isend(msgin=sbuff(j)%iv, &
3446 : dest=mixed_cdft%dest_list(j) + (ispecial - 1)*force_env%para_env%num_pe/2, &
3447 0 : request=req_total(ind + ispecial))
3448 : END DO
3449 0 : ind = ind + my_special_work
3450 : END IF
3451 : END DO
3452 : END IF
3453 0 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3454 0 : ALLOCATE (mixed_cdft%dlb_control%recv_info(SIZE(mixed_cdft%source_list)))
3455 0 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(SIZE(mixed_cdft%source_list)))
3456 0 : ind = 1
3457 0 : DO j = 1, SIZE(mixed_cdft%source_list)
3458 : NULLIFY (mixed_cdft%dlb_control%recv_info(j)%target_list, &
3459 0 : mixed_cdft%dlb_control%recvbuff(j)%buffs)
3460 0 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3461 0 : ALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info(4 + 3*nrecv(j)))
3462 : CALL force_env%para_env%irecv(mixed_cdft%dlb_control%recv_info(j)%matrix_info, &
3463 : source=mixed_cdft%source_list(j), &
3464 0 : request=req_total(ind))
3465 0 : ind = ind + 1
3466 : END IF
3467 : END DO
3468 : END IF
3469 0 : IF (ASSOCIATED(req_total)) THEN
3470 0 : CALL mp_waitall(req_total)
3471 0 : DEALLOCATE (req_total)
3472 : END IF
3473 0 : IF (ANY(mask_send)) THEN
3474 : ALLOCATE (tmp(SIZE(mixed_cdft%dest_list) - COUNT(mask_send)), &
3475 0 : tmp_bo(2, SIZE(mixed_cdft%dest_list) - COUNT(mask_send)))
3476 0 : i = 1
3477 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3478 0 : IF (.NOT. mask_send(j)) THEN
3479 0 : tmp(i) = mixed_cdft%dest_list(j)
3480 0 : tmp_bo(1:2, i) = mixed_cdft%dest_list_bo(1:2, j)
3481 0 : i = i + 1
3482 : END IF
3483 : END DO
3484 0 : DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo)
3485 0 : ALLOCATE (mixed_cdft%dest_list(SIZE(tmp)), mixed_cdft%dest_list_bo(2, SIZE(tmp)))
3486 0 : mixed_cdft%dest_list = tmp
3487 0 : mixed_cdft%dest_list_bo = tmp_bo
3488 0 : DEALLOCATE (tmp, tmp_bo)
3489 : END IF
3490 0 : IF (ANY(mask_recv)) THEN
3491 : ALLOCATE (tmp(SIZE(mixed_cdft%source_list) - COUNT(mask_recv)), &
3492 0 : tmp_bo(4, SIZE(mixed_cdft%source_list) - COUNT(mask_recv)))
3493 0 : i = 1
3494 0 : DO j = 1, SIZE(mixed_cdft%source_list)
3495 0 : IF (.NOT. mask_recv(j)) THEN
3496 0 : tmp(i) = mixed_cdft%source_list(j)
3497 0 : tmp_bo(1:4, i) = mixed_cdft%source_list_bo(1:4, j)
3498 0 : i = i + 1
3499 : END IF
3500 : END DO
3501 0 : DEALLOCATE (mixed_cdft%source_list, mixed_cdft%source_list_bo)
3502 0 : ALLOCATE (mixed_cdft%source_list(SIZE(tmp)), mixed_cdft%source_list_bo(4, SIZE(tmp)))
3503 0 : mixed_cdft%source_list = tmp
3504 0 : mixed_cdft%source_list_bo = tmp_bo
3505 0 : DEALLOCATE (tmp, tmp_bo)
3506 : END IF
3507 0 : DEALLOCATE (mask_recv, mask_send)
3508 0 : DEALLOCATE (nsend_proc, nrecv)
3509 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3510 0 : DO j = 1, SIZE(mixed_cdft%dest_list)
3511 0 : IF (touched(j)) DEALLOCATE (sbuff(j)%iv)
3512 : END DO
3513 0 : IF (ASSOCIATED(touched)) DEALLOCATE (touched)
3514 : END IF
3515 : END IF
3516 8 : DEALLOCATE (sbuff)
3517 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
3518 8 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3519 8 : CALL timestop(handle)
3520 :
3521 16 : END SUBROUTINE mixed_becke_constraint_dlb
3522 :
3523 : ! **************************************************************************************************
3524 : !> \brief Low level routine to build mixed Becke constraint and gradients
3525 : !> \param force_env the force_env that holds the CDFT states
3526 : !> \param mixed_cdft container for structures needed to build the mixed CDFT constraint
3527 : !> \param in_memory decides whether to build the weight function gradients in parallel before solving
3528 : !> the CDFT states or later during the SCF procedure of the individual states
3529 : !> \param is_constraint a list used to determine which atoms in the system define the constraint
3530 : !> \param store_vectors should temporary arrays be stored in memory to accelerate the calculation
3531 : !> \param R12 temporary array holding the pairwise atomic distances
3532 : !> \param position_vecs temporary array holding the pbc corrected atomic position vectors
3533 : !> \param pair_dist_vecs temporary array holding the pairwise displament vectors
3534 : !> \param coefficients array that determines how atoms should be summed to form the constraint
3535 : !> \param catom temporary array to map the global index of constraint atoms to their position
3536 : !> in a list that holds only constraint atoms
3537 : !> \par History
3538 : !> 03.2016 created [Nico Holmberg]
3539 : ! **************************************************************************************************
3540 36 : SUBROUTINE mixed_becke_constraint_low(force_env, mixed_cdft, in_memory, &
3541 : is_constraint, store_vectors, R12, position_vecs, &
3542 : pair_dist_vecs, coefficients, catom)
3543 : TYPE(force_env_type), POINTER :: force_env
3544 : TYPE(mixed_cdft_type), POINTER :: mixed_cdft
3545 : LOGICAL, INTENT(IN) :: in_memory
3546 : LOGICAL, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: is_constraint
3547 : LOGICAL, INTENT(IN) :: store_vectors
3548 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
3549 : INTENT(INOUT) :: R12, position_vecs
3550 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
3551 : INTENT(INOUT) :: pair_dist_vecs
3552 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
3553 : INTENT(INOUT) :: coefficients
3554 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: catom
3555 :
3556 : CHARACTER(len=*), PARAMETER :: routineN = 'mixed_becke_constraint_low'
3557 : REAL(KIND=dp), PARAMETER :: eps_sum_cell_f_all = 1.0E-06_dp
3558 :
3559 : INTEGER :: handle, i, iatom, icomm, iforce_eval, index, iounit, ip, ispecial, iwork, j, &
3560 : jatom, jcomm, k, my_special_work, my_work, natom, nbuffs, nforce_eval, np(3), &
3561 : nsent_total, nskipped, nwork, offset, offset_repl
3562 36 : INTEGER, DIMENSION(:), POINTER :: work, work_dlb
3563 36 : INTEGER, DIMENSION(:, :), POINTER :: nsent
3564 : LOGICAL :: completed_recv, should_communicate
3565 36 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: skip_me
3566 36 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: completed
3567 : REAL(kind=dp) :: dist1, dist2, dmyexp, my1, my1_homo, &
3568 : myexp, sum_cell_f_all, &
3569 : sum_cell_f_constr, th, tmp_const
3570 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, distances, ds_dR_i, &
3571 36 : ds_dR_j
3572 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_const_dR, d_sum_Pm_dR, &
3573 36 : distance_vecs, dP_i_dRi
3574 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dP_i_dRj
3575 : REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dR_i, dmy_dR_j, &
3576 : dr, dr1_r2, dr_i_dR, dr_ij_dR, &
3577 : dr_j_dR, grid_p, r, r1, shift
3578 36 : REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
3579 36 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cavity, weight
3580 36 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: gradients
3581 : TYPE(cdft_control_type), POINTER :: cdft_control
3582 : TYPE(cell_type), POINTER :: cell
3583 : TYPE(cp_logger_type), POINTER :: logger
3584 : TYPE(cp_subsys_type), POINTER :: subsys_mix
3585 : TYPE(force_env_type), POINTER :: force_env_qs
3586 36 : TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
3587 36 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_send
3588 : TYPE(particle_list_type), POINTER :: particles
3589 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3590 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3591 : TYPE(section_vals_type), POINTER :: force_env_section, print_section
3592 :
3593 72 : logger => cp_get_default_logger()
3594 36 : NULLIFY (work, req_recv, req_send, work_dlb, nsent, cutoffs, cavity, &
3595 36 : weight, gradients, cell, subsys_mix, force_env_qs, &
3596 36 : particle_set, particles, auxbas_pw_pool, force_env_section, &
3597 36 : print_section, cdft_control)
3598 36 : CALL timeset(routineN, handle)
3599 36 : nforce_eval = SIZE(force_env%sub_force_env)
3600 36 : CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
3601 36 : print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3602 36 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
3603 36 : IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
3604 : CALL force_env_get(force_env=force_env, &
3605 : subsys=subsys_mix, &
3606 28 : cell=cell)
3607 : CALL cp_subsys_get(subsys=subsys_mix, &
3608 : particles=particles, &
3609 28 : particle_set=particle_set)
3610 : ELSE
3611 24 : DO iforce_eval = 1, nforce_eval
3612 16 : IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
3613 24 : force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
3614 : END DO
3615 : CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
3616 : cp_subsys=subsys_mix, &
3617 8 : cell=cell)
3618 : CALL cp_subsys_get(subsys=subsys_mix, &
3619 : particles=particles, &
3620 8 : particle_set=particle_set)
3621 : END IF
3622 36 : natom = SIZE(particles%els)
3623 36 : cdft_control => mixed_cdft%cdft_control
3624 36 : CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
3625 144 : np = auxbas_pw_pool%pw_grid%npts
3626 144 : dr = auxbas_pw_pool%pw_grid%dr
3627 144 : shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp
3628 180 : ALLOCATE (cell_functions(natom), skip_me(natom))
3629 36 : IF (store_vectors) THEN
3630 72 : ALLOCATE (distances(natom))
3631 108 : ALLOCATE (distance_vecs(3, natom))
3632 : END IF
3633 36 : IF (in_memory) THEN
3634 24 : ALLOCATE (ds_dR_j(3))
3635 24 : ALLOCATE (ds_dR_i(3))
3636 72 : ALLOCATE (d_sum_Pm_dR(3, natom))
3637 48 : ALLOCATE (d_sum_const_dR(3, natom))
3638 96 : ALLOCATE (dP_i_dRj(3, natom, natom))
3639 48 : ALLOCATE (dP_i_dRi(3, natom))
3640 24 : th = 1.0e-8_dp
3641 : END IF
3642 36 : IF (mixed_cdft%dlb) THEN
3643 32 : ALLOCATE (work(force_env%para_env%num_pe), work_dlb(force_env%para_env%num_pe))
3644 24 : work = 0
3645 24 : work_dlb = 0
3646 : END IF
3647 36 : my_work = 1
3648 36 : my_special_work = 1
3649 : ! Load balancing: allocate storage for receiving buffers and post recv requests
3650 36 : IF (mixed_cdft%dlb) THEN
3651 8 : IF (mixed_cdft%dlb_control%recv_work) THEN
3652 4 : my_work = 2
3653 4 : IF (.NOT. mixed_cdft%is_special) THEN
3654 40 : ALLOCATE (req_send(2, 3))
3655 : ELSE
3656 0 : ALLOCATE (req_send(2, 3*SIZE(mixed_cdft%dlb_control%sendbuff)))
3657 : END IF
3658 : END IF
3659 16 : IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
3660 8 : IF (.NOT. mixed_cdft%is_special) THEN
3661 8 : offset_repl = 0
3662 8 : IF (mixed_cdft%dlb_control%recv_work_repl(1) .AND. mixed_cdft%dlb_control%recv_work_repl(2)) THEN
3663 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2) + &
3664 0 : SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3665 0 : offset_repl = 3*SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2)
3666 8 : ELSE IF (mixed_cdft%dlb_control%recv_work_repl(1)) THEN
3667 0 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2))))
3668 : ELSE
3669 48 : ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3670 : END IF
3671 : ELSE
3672 0 : nbuffs = 0
3673 0 : offset_repl = 1
3674 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3675 0 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3676 0 : nbuffs = nbuffs + (SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3
3677 : END IF
3678 : END DO
3679 0 : ALLOCATE (req_recv(3*nbuffs))
3680 : END IF
3681 24 : DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3682 24 : IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3683 8 : IF (.NOT. mixed_cdft%is_special) THEN
3684 8 : offset = 0
3685 8 : index = j + (j/2)
3686 64 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)))
3687 16 : DO i = 1, SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)
3688 8 : IF (mixed_cdft%is_pencil) THEN
3689 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3690 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3691 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3692 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3693 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3694 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3695 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3696 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3697 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3698 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3699 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3700 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3701 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3702 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3703 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3704 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3705 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3706 : gradients(3*natom, &
3707 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3708 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3709 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3710 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3711 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3712 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3713 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3714 : ELSE
3715 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3716 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3717 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3718 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3719 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3720 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3721 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3722 40 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3723 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3724 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3725 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3726 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3727 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3728 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3729 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3730 40 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3731 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3732 : gradients(3*natom, &
3733 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(5): &
3734 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(6), &
3735 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset: &
3736 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4) + offset + &
3737 : (mixed_cdft%dlb_control%recv_info(j)%target_list(2, i) - 1), &
3738 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2): &
3739 48 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3)))
3740 : END IF
3741 :
3742 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, &
3743 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3744 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 1), &
3745 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i))
3746 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, &
3747 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3748 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 2), &
3749 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i) + 1)
3750 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, &
3751 : source=mixed_cdft%dlb_control%recv_info(j)%target_list(1, i), &
3752 : request=req_recv(3*(i - 1) + (j - 1)*offset_repl + 3), &
3753 8 : tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i) + 2)
3754 16 : offset = offset + mixed_cdft%dlb_control%recv_info(j)%target_list(2, i)
3755 : END DO
3756 8 : DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info)
3757 : ELSE
3758 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)% &
3759 0 : buffs((SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3))
3760 0 : index = 6
3761 0 : DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
3762 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3763 : weight(mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3764 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3765 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3766 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3767 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3768 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3769 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3770 : cavity(mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3771 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3772 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3773 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3774 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3775 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3776 : ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(i)% &
3777 : gradients(3*natom, mixed_cdft%dlb_control%recv_info(j)%matrix_info(index): &
3778 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(index + 1), &
3779 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(1): &
3780 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(2), &
3781 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(3): &
3782 0 : mixed_cdft%dlb_control%recv_info(j)%matrix_info(4)))
3783 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity, &
3784 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3785 0 : request=req_recv(offset_repl), tag=1)
3786 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, &
3787 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3788 0 : request=req_recv(offset_repl + 1), tag=2)
3789 : CALL force_env%para_env%irecv(msgout=mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients, &
3790 : source=mixed_cdft%dlb_control%recv_info(j)%matrix_info(index - 1), &
3791 0 : request=req_recv(offset_repl + 2), tag=3)
3792 0 : index = index + 3
3793 0 : offset_repl = offset_repl + 3
3794 : END DO
3795 0 : DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info)
3796 : END IF
3797 : END IF
3798 : END DO
3799 : END IF
3800 : END IF
3801 36 : cutoffs => cdft_control%becke_control%cutoffs
3802 36 : should_communicate = .FALSE.
3803 144 : DO i = 1, 3
3804 144 : cell_v(i) = cell%hmat(i, i)
3805 : END DO
3806 76 : DO iwork = my_work, 1, -1
3807 40 : IF (iwork == 2) THEN
3808 4 : IF (.NOT. mixed_cdft%is_special) THEN
3809 4 : cavity => mixed_cdft%dlb_control%cavity
3810 4 : weight => mixed_cdft%dlb_control%weight
3811 4 : gradients => mixed_cdft%dlb_control%gradients
3812 4 : ALLOCATE (completed(2, 3), nsent(2, 3))
3813 : ELSE
3814 0 : my_special_work = SIZE(mixed_cdft%dlb_control%sendbuff)
3815 0 : ALLOCATE (completed(2, 3*my_special_work), nsent(2, 3*my_special_work))
3816 : END IF
3817 4 : completed = .FALSE.
3818 40 : nsent = 0
3819 : ELSE
3820 36 : IF (.NOT. mixed_cdft%is_special) THEN
3821 36 : weight => mixed_cdft%weight
3822 36 : cavity => mixed_cdft%cavity
3823 36 : gradients => cdft_control%group(1)%gradients
3824 : ELSE
3825 0 : my_special_work = SIZE(mixed_cdft%dest_list)
3826 : END IF
3827 : END IF
3828 116 : DO ispecial = 1, my_special_work
3829 40 : nwork = 0
3830 40 : IF (mixed_cdft%is_special) THEN
3831 0 : IF (iwork == 1) THEN
3832 0 : weight => mixed_cdft%sendbuff(ispecial)%weight
3833 0 : cavity => mixed_cdft%sendbuff(ispecial)%cavity
3834 0 : gradients => mixed_cdft%sendbuff(ispecial)%gradients
3835 : ELSE
3836 0 : weight => mixed_cdft%dlb_control%sendbuff(ispecial)%weight
3837 0 : cavity => mixed_cdft%dlb_control%sendbuff(ispecial)%cavity
3838 0 : gradients => mixed_cdft%dlb_control%sendbuff(ispecial)%gradients
3839 : END IF
3840 : END IF
3841 1021 : DO k = LBOUND(weight, 1), UBOUND(weight, 1)
3842 901 : IF (mixed_cdft%dlb .AND. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3843 0 : IF (mixed_cdft%dlb_control%send_work) THEN
3844 0 : IF (k >= mixed_cdft%dlb_control%distributed(1) .AND. &
3845 : k <= mixed_cdft%dlb_control%distributed(2)) THEN
3846 : CYCLE
3847 : END IF
3848 : END IF
3849 : END IF
3850 41952 : DO j = LBOUND(weight, 2), UBOUND(weight, 2)
3851 39209 : IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3852 6400 : IF (mixed_cdft%dlb_control%send_work) THEN
3853 3120 : IF (j >= mixed_cdft%dlb_control%distributed(1) .AND. &
3854 : j <= mixed_cdft%dlb_control%distributed(2)) THEN
3855 : CYCLE
3856 : END IF
3857 : END IF
3858 : END IF
3859 : ! Check if any of the buffers have become available for deallocation
3860 39209 : IF (should_communicate) THEN
3861 2048 : DO icomm = 1, SIZE(nsent, 2)
3862 5120 : DO jcomm = 1, SIZE(nsent, 1)
3863 3072 : IF (nsent(jcomm, icomm) == 1) CYCLE
3864 540 : completed(jcomm, icomm) = req_send(jcomm, icomm)%test()
3865 540 : IF (completed(jcomm, icomm)) THEN
3866 24 : nsent(jcomm, icomm) = nsent(jcomm, icomm) + 1
3867 24 : nsent_total = nsent_total + 1
3868 24 : IF (nsent_total == SIZE(nsent, 1)*SIZE(nsent, 2)) should_communicate = .FALSE.
3869 : END IF
3870 2612 : IF (ALL(completed(:, icomm))) THEN
3871 12 : IF (MODULO(icomm, 3) == 1) THEN
3872 4 : IF (.NOT. mixed_cdft%is_special) THEN
3873 4 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
3874 : ELSE
3875 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%cavity)
3876 : END IF
3877 8 : ELSE IF (MODULO(icomm, 3) == 2) THEN
3878 4 : IF (.NOT. mixed_cdft%is_special) THEN
3879 4 : DEALLOCATE (mixed_cdft%dlb_control%weight)
3880 : ELSE
3881 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%weight)
3882 : END IF
3883 : ELSE
3884 4 : IF (.NOT. mixed_cdft%is_special) THEN
3885 4 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
3886 : ELSE
3887 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%gradients)
3888 : END IF
3889 : END IF
3890 : END IF
3891 : END DO
3892 : END DO
3893 : END IF
3894 : ! Poll to prevent starvation
3895 39209 : IF (ASSOCIATED(req_recv)) THEN
3896 6400 : completed_recv = mp_testall(req_recv)
3897 : END IF
3898 : !
3899 1926389 : DO i = LBOUND(weight, 3), UBOUND(weight, 3)
3900 1807861 : IF (cdft_control%becke_control%cavity_confine) THEN
3901 1716736 : IF (cavity(k, j, i) < cdft_control%becke_control%eps_cavity) CYCLE
3902 : END IF
3903 988401 : grid_p(1) = k*dr(1) + shift(1)
3904 988401 : grid_p(2) = j*dr(2) + shift(2)
3905 988401 : grid_p(3) = i*dr(3) + shift(3)
3906 988401 : nskipped = 0
3907 2965203 : cell_functions = 1.0_dp
3908 988401 : skip_me = .FALSE.
3909 988401 : IF (store_vectors) distances = 0.0_dp
3910 988401 : IF (in_memory) THEN
3911 675116 : d_sum_Pm_dR = 0.0_dp
3912 675116 : d_sum_const_dR = 0.0_dp
3913 675116 : dP_i_dRi = 0.0_dp
3914 : END IF
3915 2554568 : DO iatom = 1, natom
3916 1976802 : IF (skip_me(iatom)) THEN
3917 65541 : cell_functions(iatom) = 0.0_dp
3918 65541 : IF (cdft_control%becke_control%should_skip) THEN
3919 37642 : IF (is_constraint(iatom)) nskipped = nskipped + 1
3920 37642 : IF (nskipped == cdft_control%natoms) THEN
3921 0 : IF (in_memory) THEN
3922 0 : IF (cdft_control%becke_control%cavity_confine) THEN
3923 0 : cavity(k, j, i) = 0.0_dp
3924 : END IF
3925 : END IF
3926 : EXIT
3927 : END IF
3928 : END IF
3929 : CYCLE
3930 : END IF
3931 1911261 : IF (store_vectors) THEN
3932 1911261 : IF (distances(iatom) == 0.0_dp) THEN
3933 6889176 : r = position_vecs(:, iatom)
3934 6889176 : dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
3935 6889176 : dist1 = NORM2(dist_vec)
3936 6889176 : distance_vecs(:, iatom) = dist_vec
3937 1722294 : distances(iatom) = dist1
3938 : ELSE
3939 755868 : dist_vec = distance_vecs(:, iatom)
3940 : dist1 = distances(iatom)
3941 : END IF
3942 : ELSE
3943 0 : r = particle_set(iatom)%r
3944 0 : DO ip = 1, 3
3945 0 : r(ip) = MODULO(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3946 : END DO
3947 0 : dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v
3948 0 : dist1 = NORM2(dist_vec)
3949 : END IF
3950 2489027 : IF (dist1 <= cutoffs(iatom)) THEN
3951 509029 : IF (in_memory) THEN
3952 : IF (dist1 <= th) dist1 = th
3953 1184700 : dr_i_dR(:) = dist_vec(:)/dist1
3954 : END IF
3955 1527087 : DO jatom = 1, natom
3956 1527087 : IF (jatom /= iatom) THEN
3957 509029 : IF (jatom < iatom) THEN
3958 254521 : IF (.NOT. skip_me(jatom)) CYCLE
3959 : END IF
3960 320062 : IF (store_vectors) THEN
3961 320062 : IF (distances(jatom) == 0.0_dp) THEN
3962 1018032 : r1 = position_vecs(:, jatom)
3963 1018032 : dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
3964 1018032 : dist2 = NORM2(dist_vec)
3965 1018032 : distance_vecs(:, jatom) = dist_vec
3966 254508 : distances(jatom) = dist2
3967 : ELSE
3968 262216 : dist_vec = distance_vecs(:, jatom)
3969 : dist2 = distances(jatom)
3970 : END IF
3971 : ELSE
3972 0 : r1 = particle_set(jatom)%r
3973 0 : DO ip = 1, 3
3974 0 : r1(ip) = MODULO(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3975 : END DO
3976 0 : dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v
3977 0 : dist2 = NORM2(dist_vec)
3978 : END IF
3979 320062 : IF (in_memory) THEN
3980 191129 : IF (store_vectors) THEN
3981 764516 : dr1_r2 = pair_dist_vecs(:, iatom, jatom)
3982 : ELSE
3983 0 : dr1_r2 = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
3984 : END IF
3985 : IF (dist2 <= th) dist2 = th
3986 191129 : tmp_const = (R12(iatom, jatom)**3)
3987 764516 : dr_ij_dR(:) = dr1_r2(:)/tmp_const
3988 : !derivativ w.r.t. Rj
3989 764516 : dr_j_dR = dist_vec(:)/dist2
3990 764516 : dmy_dR_j(:) = -(dr_j_dR(:)/R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:))
3991 : !derivativ w.r.t. Ri
3992 764516 : dmy_dR_i(:) = dr_i_dR(:)/R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:)
3993 : END IF
3994 320062 : my1 = (dist1 - dist2)/R12(iatom, jatom)
3995 320062 : IF (cdft_control%becke_control%adjust) THEN
3996 166839 : my1_homo = my1
3997 : my1 = my1 + &
3998 166839 : cdft_control%becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
3999 : END IF
4000 320062 : myexp = 1.5_dp*my1 - 0.5_dp*my1**3
4001 320062 : IF (in_memory) THEN
4002 191129 : dmyexp = 1.5_dp - 1.5_dp*my1**2
4003 : tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
4004 191129 : (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
4005 :
4006 764516 : ds_dR_i(:) = -0.5_dp*tmp_const*dmy_dR_i(:)
4007 764516 : ds_dR_j(:) = -0.5_dp*tmp_const*dmy_dR_j(:)
4008 191129 : IF (cdft_control%becke_control%adjust) THEN
4009 102514 : tmp_const = 1.0_dp - 2.0_dp*my1_homo*cdft_control%becke_control%aij(iatom, jatom)
4010 410056 : ds_dR_i(:) = ds_dR_i(:)*tmp_const
4011 410056 : ds_dR_j(:) = ds_dR_j(:)*tmp_const
4012 : END IF
4013 : END IF
4014 320062 : myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4015 320062 : myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4016 320062 : tmp_const = 0.5_dp*(1.0_dp - myexp)
4017 320062 : cell_functions(iatom) = cell_functions(iatom)*tmp_const
4018 320062 : IF (in_memory) THEN
4019 191129 : IF (ABS(tmp_const) <= th) tmp_const = tmp_const + th
4020 764516 : dP_i_dRi(:, iatom) = dP_i_dRi(:, iatom) + ds_dR_i(:)/tmp_const
4021 764516 : dP_i_dRj(:, iatom, jatom) = ds_dR_j(:)/tmp_const
4022 : END IF
4023 :
4024 320062 : IF (dist2 <= cutoffs(jatom)) THEN
4025 188967 : tmp_const = 0.5_dp*(1.0_dp + myexp)
4026 188967 : cell_functions(jatom) = cell_functions(jatom)*tmp_const
4027 188967 : IF (in_memory) THEN
4028 105046 : IF (ABS(tmp_const) <= th) tmp_const = tmp_const + th
4029 420184 : dP_i_dRj(:, jatom, iatom) = -ds_dR_i(:)/tmp_const
4030 420184 : dP_i_dRi(:, jatom) = dP_i_dRi(:, jatom) - ds_dR_j(:)/tmp_const
4031 : END IF
4032 : ELSE
4033 131095 : skip_me(jatom) = .TRUE.
4034 : END IF
4035 : END IF
4036 : END DO
4037 509029 : IF (in_memory) THEN
4038 1184700 : dP_i_dRi(:, iatom) = cell_functions(iatom)*dP_i_dRi(:, iatom)
4039 1184700 : d_sum_Pm_dR(:, iatom) = d_sum_Pm_dR(:, iatom) + dP_i_dRi(:, iatom)
4040 296175 : IF (is_constraint(iatom)) THEN
4041 : d_sum_const_dR(:, iatom) = d_sum_const_dR(:, iatom) + dP_i_dRi(:, iatom)* &
4042 1184700 : coefficients(iatom)
4043 : END IF
4044 888525 : DO jatom = 1, natom
4045 888525 : IF (jatom /= iatom) THEN
4046 296175 : IF (jatom < iatom) THEN
4047 148094 : IF (.NOT. skip_me(jatom)) THEN
4048 420184 : dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom)
4049 420184 : d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)
4050 105046 : IF (is_constraint(iatom)) THEN
4051 : d_sum_const_dR(:, jatom) = d_sum_const_dR(:, jatom) + &
4052 : dP_i_dRj(:, iatom, jatom)* &
4053 420184 : coefficients(iatom)
4054 : END IF
4055 : CYCLE
4056 : END IF
4057 : END IF
4058 764516 : dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom)
4059 764516 : d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)
4060 191129 : IF (is_constraint(iatom)) THEN
4061 : d_sum_const_dR(:, jatom) = d_sum_const_dR(:, jatom) + dP_i_dRj(:, iatom, jatom)* &
4062 764516 : coefficients(iatom)
4063 : END IF
4064 : END IF
4065 : END DO
4066 : END IF
4067 : ELSE
4068 1402232 : cell_functions(iatom) = 0.0_dp
4069 1402232 : skip_me(iatom) = .TRUE.
4070 1402232 : IF (cdft_control%becke_control%should_skip) THEN
4071 858928 : IF (is_constraint(iatom)) nskipped = nskipped + 1
4072 858928 : IF (nskipped == cdft_control%natoms) THEN
4073 410635 : IF (in_memory) THEN
4074 252800 : IF (cdft_control%becke_control%cavity_confine) THEN
4075 252800 : cavity(k, j, i) = 0.0_dp
4076 : END IF
4077 : END IF
4078 : EXIT
4079 : END IF
4080 : END IF
4081 : END IF
4082 : END DO
4083 988401 : IF (nskipped == cdft_control%natoms) CYCLE
4084 : sum_cell_f_constr = 0.0_dp
4085 1733298 : DO ip = 1, cdft_control%natoms
4086 : sum_cell_f_constr = sum_cell_f_constr + cell_functions(catom(ip))* &
4087 1733298 : cdft_control%group(1)%coeff(ip)
4088 : END DO
4089 577766 : sum_cell_f_all = 0.0_dp
4090 577766 : nwork = nwork + 1
4091 1733298 : DO ip = 1, natom
4092 1733298 : sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
4093 : END DO
4094 577766 : IF (in_memory) THEN
4095 1266948 : DO iatom = 1, natom
4096 1266948 : IF (ABS(sum_cell_f_all) > 0.0_dp) THEN
4097 : gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = &
4098 : d_sum_const_dR(:, iatom)/sum_cell_f_all - sum_cell_f_constr* &
4099 1529032 : d_sum_Pm_dR(:, iatom)/(sum_cell_f_all**2)
4100 : END IF
4101 : END DO
4102 : END IF
4103 616975 : IF (ABS(sum_cell_f_all) > eps_sum_cell_f_all) THEN
4104 320062 : weight(k, j, i) = sum_cell_f_constr/sum_cell_f_all
4105 : END IF
4106 : END DO ! i
4107 : END DO ! j
4108 : END DO ! k
4109 : ! Load balancing: post send requests
4110 80 : IF (iwork == 2) THEN
4111 4 : IF (.NOT. mixed_cdft%is_special) THEN
4112 12 : DO i = 1, SIZE(req_send, 1)
4113 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%cavity, &
4114 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4115 : request=req_send(i, 1), &
4116 8 : tag=mixed_cdft%dlb_control%dest_tags_repl(i))
4117 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%weight, &
4118 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4119 : request=req_send(i, 2), &
4120 8 : tag=mixed_cdft%dlb_control%dest_tags_repl(i) + 1)
4121 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%gradients, &
4122 : dest=mixed_cdft%dlb_control%my_dest_repl(i), &
4123 : request=req_send(i, 3), &
4124 12 : tag=mixed_cdft%dlb_control%dest_tags_repl(i) + 2)
4125 : END DO
4126 : should_communicate = .TRUE.
4127 : nsent_total = 0
4128 : ELSE
4129 0 : DO i = 1, SIZE(req_send, 1)
4130 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%cavity, &
4131 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4132 0 : request=req_send(i, 3*(ispecial - 1) + 1), tag=1)
4133 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%weight, &
4134 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4135 0 : request=req_send(i, 3*(ispecial - 1) + 2), tag=2)
4136 : CALL force_env%para_env%isend(msgin=mixed_cdft%dlb_control%sendbuff(ispecial)%gradients, &
4137 : dest=mixed_cdft%dlb_control%sendbuff(ispecial)%rank(i), &
4138 0 : request=req_send(i, 3*(ispecial - 1) + 3), tag=3)
4139 : END DO
4140 0 : IF (ispecial == my_special_work) THEN
4141 0 : should_communicate = .TRUE.
4142 0 : nsent_total = 0
4143 : END IF
4144 : END IF
4145 4 : work(mixed_cdft%dlb_control%my_source + 1) = work(mixed_cdft%dlb_control%my_source + 1) + nwork
4146 4 : work_dlb(force_env%para_env%mepos + 1) = work_dlb(force_env%para_env%mepos + 1) + nwork
4147 : ELSE
4148 36 : IF (mixed_cdft%dlb) work(force_env%para_env%mepos + 1) = work(force_env%para_env%mepos + 1) + nwork
4149 36 : IF (mixed_cdft%dlb) work_dlb(force_env%para_env%mepos + 1) = work_dlb(force_env%para_env%mepos + 1) + nwork
4150 : END IF
4151 : END DO ! ispecial
4152 : END DO ! iwork
4153 : ! Load balancing: wait for communication and deallocate sending buffers
4154 36 : IF (mixed_cdft%dlb) THEN
4155 16 : IF (mixed_cdft%dlb_control%recv_work .AND. &
4156 : ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
4157 48 : ALLOCATE (req_total(SIZE(req_recv) + SIZE(req_send, 1)*SIZE(req_send, 2)))
4158 4 : index = SIZE(req_recv)
4159 28 : req_total(1:index) = req_recv
4160 16 : DO i = 1, SIZE(req_send, 2)
4161 40 : DO j = 1, SIZE(req_send, 1)
4162 24 : index = index + 1
4163 36 : req_total(index) = req_send(j, i)
4164 : END DO
4165 : END DO
4166 4 : CALL mp_waitall(req_total)
4167 4 : DEALLOCATE (req_total)
4168 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4169 0 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
4170 : END IF
4171 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4172 0 : DEALLOCATE (mixed_cdft%dlb_control%weight)
4173 : END IF
4174 4 : IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4175 0 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
4176 : END IF
4177 4 : IF (mixed_cdft%is_special) THEN
4178 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4179 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4180 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4181 : END IF
4182 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4183 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4184 : END IF
4185 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4186 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4187 : END IF
4188 : END DO
4189 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4190 : END IF
4191 4 : DEALLOCATE (req_send, req_recv)
4192 4 : ELSE IF (mixed_cdft%dlb_control%recv_work) THEN
4193 0 : IF (should_communicate) THEN
4194 0 : CALL mp_waitall(req_send)
4195 : END IF
4196 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4197 0 : DEALLOCATE (mixed_cdft%dlb_control%cavity)
4198 : END IF
4199 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4200 0 : DEALLOCATE (mixed_cdft%dlb_control%weight)
4201 : END IF
4202 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4203 0 : DEALLOCATE (mixed_cdft%dlb_control%gradients)
4204 : END IF
4205 0 : IF (mixed_cdft%is_special) THEN
4206 0 : DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4207 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4208 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4209 : END IF
4210 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4211 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4212 : END IF
4213 0 : IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4214 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4215 : END IF
4216 : END DO
4217 0 : DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4218 : END IF
4219 0 : DEALLOCATE (req_send)
4220 8 : ELSE IF (ANY(mixed_cdft%dlb_control%recv_work_repl)) THEN
4221 4 : CALL mp_waitall(req_recv)
4222 4 : DEALLOCATE (req_recv)
4223 : END IF
4224 : END IF
4225 36 : IF (mixed_cdft%dlb) THEN
4226 40 : CALL force_env%para_env%sum(work)
4227 40 : CALL force_env%para_env%sum(work_dlb)
4228 8 : IF (.NOT. ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
4229 12 : ALLOCATE (mixed_cdft%dlb_control%prediction_error(force_env%para_env%num_pe))
4230 : END IF
4231 40 : mixed_cdft%dlb_control%prediction_error = mixed_cdft%dlb_control%expected_work - work
4232 : IF (debug_this_module .AND. iounit > 0) THEN
4233 : DO i = 1, SIZE(work, 1)
4234 : WRITE (iounit, '(A,I10,I10,I10)') &
4235 : 'Work', work(i), work_dlb(i), mixed_cdft%dlb_control%expected_work(i)
4236 : END DO
4237 : END IF
4238 8 : DEALLOCATE (work, work_dlb, mixed_cdft%dlb_control%expected_work)
4239 : END IF
4240 36 : NULLIFY (gradients, weight, cavity)
4241 36 : IF (ALLOCATED(coefficients)) THEN
4242 36 : DEALLOCATE (coefficients)
4243 : END IF
4244 36 : IF (in_memory) THEN
4245 24 : DEALLOCATE (ds_dR_j)
4246 24 : DEALLOCATE (ds_dR_i)
4247 24 : DEALLOCATE (d_sum_Pm_dR)
4248 24 : DEALLOCATE (d_sum_const_dR)
4249 24 : DEALLOCATE (dP_i_dRj)
4250 24 : DEALLOCATE (dP_i_dRi)
4251 24 : NULLIFY (gradients)
4252 24 : IF (store_vectors) THEN
4253 24 : DEALLOCATE (pair_dist_vecs)
4254 : END IF
4255 : END IF
4256 36 : NULLIFY (cutoffs)
4257 36 : IF (ALLOCATED(is_constraint)) THEN
4258 36 : DEALLOCATE (is_constraint)
4259 : END IF
4260 36 : DEALLOCATE (catom)
4261 36 : DEALLOCATE (R12)
4262 36 : DEALLOCATE (cell_functions)
4263 36 : DEALLOCATE (skip_me)
4264 36 : IF (ALLOCATED(completed)) THEN
4265 4 : DEALLOCATE (completed)
4266 : END IF
4267 36 : IF (ASSOCIATED(nsent)) THEN
4268 4 : DEALLOCATE (nsent)
4269 : END IF
4270 36 : IF (store_vectors) THEN
4271 36 : DEALLOCATE (distances)
4272 36 : DEALLOCATE (distance_vecs)
4273 36 : DEALLOCATE (position_vecs)
4274 : END IF
4275 36 : IF (ASSOCIATED(req_send)) THEN
4276 0 : DEALLOCATE (req_send)
4277 : END IF
4278 36 : IF (ASSOCIATED(req_recv)) THEN
4279 0 : DEALLOCATE (req_recv)
4280 : END IF
4281 : CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
4282 36 : "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
4283 36 : CALL timestop(handle)
4284 :
4285 72 : END SUBROUTINE mixed_becke_constraint_low
4286 :
4287 0 : END MODULE mixed_cdft_methods
|
