(git:f2099e5)
Loading...
Searching...
No Matches
mixed_cdft_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods for mixed CDFT calculations
10!> \par History
11!> Separated CDFT routines from mixed_environment_utils
12!> \author Nico Holmberg [01.2017]
13! **************************************************************************************************
14MODULE 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
144CONTAINS
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 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 NULLIFY (subsys_mix, force_env_qs, force_env_section, print_section, &
174 root_section, mixed_section, md_section, mixed_env, mixed_cdft, &
175 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 is_qmmm = .false.
183 logger => cp_get_default_logger()
184 cpassert(ASSOCIATED(force_env))
185 nforce_eval = SIZE(force_env%sub_force_env)
186 CALL timeset(routinen, handle)
187 CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
188 mixed_env => force_env%mixed_env
189 mixed_section => section_vals_get_subs_vals(force_env_section, "MIXED")
190 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
191 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
192 ! Check if a mixed CDFT calculation is requested
193 CALL section_vals_val_get(mixed_section, "MIXING_TYPE", i_val=mixing_type)
194 IF (mixing_type == mix_cdft .AND. .NOT. ASSOCIATED(mixed_env%cdft_control)) THEN
195 mixed_env%do_mixed_cdft = .true.
196 IF (mixed_env%do_mixed_cdft) THEN
197 ! Sanity check
198 IF (nforce_eval < 2) THEN
199 CALL cp_abort(__location__, &
200 "Mixed CDFT calculation requires at least 2 force_evals.")
201 END IF
202 mapping_section => section_vals_get_subs_vals(mixed_section, "MAPPING")
203 CALL section_vals_get(mapping_section, explicit=explicit)
204 ! The sub_force_envs must share the same geometrical structure
205 IF (explicit) THEN
206 cpabort("Please disable section &MAPPING for mixed CDFT calculations")
207 END IF
208 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%COUPLING", i_val=et_freq)
209 IF (et_freq < 0) THEN
210 mixed_env%do_mixed_et = .false.
211 ELSE
212 mixed_env%do_mixed_et = .true.
213 IF (et_freq == 0) THEN
214 mixed_env%et_freq = 1
215 ELSE
216 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 DO iforce_eval = 1, nforce_eval
222 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
223 SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
224 CASE (use_qs_force)
225 force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
226 CASE (use_qmmm)
227 is_qmmm = .true.
228 ! This is really the container for QMMM
229 force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
230 CASE (use_qmmmx)
231 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 "supported for mixed_cdft_init")
236 END SELECT
237 cpassert(ASSOCIATED(force_env_qs))
238 END DO
239 ! Get infos about the mixed subsys
240 IF (.NOT. is_qmmm) THEN
241 CALL force_env_get(force_env=force_env, &
242 subsys=subsys_mix)
243 CALL cp_subsys_get(subsys=subsys_mix, &
244 particles=particles_mix)
245 ELSE
246 CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
247 cp_subsys=subsys_mix)
248 CALL cp_subsys_get(subsys=subsys_mix, &
249 particles=particles_mix)
250 END IF
251 ! Init mixed_cdft_type
252 ALLOCATE (mixed_cdft)
253 CALL mixed_cdft_type_create(mixed_cdft)
254 mixed_cdft%first_iteration = .true.
255 ! Determine what run type to use
256 IF (mixed_env%ngroups == 1) THEN
257 ! States treated in serial, possibly copying CDFT weight function and gradients from state to state
258 mixed_cdft%run_type = mixed_cdft_serial
259 ELSE IF (mixed_env%ngroups == 2) THEN
260 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%PARALLEL_BUILD", l_val=is_parallel)
261 IF (is_parallel) THEN
262 ! Treat states in parallel, build weight function and gradients in parallel before SCF process
263 mixed_cdft%run_type = mixed_cdft_parallel
264 IF (.NOT. nforce_eval == 2) THEN
265 CALL cp_abort(__location__, &
266 "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 mixed_cdft%run_type = mixed_cdft_parallel_nobuild
271 END IF
272 ELSE
273 mixed_cdft%run_type = mixed_cdft_parallel_nobuild
274 END IF
275 ! Store QMMM flag
276 mixed_env%do_mixed_qmmm_cdft = is_qmmm
277 ! Setup dynamic load balancing
278 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%DLB", l_val=mixed_cdft%dlb)
279 mixed_cdft%dlb = mixed_cdft%dlb .AND. calculate_forces ! disable if forces are not needed
280 mixed_cdft%dlb = mixed_cdft%dlb .AND. (mixed_cdft%run_type == mixed_cdft_parallel) ! disable if not parallel
281 IF (mixed_cdft%dlb) THEN
282 ALLOCATE (mixed_cdft%dlb_control)
283 NULLIFY (mixed_cdft%dlb_control%weight, mixed_cdft%dlb_control%gradients, &
284 mixed_cdft%dlb_control%cavity, mixed_cdft%dlb_control%target_list, &
285 mixed_cdft%dlb_control%bo, mixed_cdft%dlb_control%expected_work, &
286 mixed_cdft%dlb_control%prediction_error, mixed_cdft%dlb_control%sendbuff, &
287 mixed_cdft%dlb_control%recvbuff, mixed_cdft%dlb_control%recv_work_repl, &
288 mixed_cdft%dlb_control%recv_info)
289 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOAD_SCALE", &
290 r_val=mixed_cdft%dlb_control%load_scale)
291 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%VERY_OVERLOADED", &
292 r_val=mixed_cdft%dlb_control%very_overloaded)
293 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%MORE_WORK", &
294 i_val=mixed_cdft%dlb_control%more_work)
295 END IF
296 ! Metric/Wavefunction overlap method/Lowdin orthogonalization/CDFT-CI
297 mixed_cdft%calculate_metric = .false.
298 mixed_cdft%wfn_overlap_method = .false.
299 mixed_cdft%use_lowdin = .false.
300 mixed_cdft%do_ci = .false.
301 mixed_cdft%nonortho_coupling = .false.
302 mixed_cdft%identical_constraints = .true.
303 mixed_cdft%block_diagonalize = .false.
304 IF (mixed_env%do_mixed_et) THEN
305 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%METRIC", &
306 l_val=mixed_cdft%calculate_metric)
307 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%WFN_OVERLAP", &
308 l_val=mixed_cdft%wfn_overlap_method)
309 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%LOWDIN", &
310 l_val=mixed_cdft%use_lowdin)
311 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%CI", &
312 l_val=mixed_cdft%do_ci)
313 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%NONORTHOGONAL_COUPLING", &
314 l_val=mixed_cdft%nonortho_coupling)
315 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%BLOCK_DIAGONALIZE", &
316 l_val=mixed_cdft%block_diagonalize)
317 END IF
318 ! Inversion method
319 CALL section_vals_val_get(mixed_section, "MIXED_CDFT%EPS_SVD", r_val=mixed_cdft%eps_svd)
320 IF (mixed_cdft%eps_svd < 0.0_dp .OR. mixed_cdft%eps_svd > 1.0_dp) THEN
321 cpabort("Illegal value for EPS_SVD. Value must be between 0.0 and 1.0.")
322 END IF
323 ! MD related settings
324 CALL force_env_get(force_env, root_section=root_section)
325 md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
326 CALL section_vals_val_get(md_section, "TIMESTEP", r_val=mixed_cdft%sim_dt)
327 CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=mixed_cdft%sim_step)
328 mixed_cdft%sim_step = mixed_cdft%sim_step - 1 ! to get the first step correct
329 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 settings, natom=SIZE(particles_mix%els))
333 ! Transfer settings to mixed_cdft
334 CALL mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
335 ! Initilize necessary structures
336 CALL mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
337 ! Write information about the mixed CDFT calculation
338 IF (iounit > 0) THEN
339 WRITE (iounit, *) ""
340 WRITE (iounit, fmt="(T2,A,T71)") &
341 "MIXED_CDFT| Activating mixed CDFT calculation"
342 WRITE (iounit, fmt="(T2,A,T71,I10)") &
343 "MIXED_CDFT| Number of CDFT states: ", nforce_eval
344 SELECT CASE (mixed_cdft%run_type)
345 CASE (mixed_cdft_parallel)
346 WRITE (iounit, fmt="(T2,A,T71)") &
347 "MIXED_CDFT| CDFT states calculation mode: parallel with build"
348 WRITE (iounit, fmt="(T2,A,T71)") &
349 "MIXED_CDFT| Becke constraint is first built using all available processors"
350 WRITE (iounit, fmt="(T2,A,T71)") &
351 " and then copied to both states with their own processor groups"
352 CASE (mixed_cdft_serial)
353 WRITE (iounit, fmt="(T2,A,T71)") &
354 "MIXED_CDFT| CDFT states calculation mode: serial"
355 IF (mixed_cdft%identical_constraints) THEN
356 WRITE (iounit, fmt="(T2,A,T71)") &
357 "MIXED_CDFT| The constraints are built before the SCF procedure of the first"
358 WRITE (iounit, fmt="(T2,A,T71)") &
359 " CDFT state and subsequently copied to the other states"
360 ELSE
361 WRITE (iounit, fmt="(T2,A,T71)") &
362 "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 "MIXED_CDFT| CDFT states calculation mode: parallel without build"
367 WRITE (iounit, fmt="(T2,A,T71)") &
368 "MIXED_CDFT| The constraints are separately built for all CDFT states"
369 CASE DEFAULT
370 cpabort("Unknown mixed CDFT run type.")
371 END SELECT
372 WRITE (iounit, fmt="(T2,A,T71,L10)") &
373 "MIXED_CDFT| Calculating electronic coupling between states: ", mixed_env%do_mixed_et
374 WRITE (iounit, fmt="(T2,A,T71,L10)") &
375 "MIXED_CDFT| Calculating electronic coupling reliability metric: ", mixed_cdft%calculate_metric
376 WRITE (iounit, fmt="(T2,A,T71,L10)") &
377 "MIXED_CDFT| Configuration interaction (CDFT-CI) was requested: ", mixed_cdft%do_ci
378 WRITE (iounit, fmt="(T2,A,T71,L10)") &
379 "MIXED_CDFT| Block diagonalizing the mixed CDFT Hamiltonian: ", mixed_cdft%block_diagonalize
380 IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
381 WRITE (iounit, fmt="(T2,A,T71,L10)") &
382 "MIXED_CDFT| Dynamic load balancing enabled: ", mixed_cdft%dlb
383 IF (mixed_cdft%dlb) THEN
384 WRITE (iounit, fmt="(T2,A,T71)") "MIXED_CDFT| Dynamic load balancing parameters:"
385 WRITE (iounit, fmt="(T2,A,T71,F10.2)") &
386 "MIXED_CDFT| load_scale:", mixed_cdft%dlb_control%load_scale
387 WRITE (iounit, fmt="(T2,A,T71,F10.2)") &
388 "MIXED_CDFT| very_overloaded:", mixed_cdft%dlb_control%very_overloaded
389 WRITE (iounit, fmt="(T2,A,T71,I10)") &
390 "MIXED_CDFT| more_work:", mixed_cdft%dlb_control%more_work
391 END IF
392 END IF
393 IF (mixed_env%do_mixed_et) THEN
394 IF (mixed_cdft%eps_svd == 0.0_dp) THEN
395 WRITE (iounit, fmt="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with LU decomposition."
396 ELSE
397 WRITE (iounit, fmt="(T2,A,T71)") "MIXED_CDFT| Matrix inversions calculated with SVD decomposition."
398 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 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 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
407 CALL timestop(handle)
408
409 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 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 NULLIFY (mixed_cdft)
427 cpassert(ASSOCIATED(force_env))
428 CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
429 cpassert(ASSOCIATED(mixed_cdft))
430 IF (.NOT. PRESENT(iforce_eval)) THEN
431 SELECT CASE (mixed_cdft%run_type)
432 CASE (mixed_cdft_parallel)
433 CALL mixed_cdft_build_weight_parallel(force_env, calculate_forces)
434 CASE (mixed_cdft_parallel_nobuild)
435 CALL mixed_cdft_set_flags(force_env)
436 CASE DEFAULT
437 ! Do nothing
438 END SELECT
439 ELSE
440 SELECT CASE (mixed_cdft%run_type)
441 CASE (mixed_cdft_serial)
442 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 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 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 INTEGER, DIMENSION(:), POINTER :: lb, sendbuffer_i, ub
467 REAL(kind=dp) :: t1, t2
468 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 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 NULLIFY (subsys_mix, force_env_qs, particles_mix, force_env_section, print_section, &
484 mixed_env, mixed_cdft, pw_env, auxbas_pw_pool, mixed_auxbas_pw_pool, &
485 qs_env, dft_control, sendbuffer_i, lb, ub, req_total, recvbuffer, &
486 cdft_control, cdft_control_target)
487
488 logger => cp_get_default_logger()
489 cpassert(ASSOCIATED(force_env))
490 nforce_eval = SIZE(force_env%sub_force_env)
491 CALL timeset(routinen, handle)
492 t1 = m_walltime()
493 ! Get infos about the mixed subsys
494 CALL force_env_get(force_env=force_env, &
495 subsys=subsys_mix, &
496 force_env_section=force_env_section)
497 CALL cp_subsys_get(subsys=subsys_mix, &
498 particles=particles_mix)
499 DO iforce_eval = 1, nforce_eval
500 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
501 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 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 "supported for mixed_cdft_build_weight_parallel")
510 END SELECT
511 END DO
512 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 subsys=subsys_mix)
516 CALL cp_subsys_get(subsys=subsys_mix, &
517 particles=particles_mix)
518 ELSE
519 qs_env => force_env_qs%qmmm_env%qs_env
520 CALL get_qs_env(qs_env, cp_subsys=subsys_mix)
521 CALL cp_subsys_get(subsys=subsys_mix, &
522 particles=particles_mix)
523 END IF
524 mixed_env => force_env%mixed_env
525 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
526 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
527 CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
528 cpassert(ASSOCIATED(mixed_cdft))
529 cdft_control => mixed_cdft%cdft_control
530 cpassert(ASSOCIATED(cdft_control))
531 ! Calculate the Becke weight function and gradient on all np processors
532 CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=mixed_auxbas_pw_pool)
533 natom = SIZE(particles_mix%els)
534 CALL mixed_becke_constraint(force_env, calculate_forces)
535 ! Start replicating the working arrays on both np/2 processor groups
536 mixed_cdft%sim_step = mixed_cdft%sim_step + 1
537 CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
538 cdft_control_target => dft_control%qs_control%cdft_control
539 cpassert(dft_control%qs_control%cdft)
540 cpassert(ASSOCIATED(cdft_control_target))
541 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
542 bo = auxbas_pw_pool%pw_grid%bounds_local
543 ! First determine the size of the arrays along the confinement dir
544 IF (mixed_cdft%is_special) THEN
545 my_special_work = 2
546 ELSE
547 my_special_work = 1
548 END IF
549 ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)))
550 ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)))
551 ALLOCATE (lb(SIZE(mixed_cdft%source_list)), ub(SIZE(mixed_cdft%source_list)))
552 IF (cdft_control%becke_control%cavity_confine) THEN
553 ! Gaussian confinement => the bounds depend on the processor and need to be communicated
554 ALLOCATE (sendbuffer_i(2))
555 sendbuffer_i = cdft_control%becke_control%confine_bounds
556 DO i = 1, SIZE(mixed_cdft%source_list)
557 ALLOCATE (recvbuffer(i)%iv(2))
558 CALL force_env%para_env%irecv(msgout=recvbuffer(i)%iv, source=mixed_cdft%source_list(i), &
559 request=req_total(i))
560 END DO
561 DO i = 1, my_special_work
562 DO j = 1, SIZE(mixed_cdft%dest_list)
563 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 request=req_total(ind))
567 END DO
568 END DO
569 CALL mp_waitall(req_total)
570 ! Find the largest/smallest value of ub/lb
571 DEALLOCATE (sendbuffer_i)
572 lb_min = huge(0)
573 ub_max = -huge(0)
574 DO i = 1, SIZE(mixed_cdft%source_list)
575 lb(i) = recvbuffer(i)%iv(1)
576 ub(i) = recvbuffer(i)%iv(2)
577 IF (lb(i) < lb_min) lb_min = lb(i)
578 IF (ub(i) > ub_max) ub_max = ub(i)
579 DEALLOCATE (recvbuffer(i)%iv)
580 END DO
581 ! Take into account the grids already communicated during dlb
582 IF (mixed_cdft%dlb) THEN
583 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
584 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
585 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
586 DO i = 1, SIZE(mixed_cdft%dlb_control%recvbuff(j)%buffs)
587 IF (lbound(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
588 < lb_min) lb_min = lbound(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3)
589 IF (ubound(mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight, 3) &
590 > 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 ub_max = bo(2, 3)
599 lb_min = bo(1, 3)
600 lb = lb_min
601 ub = ub_max
602 END IF
603 ! Determine the sender specific indices of grid slices that are to be received
604 CALL timeset(routinen//"_comm", handle2)
605 DO j = 1, SIZE(recvbuffer)
606 ind = j + (j/2)
607 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 lb(j), ub(j)]
611 ELSE IF (mixed_cdft%is_pencil) THEN
612 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 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 IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_special) THEN
618 IF (mixed_cdft%dlb_control%recv_work_repl(1) .OR. mixed_cdft%dlb_control%recv_work_repl(2)) THEN
619 DO j = 1, 2
620 recv_offset = 0
621 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
622 recv_offset = sum(mixed_cdft%dlb_control%recv_info(j)%target_list(2, :))
623 END IF
624 IF (mixed_cdft%is_pencil) THEN
625 recvbuffer(j)%imap(1) = recvbuffer(j)%imap(1) + recv_offset
626 ELSE
627 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 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 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 request=req_total(j))
641 END DO
642 DO i = 1, my_special_work
643 DO j = 1, SIZE(mixed_cdft%dest_list)
644 ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
645 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 request=req_total(ind))
649 ELSE
650 CALL force_env%para_env%isend(msgin=mixed_cdft%weight, dest=mixed_cdft%dest_list(j), &
651 request=req_total(ind))
652 END IF
653 END DO
654 END DO
655 CALL mp_waitall(req_total)
656 IF (mixed_cdft%is_special) THEN
657 DO j = 1, SIZE(mixed_cdft%dest_list)
658 DEALLOCATE (mixed_cdft%sendbuff(j)%weight)
659 END DO
660 ELSE
661 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 ALLOCATE (cdft_control_target%group(1)%weight)
666 CALL auxbas_pw_pool%create_pw(cdft_control_target%group(1)%weight)
667 CALL pw_zero(cdft_control_target%group(1)%weight)
668 ! Assemble the recved slices
669 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 recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
673 END DO
674 ! Do the same for slices sent during dlb
675 IF (mixed_cdft%dlb) THEN
676 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
677 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
678 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
679 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 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 mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%weight
690 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 IF (cdft_control%becke_control%cavity_confine) THEN
698 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 request=req_total(j))
701 END DO
702 DO i = 1, my_special_work
703 DO j = 1, SIZE(mixed_cdft%dest_list)
704 ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
705 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 request=req_total(ind))
709 ELSE
710 CALL force_env%para_env%isend(msgin=mixed_cdft%cavity, dest=mixed_cdft%dest_list(j), &
711 request=req_total(ind))
712 END IF
713 END DO
714 END DO
715 CALL mp_waitall(req_total)
716 IF (mixed_cdft%is_special) THEN
717 DO j = 1, SIZE(mixed_cdft%dest_list)
718 DEALLOCATE (mixed_cdft%sendbuff(j)%cavity)
719 END DO
720 ELSE
721 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 lb_min:ub_max))
727 cdft_control_target%becke_control%cavity_mat = 0.0_dp
728
729 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 recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r3
733 END DO
734 IF (mixed_cdft%dlb) THEN
735 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
736 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
737 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
738 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 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 mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%cavity
749 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 DO j = 1, SIZE(mixed_cdft%source_list)
757 DEALLOCATE (recvbuffer(j)%r3)
758 END DO
759 IF (calculate_forces) THEN
760 ! Gradients of the weight function
761 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 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 request=req_total(j))
767 END DO
768 DO i = 1, my_special_work
769 DO j = 1, SIZE(mixed_cdft%dest_list)
770 ind = j + (i - 1)*SIZE(mixed_cdft%dest_list) + SIZE(mixed_cdft%source_list)
771 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 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 request=req_total(ind))
778 END IF
779 END DO
780 END DO
781 CALL mp_waitall(req_total)
782 IF (mixed_cdft%is_special) THEN
783 DO j = 1, SIZE(mixed_cdft%dest_list)
784 DEALLOCATE (mixed_cdft%sendbuff(j)%gradients)
785 END DO
786 DEALLOCATE (mixed_cdft%sendbuff)
787 ELSE
788 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 bo(1, 2):bo(2, 2), lb_min:ub_max))
792 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 recvbuffer(j)%imap(5):recvbuffer(j)%imap(6)) = recvbuffer(j)%r4
796 DEALLOCATE (recvbuffer(j)%r4)
797 END DO
798 IF (mixed_cdft%dlb) THEN
799 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
800 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
801 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
802 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 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 mixed_cdft%dlb_control%recvbuff(j)%buffs(i)%gradients
813 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 IF (mixed_cdft%dlb) THEN
822 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
823 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
824 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
825 IF (ASSOCIATED(mixed_cdft%dlb_control%recv_info(j)%target_list)) THEN
826 DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%target_list)
827 END IF
828 DEALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs)
829 END IF
830 END DO
831 DEALLOCATE (mixed_cdft%dlb_control%recv_info, mixed_cdft%dlb_control%recvbuff)
832 END IF
833 IF (ASSOCIATED(mixed_cdft%dlb_control%target_list)) THEN
834 DEALLOCATE (mixed_cdft%dlb_control%target_list)
835 END IF
836 DEALLOCATE (mixed_cdft%dlb_control%recv_work_repl)
837 END IF
838 DEALLOCATE (recvbuffer)
839 DEALLOCATE (req_total)
840 DEALLOCATE (lb)
841 DEALLOCATE (ub)
842 CALL timestop(handle2)
843 ! Set some flags so the weight is not rebuilt during SCF
844 cdft_control_target%external_control = .true.
845 cdft_control_target%need_pot = .false.
846 cdft_control_target%transfer_pot = .false.
847 ! Store the bound indices for force calculation
848 IF (calculate_forces) THEN
849 cdft_control_target%becke_control%confine_bounds(2) = ub_max
850 cdft_control_target%becke_control%confine_bounds(1) = lb_min
851 END IF
852 CALL pw_scale(cdft_control_target%group(1)%weight, &
853 cdft_control_target%group(1)%weight%pw_grid%dvol)
854 ! Set flags for ET coupling calculation
855 IF (mixed_env%do_mixed_et) THEN
856 IF (modulo(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
857 dft_control%qs_control%cdft_control%do_et = .true.
858 dft_control%qs_control%cdft_control%calculate_metric = mixed_cdft%calculate_metric
859 ELSE
860 dft_control%qs_control%cdft_control%do_et = .false.
861 dft_control%qs_control%cdft_control%calculate_metric = .false.
862 END IF
863 END IF
864 t2 = m_walltime()
865 IF (iounit > 0) THEN
866 WRITE (iounit, '(A)') ' '
867 WRITE (iounit, '(T2,A,F6.1,A)') 'MIXED_CDFT| Becke constraint built in ', t2 - t1, ' seconds'
868 WRITE (iounit, '(A)') ' '
869 END IF
870 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
871 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
872 CALL timestop(handle)
873
874 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 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 NULLIFY (mixed_cdft, dft_control_source, dft_control_target, force_env_qs_source, &
905 force_env_qs_target, pw_env_source, pw_env_target, auxbas_pw_pool_source, &
906 auxbas_pw_pool_target, qs_env_source, qs_env_target, mixed_env, &
907 cdft_control_source, cdft_control_target)
908 mixed_env => force_env%mixed_env
909 CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
910 CALL timeset(routinen, handle)
911 IF (iforce_eval == 1) THEN
912 jforce_eval = 1
913 ELSE
914 jforce_eval = iforce_eval - 1
915 END IF
916 nforce_eval = SIZE(force_env%sub_force_env)
917 SELECT CASE (force_env%sub_force_env(jforce_eval)%force_env%in_use)
918 CASE (use_qs_force, use_qmmm)
919 force_env_qs_source => force_env%sub_force_env(jforce_eval)%force_env
920 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 "supported for mixed_cdft_transfer_weight")
925 END SELECT
926 IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
927 CALL force_env_get(force_env=force_env_qs_source, &
928 qs_env=qs_env_source)
929 CALL force_env_get(force_env=force_env_qs_target, &
930 qs_env=qs_env_target)
931 ELSE
932 qs_env_source => force_env_qs_source%qmmm_env%qs_env
933 qs_env_target => force_env_qs_target%qmmm_env%qs_env
934 END IF
935 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 CALL get_qs_env(qs_env_source, dft_control=dft_control_source)
939 cdft_control_source => dft_control_source%qs_control%cdft_control
940 cdft_control_source%external_control = .false.
941 cdft_control_source%need_pot = .true.
942 IF (mixed_cdft%identical_constraints) THEN
943 cdft_control_source%transfer_pot = .true.
944 ELSE
945 cdft_control_source%transfer_pot = .false.
946 END IF
947 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 pw_env=pw_env_source)
952 CALL pw_env_get(pw_env_source, auxbas_pw_pool=auxbas_pw_pool_source)
953 cdft_control_source => dft_control_source%qs_control%cdft_control
954 CALL get_qs_env(qs_env_target, dft_control=dft_control_target, &
955 pw_env=pw_env_target)
956 CALL pw_env_get(pw_env_target, auxbas_pw_pool=auxbas_pw_pool_target)
957 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 IF (mixed_cdft%identical_constraints) THEN
960 ! Weight function
961 DO igroup = 1, SIZE(cdft_control_target%group)
962 ALLOCATE (cdft_control_target%group(igroup)%weight)
963 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 CALL pw_copy(cdft_control_source%group(igroup)%weight, cdft_control_target%group(igroup)%weight)
966 CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%group(igroup)%weight)
967 DEALLOCATE (cdft_control_source%group(igroup)%weight)
968 END DO
969 ! Cavity
970 IF (cdft_control_source%type == outer_scf_becke_constraint) THEN
971 IF (cdft_control_source%becke_control%cavity_confine) THEN
972 CALL auxbas_pw_pool_target%create_pw(cdft_control_target%becke_control%cavity)
973 CALL pw_copy(cdft_control_source%becke_control%cavity, cdft_control_target%becke_control%cavity)
974 CALL auxbas_pw_pool_source%give_back_pw(cdft_control_source%becke_control%cavity)
975 END IF
976 END IF
977 ! Gradients
978 IF (calculate_forces) THEN
979 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 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 bounds_of(5):bounds_of(6), bounds_of(7):bounds_of(8)))
991 cdft_control_target%group(igroup)%gradients = cdft_control_source%group(igroup)%gradients
992 DEALLOCATE (cdft_control_source%group(igroup)%gradients)
993 END DO
994 END IF
995 ! Atomic weight functions needed for CDFT charges
996 IF (cdft_control_source%atomic_charges) THEN
997 IF (.NOT. ASSOCIATED(cdft_control_target%charge)) THEN
998 ALLOCATE (cdft_control_target%charge(cdft_control_target%natoms))
999 END IF
1000 DO iatom = 1, cdft_control_target%natoms
1001 CALL auxbas_pw_pool_target%create_pw(cdft_control_target%charge(iatom))
1002 CALL pw_copy(cdft_control_source%charge(iatom), cdft_control_target%charge(iatom))
1003 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 cdft_control_target%external_control = .false.
1008 cdft_control_target%need_pot = .false.
1009 ! For states i+1 < nforce_eval, prevent deallocation of constraint
1010 IF (iforce_eval == nforce_eval) THEN
1011 cdft_control_target%transfer_pot = .false.
1012 ELSE
1013 cdft_control_target%transfer_pot = .true.
1014 END IF
1015 cdft_control_target%first_iteration = .false.
1016 ELSE
1017 ! Force rebuild of constraint and dont save it
1018 cdft_control_target%external_control = .false.
1019 cdft_control_target%need_pot = .true.
1020 cdft_control_target%transfer_pot = .false.
1021 IF (first_call) THEN
1022 cdft_control_target%first_iteration = .true.
1023 ELSE
1024 cdft_control_target%first_iteration = .false.
1025 END IF
1026 END IF
1027 END IF
1028 ! Set flags for ET coupling calculation
1029 IF (mixed_env%do_mixed_et) THEN
1030 IF (modulo(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
1031 IF (iforce_eval == 1) THEN
1032 cdft_control_source%do_et = .true.
1033 cdft_control_source%calculate_metric = mixed_cdft%calculate_metric
1034 ELSE
1035 cdft_control_target%do_et = .true.
1036 cdft_control_target%calculate_metric = mixed_cdft%calculate_metric
1037 END IF
1038 ELSE
1039 IF (iforce_eval == 1) THEN
1040 cdft_control_source%do_et = .false.
1041 cdft_control_source%calculate_metric = .false.
1042 ELSE
1043 cdft_control_target%do_et = .false.
1044 cdft_control_target%calculate_metric = .false.
1045 END IF
1046 END IF
1047 END IF
1048 IF (iforce_eval == nforce_eval .AND. first_call) first_call = .false.
1049 CALL timestop(handle)
1050
1051 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 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 NULLIFY (mixed_cdft, dft_control, force_env_qs, qs_env, mixed_env, cdft_control)
1075 mixed_env => force_env%mixed_env
1076 CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1077 CALL timeset(routinen, handle)
1078 nforce_eval = SIZE(force_env%sub_force_env)
1079 DO iforce_eval = 1, nforce_eval
1080 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1081 SELECT CASE (force_env%sub_force_env(iforce_eval)%force_env%in_use)
1082 CASE (use_qs_force, use_qmmm)
1083 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 "supported for mixed_cdft_set_flags")
1088 END SELECT
1089 IF (.NOT. force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1090 CALL force_env_get(force_env=force_env_qs, qs_env=qs_env)
1091 ELSE
1092 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 CALL get_qs_env(qs_env, dft_control=dft_control)
1097 cdft_control => dft_control%qs_control%cdft_control
1098 cdft_control%external_control = .false.
1099 cdft_control%need_pot = .true.
1100 cdft_control%transfer_pot = .false.
1101 IF (first_call) THEN
1102 cdft_control%first_iteration = .true.
1103 ELSE
1104 cdft_control%first_iteration = .false.
1105 END IF
1106 ! Set flags for ET coupling calculation
1107 IF (mixed_env%do_mixed_et) THEN
1108 IF (modulo(mixed_cdft%sim_step, mixed_env%et_freq) == 0) THEN
1109 cdft_control%do_et = .true.
1110 cdft_control%calculate_metric = mixed_cdft%calculate_metric
1111 ELSE
1112 cdft_control%do_et = .false.
1113 cdft_control%calculate_metric = .false.
1114 END IF
1115 END IF
1116 END DO
1117 mixed_cdft%sim_step = mixed_cdft%sim_step + 1
1118 IF (first_call) first_call = .false.
1119 CALL timestop(handle)
1120
1121 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 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 cpassert(ASSOCIATED(force_env))
1137 CALL timeset(routinen, handle)
1138 ! Move needed arrays from individual CDFT states to the mixed CDFT env
1139 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 CALL mixed_cdft_interaction_matrices(force_env)
1144 ! Calculate eletronic couplings between states (Lowdin/rotation)
1145 CALL mixed_cdft_calculate_coupling_low(force_env)
1146 ! Print out couplings
1147 CALL mixed_cdft_print_couplings(force_env)
1148 ! Block diagonalize the mixed CDFT Hamiltonian matrix
1149 CALL mixed_cdft_block_diag(force_env)
1150 ! CDFT Configuration Interaction
1151 CALL mixed_cdft_configuration_interaction(force_env)
1152 ! Clean up
1153 CALL mixed_cdft_release_work(force_env)
1154 CALL timestop(handle)
1155
1156 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coupling_nonortho, eigenv, energy, sda
1180 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: h_mat, s_det, s_mat, strength, tmp_mat, &
1181 w_diagonal, wad, wda
1182 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: a, b
1183 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_overlap
1187 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: w_matrix_mo
1188 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: mixed_mo_coeff
1189 TYPE(cp_logger_type), POINTER :: logger
1190 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix, density_matrix_diff, &
1191 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 NULLIFY (force_env_section, print_section, mixed_cdft_section, &
1202 mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
1203 density_matrix_diff, mo_mo_fmstruct, &
1204 mixed_mo_coeff, mixed_matrix_s, &
1205 density_matrix, energy_qs, w_matrix, mo_eigval)
1206 logger => cp_get_default_logger()
1207 cpassert(ASSOCIATED(force_env))
1208 CALL timeset(routinen, handle)
1209 CALL force_env_get(force_env=force_env, &
1210 force_env_section=force_env_section)
1211 mixed_env => force_env%mixed_env
1212 nforce_eval = SIZE(force_env%sub_force_env)
1213 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1214 IF (section_get_lval(print_section, "MO_OVERLAP_MATRIX")) THEN
1215 print_mo = .true.
1216 mounit = cp_print_key_unit_nr(logger, print_section, extension='.moOverlap', on_file=.true.)
1217 ELSE
1218 print_mo = .false.
1219 END IF
1220 IF (section_get_lval(print_section, "MO_OVERLAP_EIGENVALUES")) THEN
1221 print_mo_eigval = .true.
1222 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 CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
1227 ! Get redistributed work matrices
1228 cpassert(ASSOCIATED(mixed_cdft))
1229 cpassert(ASSOCIATED(mixed_cdft%matrix%mixed_mo_coeff))
1230 cpassert(ASSOCIATED(mixed_cdft%matrix%w_matrix))
1231 cpassert(ASSOCIATED(mixed_cdft%matrix%mixed_matrix_s))
1232 mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
1233 w_matrix => mixed_cdft%matrix%w_matrix
1234 mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
1235 IF (mixed_cdft%calculate_metric) THEN
1236 cpassert(ASSOCIATED(mixed_cdft%matrix%density_matrix))
1237 density_matrix => mixed_cdft%matrix%density_matrix
1238 END IF
1239 ! Get number of weight functions per state
1240 nvar = SIZE(w_matrix, 2)
1241 nspins = SIZE(mixed_mo_coeff, 2)
1242 ! Check that the number of MOs/AOs is equal in every CDFT state
1243 ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
1244 DO ispin = 1, nspins
1245 CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=check_mo(1), nrow_global=check_ao(1))
1246 DO iforce_eval = 2, nforce_eval
1247 CALL cp_fm_get_info(mixed_mo_coeff(iforce_eval, ispin), ncol_global=check_mo(2), nrow_global=check_ao(2))
1248 IF (check_ao(1) /= check_ao(2)) THEN
1249 CALL cp_abort(__location__, &
1250 "The number of atomic orbitals must be the same in every CDFT state.")
1251 END IF
1252 IF (check_mo(1) /= check_mo(2)) THEN
1253 CALL cp_abort(__location__, &
1254 "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 npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
1260 ALLOCATE (w_matrix_mo(nforce_eval, nforce_eval, nvar))
1261 ALLOCATE (mo_overlap(npermutations), s_det(npermutations, nspins))
1262 ALLOCATE (a(nspins, nvar, npermutations), b(nspins, nvar, npermutations))
1263 a = 0.0_dp
1264 b = 0.0_dp
1265 IF (mixed_cdft%calculate_metric) THEN
1266 ALLOCATE (density_matrix_diff(npermutations, nspins))
1267 DO ispin = 1, nspins
1268 DO ipermutation = 1, npermutations
1269 NULLIFY (density_matrix_diff(ipermutation, ispin)%matrix)
1270 CALL dbcsr_init_p(density_matrix_diff(ipermutation, ispin)%matrix)
1271 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1272 CALL dbcsr_copy(density_matrix_diff(ipermutation, ispin)%matrix, &
1273 density_matrix(istate, ispin)%matrix, name="DENSITY_MATRIX")
1274 END DO
1275 END DO
1276 END IF
1277 ! Check for uniform occupations
1278 uniform_occupation = .NOT. ALLOCATED(mixed_cdft%occupations)
1279 should_scale = .false.
1280 IF (.NOT. uniform_occupation) THEN
1281 ALLOCATE (homo(nforce_eval, nspins))
1282 mixed_cdft_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT")
1283 CALL section_vals_val_get(mixed_cdft_section, "EPS_OCCUPIED", r_val=eps_occupied)
1284 IF (eps_occupied > 1.0_dp .OR. eps_occupied < 0.0_dp) THEN
1285 CALL cp_abort(__location__, &
1286 "Keyword EPS_OCCUPIED only accepts values between 0.0 and 1.0")
1287 END IF
1288 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 "orbitals is strongly recommended to screen nearly orthogonal states.")
1292 END IF
1293 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 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 NULLIFY (fm_struct_mo, mo_mo_fmstruct)
1300 CALL cp_fm_get_info(mixed_mo_coeff(1, ispin), ncol_global=ncol_mo(ispin), nrow_global=nrow_mo(ispin))
1301 nao = nrow_mo(ispin)
1302 IF (uniform_occupation) THEN
1303 nmo = ncol_mo(ispin)
1304 ELSE
1305 nmo = ncol_mo(ispin)
1306 ! Find indices of highest (fractionally) occupied molecular orbital
1307 homo(:, ispin) = nmo
1308 DO istate = 1, nforce_eval
1309 DO j = nmo, 1, -1
1310 IF (mixed_cdft%occupations(istate, ispin)%array(j) >= eps_occupied) THEN
1311 homo(istate, ispin) = j
1312 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 nmo = maxval(homo(:, ispin))
1320 ! Also check that the number of electrons is conserved (using a fixed sensible threshold)
1321 nelectron_mismatch = .false.
1322 nelectron_tot = sum(mixed_cdft%occupations(1, ispin)%array(1:nmo))
1323 DO istate = 2, nforce_eval
1324 IF (abs(sum(mixed_cdft%occupations(istate, ispin)%array(1:nmo)) - nelectron_tot) > 1.0e-4_dp) THEN
1325 nelectron_mismatch = .true.
1326 END IF
1327 END DO
1328 IF (any(homo(:, ispin) /= nmo)) THEN
1329 IF (ispin == 1) THEN
1330 CALL cp_warn(__location__, &
1331 "The number of occupied alpha MOs is not constant across all CDFT states. "// &
1332 "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 "Calculation proceeds but the results will likely be meaningless.")
1337 END IF
1338 ELSE IF (nelectron_mismatch) THEN
1339 IF (ispin == 1) THEN
1340 CALL cp_warn(__location__, &
1341 "The number of alpha electrons is not constant across all CDFT states. "// &
1342 "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 "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 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 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 name="ET_TMP_"//trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1357 CALL cp_fm_struct_release(fm_struct_mo)
1358 CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
1359 name="INVERSE_"//trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1360 CALL cp_fm_create(matrix=tinverse, matrix_struct=mo_mo_fmstruct, &
1361 name="T_INVERSE_"//trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1362 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 trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1366 END DO
1367 DO ivar = 1, nvar
1368 DO istate = 1, nforce_eval
1369 DO jstate = 1, nforce_eval
1370 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 trim(adjustl(cp_to_string(ivar)))//"_MATRIX")
1375 END DO
1376 END DO
1377 END DO
1378 CALL cp_fm_struct_release(mo_mo_fmstruct)
1379 ! Remove empty MOs and (possibly) scale rest with occupation numbers
1380 IF (.NOT. uniform_occupation) THEN
1381 DO iforce_eval = 1, nforce_eval
1382 CALL cp_fm_to_fm(mixed_mo_coeff(iforce_eval, ispin), tmp2, nmo, 1, 1)
1383 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 //trim(adjustl(cp_to_string(ispin)))//"_MATRIX")
1388 CALL cp_fm_to_fm(tmp2, mixed_mo_coeff(iforce_eval, ispin))
1389 IF (should_scale) THEN
1390 CALL cp_fm_column_scale(mixed_mo_coeff(iforce_eval, ispin), &
1391 mixed_cdft%occupations(iforce_eval, ispin)%array(1:nmo))
1392 END IF
1393 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 ipermutation = 0
1398 DO istate = 1, nforce_eval
1399 DO jstate = istate + 1, nforce_eval
1400 ipermutation = ipermutation + 1
1401 CALL cp_dbcsr_sm_fm_multiply(mixed_matrix_s, mixed_mo_coeff(istate, ispin), &
1402 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 tmp2, 0.0_dp, mo_overlap(ipermutation))
1406 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 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 DO ivar = 1, nvar
1417 DO jstate = 1, nforce_eval
1418 DO istate = 1, nforce_eval
1419 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 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 tmp2, 0.0_dp, w_matrix_mo(istate, jstate, ivar))
1427 END DO
1428 END DO
1429 END DO
1430 DO ipermutation = 1, npermutations
1431 ! Invert and calculate determinant of MO overlaps
1432 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1433 IF (print_mo_eigval) THEN
1434 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 eigval=mo_eigval)
1438 IF (moeigvalunit > 0) THEN
1439 IF (mixed_cdft%eps_svd == 0.0_dp) THEN
1440 WRITE (moeigvalunit, '(A,I2,A,I2,A,I1,A)') &
1441 "# MO Overlap matrix eigenvalues for CDFT states ", istate, " and ", jstate, &
1442 " (spin ", ispin, ")"
1443 ELSE
1444 WRITE (moeigvalunit, '(A,I2,A,I2,A,I1,A)') &
1445 "# MO Overlap matrix singular values for CDFT states ", istate, " and ", jstate, &
1446 " (spin ", ispin, ")"
1447 END IF
1448 WRITE (moeigvalunit, '(A1, A9, A12)') "#", "Index", adjustl("Value")
1449 DO j = 1, SIZE(mo_eigval)
1450 WRITE (moeigvalunit, '(I10, F12.8)') j, mo_eigval(j)
1451 END DO
1452 END IF
1453 DEALLOCATE (mo_eigval)
1454 ELSE
1455 CALL cp_fm_invert(mo_overlap(ipermutation), inverse_mat, &
1456 s_det(ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
1457 END IF
1458 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 DO j = 1, ncol_local
1461 DO k = 1, nrow_local
1462 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 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 CALL cp_fm_transpose(inverse_mat, tinverse)
1471 DO j = 1, ncol_local
1472 DO k = 1, nrow_local
1473 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 tinverse%local_data(k, j)
1477 END DO
1478 END DO
1479 END DO
1480 ! Handle different constraint types
1481 DO ivar = 1, nvar
1482 SELECT CASE (mixed_cdft%constraint_type(ivar, istate))
1483 CASE (cdft_charge_constraint)
1484 ! No action needed
1485 CASE (cdft_magnetization_constraint)
1486 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 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 IF (ispin == 1) a(ispin, ivar, ipermutation) = 0.0_dp
1493 CASE DEFAULT
1494 cpabort("Unknown constraint type.")
1495 END SELECT
1496 SELECT CASE (mixed_cdft%constraint_type(ivar, jstate))
1497 CASE (cdft_charge_constraint)
1498 ! No action needed
1499 CASE (cdft_magnetization_constraint)
1500 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 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 IF (ispin == 1) b(ispin, ivar, ipermutation) = 0.0_dp
1507 CASE DEFAULT
1508 cpabort("Unknown constraint type.")
1509 END SELECT
1510 END DO
1511 ! Compute density matrix difference P = P_j - P_i
1512 IF (mixed_cdft%calculate_metric) THEN
1513 CALL dbcsr_add(density_matrix_diff(ipermutation, ispin)%matrix, &
1514 density_matrix(jstate, ispin)%matrix, -1.0_dp, 1.0_dp)
1515 END IF
1516 !
1517 CALL force_env%para_env%sum(a(ispin, :, ipermutation))
1518 CALL force_env%para_env%sum(b(ispin, :, ipermutation))
1519 END DO
1520 ! Release work
1521 CALL cp_fm_release(tmp2)
1522 DO ivar = 1, nvar
1523 DO jstate = 1, nforce_eval
1524 DO istate = 1, nforce_eval
1525 IF (istate == jstate) cycle
1526 CALL cp_fm_release(w_matrix_mo(istate, jstate, ivar))
1527 END DO
1528 END DO
1529 END DO
1530 DO ipermutation = 1, npermutations
1531 CALL cp_fm_release(mo_overlap(ipermutation))
1532 END DO
1533 CALL cp_fm_release(tinverse)
1534 CALL cp_fm_release(inverse_mat)
1535 END DO
1536 DEALLOCATE (mo_overlap)
1537 DEALLOCATE (w_matrix_mo)
1538 IF (.NOT. uniform_occupation) THEN
1539 DEALLOCATE (homo)
1540 DEALLOCATE (mixed_cdft%occupations)
1541 END IF
1542 IF (print_mo) THEN
1543 CALL cp_print_key_finished_output(mounit, logger, force_env_section, &
1544 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO", on_file=.true.)
1545 END IF
1546 IF (print_mo_eigval) THEN
1547 CALL cp_print_key_finished_output(moeigvalunit, logger, force_env_section, &
1548 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO", on_file=.true.)
1549 END IF
1550 ! solve eigenstates for the projector matrix
1551 ALLOCATE (wda(nvar, npermutations))
1552 ALLOCATE (sda(npermutations))
1553 IF (.NOT. mixed_cdft%identical_constraints) ALLOCATE (wad(nvar, npermutations))
1554 DO ipermutation = 1, npermutations
1555 IF (nspins == 2) THEN
1556 sda(ipermutation) = abs(s_det(ipermutation, 1)*s_det(ipermutation, 2))
1557 ELSE
1558 sda(ipermutation) = s_det(ipermutation, 1)**2
1559 END IF
1560 ! Finalize <Psi_j | w_i(r) | Psi_i> by multiplication with Sda
1561 DO ivar = 1, nvar
1562 IF (mixed_cdft%identical_constraints) THEN
1563 wda(ivar, ipermutation) = (sum(a(:, ivar, ipermutation)) + sum(b(:, ivar, ipermutation)))* &
1564 sda(ipermutation)/2.0_dp
1565 ELSE
1566 wda(ivar, ipermutation) = sum(a(:, ivar, ipermutation))*sda(ipermutation)
1567 wad(ivar, ipermutation) = sum(b(:, ivar, ipermutation))*sda(ipermutation)
1568 END IF
1569 END DO
1570 END DO
1571 DEALLOCATE (a, b, s_det)
1572 ! Transfer info about the constraint calculations
1573 ALLOCATE (w_diagonal(nvar, nforce_eval), strength(nvar, nforce_eval), energy(nforce_eval))
1574 w_diagonal = 0.0_dp
1575 DO iforce_eval = 1, nforce_eval
1576 strength(:, iforce_eval) = mixed_env%strength(iforce_eval, :)
1577 END DO
1578 energy = 0.0_dp
1579 DO iforce_eval = 1, nforce_eval
1580 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
1581 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
1582 qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
1583 ELSE
1584 CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
1585 END IF
1586 CALL get_qs_env(qs_env, energy=energy_qs, dft_control=dft_control)
1587 IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
1588 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 energy(iforce_eval) = energy_qs%total - energy_qs%cdft
1593 END IF
1594 END DO
1595 CALL force_env%para_env%sum(w_diagonal)
1596 CALL force_env%para_env%sum(energy)
1597 CALL mixed_cdft_result_type_set(mixed_cdft%results, wda=wda, w_diagonal=w_diagonal, &
1598 energy=energy, strength=strength)
1599 IF (.NOT. mixed_cdft%identical_constraints) CALL mixed_cdft_result_type_set(mixed_cdft%results, wad=wad)
1600 ! Construct S
1601 ALLOCATE (s_mat(nforce_eval, nforce_eval))
1602 DO istate = 1, nforce_eval
1603 s_mat(istate, istate) = 1.0_dp
1604 END DO
1605 DO ipermutation = 1, npermutations
1606 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1607 s_mat(istate, jstate) = sda(ipermutation)
1608 s_mat(jstate, istate) = sda(ipermutation)
1609 END DO
1610 CALL mixed_cdft_result_type_set(mixed_cdft%results, s=s_mat)
1611 ! Invert S via eigendecomposition and compute S^-(1/2)
1612 ALLOCATE (eigenv(nforce_eval), tmp_mat(nforce_eval, nforce_eval))
1613 CALL diamat_all(s_mat, eigenv, .true.)
1614 tmp_mat = 0.0_dp
1615 DO istate = 1, nforce_eval
1616 IF (eigenv(istate) < 1.0e-14_dp) THEN
1617 ! Safeguard against division with 0 and negative numbers
1618 eigenv(istate) = 1.0e-14_dp
1619 CALL cp_warn(__location__, &
1620 "The overlap matrix is numerically nearly singular. "// &
1621 "Calculation proceeds but the results might be meaningless.")
1622 END IF
1623 tmp_mat(istate, istate) = 1.0_dp/sqrt(eigenv(istate))
1624 END DO
1625 tmp_mat(:, :) = matmul(tmp_mat, transpose(s_mat))
1626 s_mat(:, :) = matmul(s_mat, tmp_mat) ! S^(-1/2)
1627 CALL mixed_cdft_result_type_set(mixed_cdft%results, s_minushalf=s_mat)
1628 DEALLOCATE (eigenv, tmp_mat, s_mat)
1629 ! Construct nonorthogonal diabatic Hamiltonian matrix H''
1630 ALLOCATE (h_mat(nforce_eval, nforce_eval))
1631 IF (mixed_cdft%nonortho_coupling) ALLOCATE (coupling_nonortho(npermutations))
1632 DO istate = 1, nforce_eval
1633 h_mat(istate, istate) = energy(istate)
1634 END DO
1635 DO ipermutation = 1, npermutations
1636 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1637 sum_a = 0.0_dp
1638 sum_b = 0.0_dp
1639 DO ivar = 1, nvar
1640 ! V_J * <Psi_J | w_J(r) | Psi_J>
1641 sum_b(1) = sum_b(1) + strength(ivar, jstate)*w_diagonal(ivar, jstate)
1642 ! V_I * <Psi_I | w_I(r) | Psi_I>
1643 sum_a(1) = sum_a(1) + strength(ivar, istate)*w_diagonal(ivar, istate)
1644 IF (mixed_cdft%identical_constraints) THEN
1645 ! V_J * W_IJ
1646 sum_b(2) = sum_b(2) + strength(ivar, jstate)*wda(ivar, ipermutation)
1647 ! V_I * W_JI
1648 sum_a(2) = sum_a(2) + strength(ivar, istate)*wda(ivar, ipermutation)
1649 ELSE
1650 ! V_J * W_IJ
1651 sum_b(2) = sum_b(2) + strength(ivar, jstate)*wad(ivar, ipermutation)
1652 ! V_I * W_JI
1653 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 c(1) = (energy(jstate) + sum_b(1))*sda(ipermutation) - sum_b(2)
1659 ! H_JI = F_I*S_JI - V_I * W_JI
1660 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 h_mat(istate, jstate) = (c(1) + c(2))*0.5_dp
1663 h_mat(jstate, istate) = h_mat(istate, jstate)
1664 IF (mixed_cdft%nonortho_coupling) coupling_nonortho(ipermutation) = h_mat(istate, jstate)
1665 END DO
1666 CALL mixed_cdft_result_type_set(mixed_cdft%results, h=h_mat)
1667 DEALLOCATE (h_mat, w_diagonal, wda, strength, energy, sda)
1668 IF (ALLOCATED(wad)) DEALLOCATE (wad)
1669 IF (mixed_cdft%nonortho_coupling) THEN
1670 CALL mixed_cdft_result_type_set(mixed_cdft%results, nonortho=coupling_nonortho)
1671 DEALLOCATE (coupling_nonortho)
1672 END IF
1673 ! Compute metric to assess reliability of coupling
1674 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 IF (mixed_cdft%wfn_overlap_method) THEN
1678 IF (.NOT. uniform_occupation) THEN
1679 CALL cp_abort(__location__, &
1680 "Wavefunction overlap method supports only uniformly occupied MOs.")
1681 END IF
1682 CALL mixed_cdft_wfn_overlap_method(force_env, mixed_cdft, ncol_mo, nrow_mo)
1683 END IF
1684 ! Release remaining work
1685 DEALLOCATE (nrow_mo, ncol_mo)
1686 CALL mixed_cdft_work_type_release(mixed_cdft%matrix)
1687 CALL timestop(handle)
1688
1689 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coupling_lowdin, coupling_rotation, &
1706 eigenv
1707 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, w_mat
1708 TYPE(mixed_cdft_type), POINTER :: mixed_cdft
1709
1710 NULLIFY (mixed_cdft)
1711 cpassert(ASSOCIATED(force_env))
1712 CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1713 CALL timeset(routinen, handle)
1714 cpassert(ASSOCIATED(mixed_cdft))
1715 cpassert(ALLOCATED(mixed_cdft%results%W_diagonal))
1716 cpassert(ALLOCATED(mixed_cdft%results%Wda))
1717 cpassert(ALLOCATED(mixed_cdft%results%S_minushalf))
1718 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 nforce_eval = SIZE(mixed_cdft%results%H, 1)
1724 nvar = SIZE(mixed_cdft%results%Wda, 1)
1725 npermutations = nforce_eval*(nforce_eval - 1)/2
1726 ALLOCATE (tmp_mat(nforce_eval, nforce_eval))
1727 IF (nvar == 1 .AND. mixed_cdft%identical_constraints) THEN
1728 use_rotation = .true.
1729 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 ALLOCATE (w_mat(nforce_eval, nforce_eval), coupling_rotation(npermutations))
1738 ALLOCATE (eigenv(nforce_eval))
1739 ! W_mat(i, i) = N_i where N_i is the value of the constraint in state i
1740 DO istate = 1, nforce_eval
1741 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 DO ipermutation = 1, npermutations
1745 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1746 w_mat(istate, jstate) = sum(mixed_cdft%results%Wda(:, ipermutation))
1747 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 tmp_mat(:, :) = matmul(w_mat, mixed_cdft%results%S_minushalf) ! W * S^(-1/2)
1752 w_mat(:, :) = matmul(mixed_cdft%results%S_minushalf, tmp_mat) ! W' = S^(-1/2) * W * S^(-1/2)
1753 CALL diamat_all(w_mat, eigenv, .true.) ! Solve W'V' = AV'
1754 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 w_mat(:, :) = matmul(mixed_cdft%results%H, tmp_mat) ! H'' * V
1757 w_mat(:, :) = matmul(transpose(tmp_mat), w_mat) ! H = V^T * H'' * V
1758 DO ipermutation = 1, npermutations
1759 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1760 coupling_rotation(ipermutation) = w_mat(istate, jstate)
1761 END DO
1762 CALL mixed_cdft_result_type_set(mixed_cdft%results, rotation=coupling_rotation)
1763 DEALLOCATE (w_mat, coupling_rotation, eigenv)
1764 END IF
1765 ! Calculate coupling by Lowdin orthogonalization
1766 IF (use_lowdin) THEN
1767 ALLOCATE (coupling_lowdin(npermutations))
1768 tmp_mat(:, :) = matmul(mixed_cdft%results%H, mixed_cdft%results%S_minushalf) ! H'' * S^(-1/2)
1769 ! Final orthogonal diabatic Hamiltonian matrix H
1770 tmp_mat(:, :) = matmul(mixed_cdft%results%S_minushalf, tmp_mat) ! H = S^(-1/2) * H'' * S^(-1/2)
1771 DO ipermutation = 1, npermutations
1772 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
1773 coupling_lowdin(ipermutation) = tmp_mat(istate, jstate)
1774 END DO
1775 CALL mixed_cdft_result_type_set(mixed_cdft%results, lowdin=coupling_lowdin)
1776 DEALLOCATE (coupling_lowdin)
1777 END IF
1778 DEALLOCATE (tmp_mat)
1779 CALL timestop(handle)
1780
1781 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenv, work
1797 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 NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1806
1807 cpassert(ASSOCIATED(force_env))
1808 CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1809 cpassert(ASSOCIATED(mixed_cdft))
1810
1811 IF (.NOT. mixed_cdft%do_ci) RETURN
1812
1813 logger => cp_get_default_logger()
1814 CALL timeset(routinen, handle)
1815 CALL force_env_get(force_env=force_env, &
1816 force_env_section=force_env_section)
1817 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1818 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1819
1820 cpassert(ALLOCATED(mixed_cdft%results%S))
1821 cpassert(ALLOCATED(mixed_cdft%results%H))
1822 nforce_eval = SIZE(mixed_cdft%results%S, 1)
1823 ALLOCATE (s_mat(nforce_eval, nforce_eval), h_mat(nforce_eval, nforce_eval))
1824 ALLOCATE (eigenv(nforce_eval))
1825 s_mat(:, :) = mixed_cdft%results%S(:, :)
1826 h_mat(:, :) = mixed_cdft%results%H(:, :)
1827 ! Workspace query
1828 ALLOCATE (work(1))
1829 info = 0
1830 ALLOCATE (h_mat_copy(nforce_eval, nforce_eval), s_mat_copy(nforce_eval, nforce_eval))
1831 h_mat_copy(:, :) = h_mat(:, :) ! Need explicit copies because dsygv destroys original values
1832 s_mat_copy(:, :) = s_mat(:, :)
1833 CALL dsygv(1, 'V', 'U', nforce_eval, h_mat_copy, nforce_eval, s_mat_copy, nforce_eval, eigenv, work, -1, info)
1834 work_array_size = nint(work(1))
1835 DEALLOCATE (h_mat_copy, s_mat_copy)
1836 ! Allocate work array
1837 DEALLOCATE (work)
1838 ALLOCATE (work(work_array_size))
1839 work = 0.0_dp
1840 ! Solve Hc = eSc
1841 info = 0
1842 CALL dsygv(1, 'V', 'U', nforce_eval, h_mat, nforce_eval, s_mat, nforce_eval, eigenv, work, work_array_size, info)
1843 IF (info /= 0) THEN
1844 IF (info > nforce_eval) THEN
1845 cpabort("Matrix S is not positive definite")
1846 ELSE
1847 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 DO ivar = 1, nforce_eval
1853 h_mat(:, ivar) = h_mat(:, ivar)/dnrm2(nforce_eval, h_mat(:, ivar), 1)
1854 END DO
1855 DEALLOCATE (work)
1856 IF (iounit > 0) THEN
1857 WRITE (iounit, '(/,T3,A)') '------------------ CDFT Configuration Interaction (CDFT-CI) ------------------'
1858 DO ivar = 1, nforce_eval
1859 IF (ivar == 1) THEN
1860 WRITE (iounit, '(T3,A,T58,(3X,F20.14))') 'Ground state energy:', eigenv(ivar)
1861 ELSE
1862 WRITE (iounit, '(/,T3,A,I2,A,T58,(3X,F20.14))') 'Excited state (', ivar - 1, ' ) energy:', eigenv(ivar)
1863 END IF
1864 DO istate = 1, nforce_eval, 2
1865 IF (istate == 1) THEN
1866 WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
1867 'Expansion coefficients:', h_mat(istate, ivar), h_mat(istate + 1, ivar)
1868 ELSE IF (istate < nforce_eval) THEN
1869 WRITE (iounit, '(T54,(3X,2F12.6))') h_mat(istate, ivar), h_mat(istate + 1, ivar)
1870 ELSE
1871 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 '------------------------------------------------------------------------------'
1877 END IF
1878 DEALLOCATE (s_mat, h_mat, eigenv)
1879 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1880 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1881 CALL timestop(handle)
1882
1883 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 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 TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: blocks
1901 TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1902 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 NULLIFY (logger, force_env_section, print_section, mixed_cdft)
1910
1911 cpassert(ASSOCIATED(force_env))
1912 CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
1913 cpassert(ASSOCIATED(mixed_cdft))
1914
1915 IF (.NOT. mixed_cdft%block_diagonalize) RETURN
1916
1917 logger => cp_get_default_logger()
1918 CALL timeset(routinen, handle)
1919
1920 cpassert(ALLOCATED(mixed_cdft%results%S))
1921 cpassert(ALLOCATED(mixed_cdft%results%H))
1922 nforce_eval = SIZE(mixed_cdft%results%S, 1)
1923
1924 CALL force_env_get(force_env=force_env, &
1925 force_env_section=force_env_section)
1926 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1927 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
1928 ! Read block definitions from input
1929 CALL mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
1930 nblk = SIZE(blocks)
1931 ! Start block diagonalization
1932 DO irecursion = 1, nrecursion
1933 ! Print block definitions
1934 IF (iounit > 0 .AND. irecursion == 1) THEN
1935 WRITE (iounit, '(/,T3,A)') '-------------------------- CDFT BLOCK DIAGONALIZATION ------------------------'
1936 WRITE (iounit, '(T3,A)') 'Block diagonalizing the mixed CDFT Hamiltonian'
1937 WRITE (iounit, '(T3,A,I3)') 'Number of blocks:', nblk
1938 WRITE (iounit, '(T3,A,L3)') 'Ignoring excited states within blocks:', ignore_excited
1939 WRITE (iounit, '(/,T3,A)') 'List of CDFT states for each block'
1940 DO i = 1, nblk
1941 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 IF (irecursion > 1) THEN
1946 nblk = nblk/2
1947 ALLOCATE (blocks(nblk))
1948 j = 1
1949 DO i = 1, nblk
1950 NULLIFY (blocks(i)%array)
1951 ALLOCATE (blocks(i)%array(2))
1952 blocks(i)%array = [j, j + 1]
1953 j = j + 2
1954 END DO
1955 ! Print info
1956 IF (iounit > 0) THEN
1957 WRITE (iounit, '(/, T3,A)') 'Recursive block diagonalization of the mixed CDFT Hamiltonian'
1958 WRITE (iounit, '(T6,A)') 'Block diagonalization is continued until only two matrix blocks remain.'
1959 WRITE (iounit, '(T6,A)') 'The new blocks are formed by collecting pairs of blocks from the previous'
1960 WRITE (iounit, '(T6,A)') 'block diagonalized matrix in ascending order.'
1961 WRITE (iounit, '(/,T3,A,I3,A,I3)') 'Recursion step:', irecursion - 1, ' of ', nrecursion - 1
1962 WRITE (iounit, '(/,T3,A)') 'List of old block indices for each new block'
1963 DO i = 1, nblk
1964 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 CALL mixed_cdft_get_blocks(mixed_cdft, blocks, h_block, s_block)
1970 ! Diagonalize blocks
1971 CALL mixed_cdft_diagonalize_blocks(blocks, h_block, s_block, eigenvalues)
1972 ! Assemble the block diagonalized matrices
1973 IF (ignore_excited) THEN
1974 n = nblk
1975 ELSE
1976 n = nforce_eval
1977 END IF
1978 CALL mixed_cdft_assemble_block_diag(mixed_cdft, blocks, h_block, eigenvalues, n, iounit)
1979 ! Deallocate work
1980 DO i = 1, nblk
1981 DEALLOCATE (h_block(i)%array)
1982 DEALLOCATE (s_block(i)%array)
1983 DEALLOCATE (eigenvalues(i)%array)
1984 DEALLOCATE (blocks(i)%array)
1985 END DO
1986 DEALLOCATE (h_block, s_block, eigenvalues, blocks)
1987 END DO ! recursion
1988 IF (iounit > 0) THEN
1989 WRITE (iounit, '(T3,A)') &
1990 '------------------------------------------------------------------------------'
1991 END IF
1992 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
1993 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
1994 CALL timestop(handle)
1995
1996 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
2018 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: metric
2019 TYPE(dbcsr_type) :: e_vectors
2020
2021 CALL timeset(routinen, handle)
2022 nforce_eval = SIZE(mixed_cdft%results%H, 1)
2023 npermutations = nforce_eval*(nforce_eval - 1)/2
2024 nspins = SIZE(density_matrix_diff, 2)
2025 ALLOCATE (metric(npermutations, nspins))
2026 metric = 0.0_dp
2027 CALL dbcsr_create(e_vectors, template=density_matrix_diff(1, 1)%matrix)
2028 DO ispin = 1, nspins
2029 ALLOCATE (evals(ncol_mo(ispin)))
2030 DO ipermutation = 1, npermutations
2031 ! Take into account doubly occupied orbitals without LSD
2032 IF (nspins == 1) THEN
2033 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 para_env=force_env%para_env, blacs_env=mixed_cdft%blacs_env)
2038 CALL dbcsr_release_p(density_matrix_diff(ipermutation, ispin)%matrix)
2039 DO j = 1, ncol_mo(ispin)
2040 metric(ipermutation, ispin) = metric(ipermutation, ispin) + (evals(j)**2 - evals(j)**4)
2041 END DO
2042 END DO
2043 DEALLOCATE (evals)
2044 END DO
2045 CALL dbcsr_release(e_vectors)
2046 DEALLOCATE (density_matrix_diff)
2047 metric(:, :) = metric(:, :)/4.0_dp
2048 CALL mixed_cdft_result_type_set(mixed_cdft%results, metric=metric)
2049 DEALLOCATE (metric)
2050 CALL timestop(handle)
2051
2052 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coupling_wfn
2077 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: overlaps
2078 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 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 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_set
2086 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2087 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2088 TYPE(section_vals_type), POINTER :: force_env_section, mixed_cdft_section
2089
2090 NULLIFY (mixed_cdft_section, subsys_mix, particle_set, qs_kind_set, atomic_kind_set, &
2091 mixed_mo_coeff, mixed_matrix_s, force_env_section)
2092 logger => cp_get_default_logger()
2093
2094 CALL timeset(routinen, handle)
2095 nforce_eval = SIZE(mixed_cdft%results%H, 1)
2096 npermutations = nforce_eval*(nforce_eval - 1)/2
2097 nspins = SIZE(nrow_mo)
2098 mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
2099 mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
2100 CALL force_env_get(force_env=force_env, &
2101 force_env_section=force_env_section)
2102 ! Create mo_set for input wfn
2103 ALLOCATE (mo_set(nspins))
2104 IF (nspins == 2) THEN
2105 maxocc = 1.0_dp
2106 ELSE
2107 maxocc = 2.0_dp
2108 END IF
2109 DO ispin = 1, nspins
2110 nao = nrow_mo(ispin)
2111 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 flexible_electron_count=0.0_dp)
2116 CALL set_mo_set(mo_set(ispin), uniform_occupation=.true., homo=nmo)
2117 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 name="GS_MO_COEFF"//trim(adjustl(cp_to_string(ispin)))//"MATRIX")
2121 ALLOCATE (mo_set(ispin)%eigenvalues(nmo))
2122 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 IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
2126 ! This really shouldnt be a problem?
2127 CALL cp_abort(__location__, &
2128 "QMMM + wavefunction overlap method not supported.")
2129 END IF
2130 CALL force_env_get(force_env=force_env, subsys=subsys_mix)
2131 mixed_cdft_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT")
2132 CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set, particle_set=particle_set)
2133 cpassert(ASSOCIATED(mixed_cdft%qs_kind_set))
2134 IF (force_env%para_env%is_source()) THEN
2135 CALL wfn_restart_file_name(file_name, exist, mixed_cdft_section, logger)
2136 END IF
2137 CALL force_env%para_env%bcast(exist)
2138 CALL force_env%para_env%bcast(file_name)
2139 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 " 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 cdft=.true.)
2150 IF (natom_mismatch) THEN
2151 CALL cp_abort(__location__, &
2152 "Restart wfn file has a wrong number of atoms")
2153 END IF
2154 ! Orthonormalize wfn
2155 DO ispin = 1, nspins
2156 IF (mixed_cdft%has_unit_metric) THEN
2157 CALL make_basis_simple(mo_set(ispin)%mo_coeff, ncol_mo(ispin))
2158 ELSE
2159 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 ALLOCATE (coupling_wfn(npermutations))
2164 ALLOCATE (overlaps(2, npermutations, nspins))
2165 overlaps = 0.0_dp
2166 DO ispin = 1, nspins
2167 ! Allocate work
2168 nao = nrow_mo(ispin)
2169 nmo = ncol_mo(ispin)
2170 CALL cp_fm_struct_create(mo_mo_fmstruct, nrow_global=nmo, ncol_global=nmo, &
2171 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 name="MO_OVERLAP_MATRIX_WFN")
2174 CALL cp_fm_create(matrix=inverse_mat, matrix_struct=mo_mo_fmstruct, &
2175 name="INVERSE_MO_OVERLAP_MATRIX_WFN")
2176 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 name="OVERLAP_MO_COEFF_WFN")
2180 DO ipermutation = 1, npermutations
2181 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 mo_tmp, nmo, 1.0_dp, 0.0_dp)
2185 ! C_i^T * (S*C_r)
2186 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 mo_tmp, 0.0_dp, mo_overlap_wfn)
2190 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 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 mo_tmp, 0.0_dp, mo_overlap_wfn)
2196 CALL cp_fm_invert(mo_overlap_wfn, inverse_mat, overlaps(2, ipermutation, ispin), eps_svd=mixed_cdft%eps_svd)
2197 END DO
2198 CALL cp_fm_release(mo_overlap_wfn)
2199 CALL cp_fm_release(inverse_mat)
2200 CALL cp_fm_release(mo_tmp)
2201 CALL deallocate_mo_set(mo_set(ispin))
2202 END DO
2203 DEALLOCATE (mo_set)
2204 DO ipermutation = 1, npermutations
2205 CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
2206 IF (nspins == 2) THEN
2207 overlaps(1, ipermutation, 1) = abs(overlaps(1, ipermutation, 1)*overlaps(1, ipermutation, 2)) ! A in eq. 12c
2208 overlaps(2, ipermutation, 1) = abs(overlaps(2, ipermutation, 1)*overlaps(2, ipermutation, 2)) ! B in eq. 12c
2209 ELSE
2210 overlaps(1, ipermutation, 1) = overlaps(1, ipermutation, 1)**2
2211 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 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 "density is fully delocalized in the unconstrained ground state.")
2220 coupling_wfn(ipermutation) = 0.0_dp
2221 ELSE
2222 energy_diff = mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate)
2223 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 sda))
2230 END IF
2231 END DO
2232 DEALLOCATE (overlaps)
2233 CALL mixed_cdft_result_type_set(mixed_cdft%results, wfn=coupling_wfn)
2234 DEALLOCATE (coupling_wfn)
2235 CALL timestop(handle)
2236
2237 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: catom
2259 LOGICAL :: in_memory, store_vectors
2260 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint
2261 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients
2262 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: position_vecs, r12
2263 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 NULLIFY (mixed_env, mixed_cdft)
2269 store_vectors = .true.
2270 logger => cp_get_default_logger()
2271 CALL timeset(routinen, handle)
2272 mixed_env => force_env%mixed_env
2273 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 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 coefficients, catom)
2282 CALL timestop(handle)
2283
2284 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 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: radii_list
2334 REAL(kind=dp), DIMENSION(:, :), POINTER :: pab
2335 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2344 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2345 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 NULLIFY (pab, cell, force_env_qs, particle_set, force_env_section, print_section, &
2350 qs_kind_set, particles, subsys_mix, rs_cavity, cavity_env, auxbas_pw_pool, &
2351 atomic_kind_set, radii_list, cdft_control)
2352 logger => cp_get_default_logger()
2353 nforce_eval = SIZE(force_env%sub_force_env)
2354 CALL timeset(routinen, handle)
2355 CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2356 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 cell=cell)
2360 CALL cp_subsys_get(subsys=subsys_mix, &
2361 particles=particles, &
2362 particle_set=particle_set)
2363 ELSE
2364 DO iforce_eval = 1, nforce_eval
2365 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
2366 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 cell=cell)
2371 CALL cp_subsys_get(subsys=subsys_mix, &
2372 particles=particles, &
2373 particle_set=particle_set)
2374 END IF
2375 natom = SIZE(particles%els)
2376 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2377 cdft_control => mixed_cdft%cdft_control
2378 cpassert(ASSOCIATED(cdft_control))
2379 IF (.NOT. ASSOCIATED(cdft_control%becke_control%cutoffs)) THEN
2380 CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2381 ALLOCATE (cdft_control%becke_control%cutoffs(natom))
2382 SELECT CASE (cdft_control%becke_control%cutoff_type)
2383 CASE (becke_cutoff_global)
2384 cdft_control%becke_control%cutoffs(:) = cdft_control%becke_control%rglobal
2385 CASE (becke_cutoff_element)
2386 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 "not match number of atomic kinds in the input coordinate file.")
2390 END IF
2391 DO ikind = 1, SIZE(atomic_kind_set)
2392 CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2393 DO iatom = 1, katom
2394 atom_a = atom_list(iatom)
2395 cdft_control%becke_control%cutoffs(atom_a) = cdft_control%becke_control%cutoffs_tmp(ikind)
2396 END DO
2397 END DO
2398 DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
2399 END SELECT
2400 END IF
2401 build = .false.
2402 IF (cdft_control%becke_control%adjust .AND. &
2403 .NOT. ASSOCIATED(cdft_control%becke_control%aij)) THEN
2404 ALLOCATE (cdft_control%becke_control%aij(natom, natom))
2405 build = .true.
2406 END IF
2407 ALLOCATE (catom(cdft_control%natoms))
2408 IF (cdft_control%save_pot .OR. &
2409 cdft_control%becke_control%cavity_confine .OR. &
2410 cdft_control%becke_control%should_skip .OR. &
2411 mixed_cdft%first_iteration) THEN
2412 ALLOCATE (is_constraint(natom))
2413 is_constraint = .false.
2414 END IF
2415 in_memory = calculate_forces .AND. cdft_control%becke_control%in_memory
2416 IF (in_memory .NEQV. calculate_forces) THEN
2417 CALL cp_abort(__location__, &
2418 "The flag BECKE_CONSTRAINT\IN_MEMORY must be activated "// &
2419 "for the calculation of mixed CDFT forces")
2420 END IF
2421 IF (in_memory .OR. mixed_cdft%first_iteration) ALLOCATE (coefficients(natom))
2422 DO i = 1, cdft_control%natoms
2423 catom(i) = cdft_control%atoms(i)
2424 IF (cdft_control%save_pot .OR. &
2425 cdft_control%becke_control%cavity_confine .OR. &
2426 cdft_control%becke_control%should_skip .OR. &
2427 mixed_cdft%first_iteration) THEN
2428 is_constraint(catom(i)) = .true.
2429 END IF
2430 IF (in_memory .OR. mixed_cdft%first_iteration) THEN
2431 coefficients(catom(i)) = cdft_control%group(1)%coeff(i)
2432 END IF
2433 END DO
2434 CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
2435 bo = auxbas_pw_pool%pw_grid%bounds_local
2436 np = auxbas_pw_pool%pw_grid%npts
2437 dr = auxbas_pw_pool%pw_grid%dr
2438 shift = -real(modulo(np, 2), dp)*dr/2.0_dp
2439 IF (store_vectors) THEN
2440 IF (in_memory) ALLOCATE (pair_dist_vecs(3, natom, natom))
2441 ALLOCATE (position_vecs(3, natom))
2442 END IF
2443 DO i = 1, 3
2444 cell_v(i) = cell%hmat(i, i)
2445 END DO
2446 ALLOCATE (r12(natom, natom))
2447 DO iatom = 1, natom - 1
2448 DO jatom = iatom + 1, natom
2449 r = particle_set(iatom)%r
2450 r1 = particle_set(jatom)%r
2451 DO i = 1, 3
2452 r(i) = modulo(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2453 r1(i) = modulo(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
2454 END DO
2455 dist_vec = (r - r1) - anint((r - r1)/cell_v)*cell_v
2456 IF (store_vectors) THEN
2457 position_vecs(:, iatom) = r(:)
2458 IF (iatom == 1 .AND. jatom == natom) position_vecs(:, jatom) = r1(:)
2459 IF (in_memory) THEN
2460 pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
2461 pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
2462 END IF
2463 END IF
2464 r12(iatom, jatom) = norm2(dist_vec)
2465 r12(jatom, iatom) = r12(iatom, jatom)
2466 IF (build) THEN
2467 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2468 kind_number=ikind)
2469 ircov = cdft_control%becke_control%radii(ikind)
2470 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
2471 kind_number=ikind)
2472 jrcov = cdft_control%becke_control%radii(ikind)
2473 IF (ircov /= jrcov) THEN
2474 chi = ircov/jrcov
2475 uij = (chi - 1.0_dp)/(chi + 1.0_dp)
2476 cdft_control%becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp)
2477 IF (cdft_control%becke_control%aij(iatom, jatom) &
2478 > 0.5_dp) THEN
2479 cdft_control%becke_control%aij(iatom, jatom) = 0.5_dp
2480 ELSE IF (cdft_control%becke_control%aij(iatom, jatom) &
2481 < -0.5_dp) THEN
2482 cdft_control%becke_control%aij(iatom, jatom) = -0.5_dp
2483 END IF
2484 ELSE
2485 cdft_control%becke_control%aij(iatom, jatom) = 0.0_dp
2486 END IF
2487 cdft_control%becke_control%aij(jatom, iatom) = &
2488 -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 IF (mixed_cdft%first_iteration) THEN
2494 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2495 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2496 IF (iounit > 0) THEN
2497 WRITE (iounit, '(/,T3,A,T66)') &
2498 '-------------------------- Becke atomic parameters ---------------------------'
2499 IF (cdft_control%becke_control%adjust) THEN
2500 WRITE (iounit, '(T3,A,A)') &
2501 'Atom Element Coefficient', ' Cutoff (angstrom) CDFT Radius (angstrom)'
2502 DO iatom = 1, natom
2503 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2504 element_symbol=element_symbol, &
2505 kind_number=ikind)
2506 ircov = cp_unit_from_cp2k(cdft_control%becke_control%radii(ikind), "angstrom")
2507 IF (is_constraint(iatom)) THEN
2508 coef = coefficients(iatom)
2509 ELSE
2510 coef = 0.0_dp
2511 END IF
2512 WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3,T73,F8.3)") &
2513 iatom, adjustr(element_symbol), coef, &
2514 cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom"), &
2515 ircov
2516 END DO
2517 ELSE
2518 WRITE (iounit, '(T3,A,A)') &
2519 'Atom Element Coefficient', ' Cutoff (angstrom)'
2520 DO iatom = 1, natom
2521 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
2522 element_symbol=element_symbol)
2523 IF (is_constraint(iatom)) THEN
2524 coef = coefficients(iatom)
2525 ELSE
2526 coef = 0.0_dp
2527 END IF
2528 WRITE (iounit, "(i6,T14,A2,T22,F8.3,T44,F8.3)") &
2529 iatom, adjustr(element_symbol), coef, &
2530 cp_unit_from_cp2k(cdft_control%becke_control%cutoffs(iatom), "angstrom")
2531 END DO
2532 END IF
2533 WRITE (iounit, '(T3,A)') &
2534 '------------------------------------------------------------------------------'
2535 END IF
2536 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
2537 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2538 mixed_cdft%first_iteration = .false.
2539 END IF
2540
2541 IF (cdft_control%becke_control%cavity_confine) THEN
2542 cpassert(ASSOCIATED(mixed_cdft%qs_kind_set))
2543 cavity_env => cdft_control%becke_control%cavity_env
2544 qs_kind_set => mixed_cdft%qs_kind_set
2545 CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
2546 nkind = SIZE(qs_kind_set)
2547 IF (.NOT. ASSOCIATED(cavity_env%kind_shape_fn)) THEN
2548 IF (ASSOCIATED(cdft_control%becke_control%radii)) THEN
2549 ALLOCATE (radii_list(SIZE(cdft_control%becke_control%radii)))
2550 DO ikind = 1, SIZE(cdft_control%becke_control%radii)
2551 IF (cavity_env%use_bohr) THEN
2552 radii_list(ikind) = cdft_control%becke_control%radii(ikind)
2553 ELSE
2554 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 radii_list=radii_list)
2561 IF (ASSOCIATED(radii_list)) THEN
2562 DEALLOCATE (radii_list)
2563 END IF
2564 END IF
2565 NULLIFY (rs_cavity)
2566 CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_rs_grid=rs_cavity, &
2567 auxbas_pw_pool=auxbas_pw_pool)
2568 ! be careful in parallel nsmax is chosen with multigrid in mind!
2569 CALL rs_grid_zero(rs_cavity)
2570 ALLOCATE (pab(1, 1))
2571 nthread = 1
2572 ithread = 0
2573 DO ikind = 1, SIZE(atomic_kind_set)
2574 numexp = cavity_env%kind_shape_fn(ikind)%numexp
2575 IF (numexp <= 0) cycle
2576 CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
2577 ALLOCATE (cores(katom))
2578 DO iex = 1, numexp
2579 alpha = cavity_env%kind_shape_fn(ikind)%zet(iex)
2580 coef = cavity_env%kind_shape_fn(ikind)%coef(iex)
2581 npme = 0
2582 cores = 0
2583 DO iatom = 1, katom
2584 IF (rs_cavity%desc%parallel .AND. .NOT. rs_cavity%desc%distributed) THEN
2585 ! replicated realspace grid, split the atoms up between procs
2586 IF (modulo(iatom, rs_cavity%desc%group_size) == rs_cavity%desc%my_pos) THEN
2587 npme = npme + 1
2588 cores(npme) = iatom
2589 END IF
2590 ELSE
2591 npme = npme + 1
2592 cores(npme) = iatom
2593 END IF
2594 END DO
2595 DO j = 1, npme
2596 iatom = cores(j)
2597 atom_a = atom_list(iatom)
2598 pab(1, 1) = coef
2599 IF (store_vectors) THEN
2600 ra(:) = position_vecs(:, atom_a) + cell_v(:)/2._dp
2601 ELSE
2602 ra(:) = pbc(particle_set(atom_a)%r, cell)
2603 END IF
2604 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 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 subpatch_pattern=0)
2617 END IF
2618 END DO
2619 END DO
2620 DEALLOCATE (cores)
2621 END DO
2622 DEALLOCATE (pab)
2623 CALL auxbas_pw_pool%create_pw(cdft_control%becke_control%cavity)
2624 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 just_zero=.false., bounds=bounds, work=my_work)
2628 IF (bounds(2) < bo(2, 3)) THEN
2629 bounds(2) = bounds(2) - 1
2630 ELSE
2631 bounds(2) = bo(2, 3)
2632 END IF
2633 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 bounds(1) = bounds(1) + 1
2638 ELSE
2639 bounds(1) = bo(1, 3)
2640 END IF
2641 IF (bounds(1) > bounds(2)) THEN
2642 my_work_size = 0
2643 ELSE
2644 my_work_size = (bounds(2) - bounds(1) + 1)
2645 IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2646 my_work_size = my_work_size*(bo(2, 2) - bo(1, 2) + 1)
2647 ELSE
2648 my_work_size = my_work_size*(bo(2, 1) - bo(1, 1) + 1)
2649 END IF
2650 END IF
2651 cdft_control%becke_control%confine_bounds = bounds
2652 IF (cdft_control%becke_control%print_cavity) THEN
2653 CALL hfun_zero(cdft_control%becke_control%cavity%array, &
2654 cdft_control%becke_control%eps_cavity, just_zero=.true.)
2655 NULLIFY (stride)
2656 ALLOCATE (stride(3))
2657 stride = [2, 2, 2]
2658 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 log_filename=.false., mpi_io=mpi_io)
2663 IF (force_env%para_env%is_source() .AND. unit_nr < 1) THEN
2664 CALL cp_abort(__location__, &
2665 "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 stride=stride, mpi_io=mpi_io)
2670 CALL cp_print_key_finished_output(unit_nr, logger, print_section, '', mpi_io=mpi_io)
2671 DEALLOCATE (stride)
2672 END IF
2673 END IF
2674 bo_conf = bo
2675 IF (cdft_control%becke_control%cavity_confine) THEN
2676 bo_conf(:, 3) = cdft_control%becke_control%confine_bounds
2677 END IF
2678 ! Load balance
2679 IF (mixed_cdft%dlb) THEN
2680 CALL mixed_becke_constraint_dlb(force_env, mixed_cdft, my_work, &
2681 my_work_size, natom, bo, bo_conf)
2682 END IF
2683 ! The bounds have been finalized => time to allocate storage for working matrices
2684 offset_dlb = 0
2685 IF (mixed_cdft%dlb) THEN
2686 IF (mixed_cdft%dlb_control%send_work .AND. .NOT. mixed_cdft%is_special) THEN
2687 offset_dlb = sum(mixed_cdft%dlb_control%target_list(2, :))
2688 END IF
2689 END IF
2690 IF (cdft_control%becke_control%cavity_confine) THEN
2691 ! Get rid of the zero part of the confinement cavity (cr3d -> real(:,:,:))
2692 IF (mixed_cdft%is_special) THEN
2693 ALLOCATE (mixed_cdft%sendbuff(SIZE(mixed_cdft%dest_list)))
2694 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 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 bo_conf(1, 3):bo_conf(2, 3))
2701 END DO
2702 ELSE IF (mixed_cdft%is_pencil) THEN
2703 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 bo_conf(1, 3):bo_conf(2, 3))
2707 ELSE
2708 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 bo_conf(1, 3):bo_conf(2, 3))
2712 END IF
2713 CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
2714 END IF
2715 IF (mixed_cdft%is_special) THEN
2716 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 bo(1, 2):bo(2, 2), bo_conf(1, 3):bo_conf(2, 3)))
2719 mixed_cdft%sendbuff(i)%weight = 0.0_dp
2720 END DO
2721 ELSE IF (mixed_cdft%is_pencil) THEN
2722 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 mixed_cdft%weight = 0.0_dp
2724 ELSE
2725 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 mixed_cdft%weight = 0.0_dp
2727 END IF
2728 IF (in_memory) THEN
2729 IF (mixed_cdft%is_special) THEN
2730 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 bo_conf(1, 3):bo_conf(2, 3)))
2735 mixed_cdft%sendbuff(i)%gradients = 0.0_dp
2736 END DO
2737 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 bo_conf(1, 3):bo_conf(2, 3)))
2741 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 bo_conf(1, 3):bo_conf(2, 3)))
2746 cdft_control%group(1)%gradients = 0.0_dp
2747 END IF
2748 END IF
2749
2750 CALL timestop(handle)
2751
2752 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 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 INTEGER, DIMENSION(:), POINTER :: buffsize, cumulative_work, expected_work, load_imbalance, &
2783 nrecv, nsend_proc, sendbuffer, should_warn, tmp, work_index, work_size
2784 INTEGER, DIMENSION(:, :), POINTER :: targets, tmp_bo
2785 LOGICAL :: consistent
2786 LOGICAL, DIMENSION(:), POINTER :: mask_recv, mask_send, touched
2787 REAL(kind=dp) :: average_work, load_scale, &
2788 very_overloaded, work_factor
2789 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cavity
2790 TYPE(buffers_bi), DIMENSION(:), POINTER :: recvbuffer, sbuff
2791 TYPE(cdft_control_type), POINTER :: cdft_control
2792 TYPE(cp_logger_type), POINTER :: logger
2793 TYPE(mp_request_type), DIMENSION(4) :: req
2794 TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
2795 TYPE(section_vals_type), POINTER :: force_env_section, print_section
2796
2797 logger => cp_get_default_logger()
2798 CALL timeset(routinen, handle)
2799 mixed_cdft%dlb_control%recv_work = .false.
2800 mixed_cdft%dlb_control%send_work = .false.
2801 NULLIFY (expected_work, work_index, load_imbalance, work_size, &
2802 cumulative_work, sendbuffer, buffsize, req_recv, req_total, &
2803 tmp, nrecv, nsend_proc, targets, tmp_bo, touched, &
2804 mask_recv, mask_send, cavity, recvbuffer, sbuff, force_env_section, &
2805 print_section, cdft_control)
2806 CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
2807 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
2808 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
2809 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 load_scale = mixed_cdft%dlb_control%load_scale
2815 very_overloaded = mixed_cdft%dlb_control%very_overloaded
2816 more_work = mixed_cdft%dlb_control%more_work
2817 max_targets = 40
2818 work_factor = 0.8_dp
2819 ! Reset targets/sources
2820 IF (mixed_cdft%is_special) THEN
2821 DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo, &
2822 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 mixed_cdft%source_list_bo(SIZE(mixed_cdft%source_bo_save, 1), SIZE(mixed_cdft%source_bo_save, 2)))
2827 mixed_cdft%dest_list = mixed_cdft%dest_list_save
2828 mixed_cdft%source_list = mixed_cdft%source_list_save
2829 mixed_cdft%dest_list_bo = mixed_cdft%dest_bo_save
2830 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 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 expected_work = 0
2840 expected_work(force_env%para_env%mepos + 1) = my_work
2841 work_size = 0
2842 work_size(force_env%para_env%mepos + 1) = my_work_size
2843 IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2844 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 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 REAL(bo(2, 2) - bo(1, 2) + 1, dp))
2854 END IF
2855 END IF
2856 CALL force_env%para_env%sum(expected_work)
2857 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 mixed_cdft%dlb_control%expected_work = expected_work
2860 ! Take into account the prediction error of the last step
2861 IF (ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
2862 expected_work = expected_work - mixed_cdft%dlb_control%prediction_error
2863 END IF
2864 !
2865 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 targets(2, force_env%para_env%num_pe))
2869 load_imbalance = expected_work - nint(average_work)
2870 no_overloaded = 0
2871 no_underloaded = 0
2872 targets = 0
2873 ! Convert the load imbalance to a multiple of the actual work size
2874 DO i = 1, force_env%para_env%num_pe
2875 IF (load_imbalance(i) > 0) THEN
2876 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 IF (expected_work(i) > nint(very_overloaded*average_work)) THEN
2879 load_imbalance(i) = (ceiling(real(load_imbalance(i), dp)/real(work_size(i), dp)) + more_work)*work_size(i)
2880 ELSE
2881 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 load_imbalance(i) = nint(load_imbalance(i)*load_scale)
2887 no_underloaded = no_underloaded + 1
2888 END IF
2889 END DO
2890 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 IF (load_imbalance(force_env%para_env%mepos + 1) > 0) THEN
2894 offset = 0
2895 mixed_cdft%dlb_control%send_work = .true.
2896 ! Build up the total amount of work that needs redistribution
2897 ALLOCATE (cumulative_work(force_env%para_env%num_pe))
2898 cumulative_work = 0
2899 DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
2900 IF (work_index(i) == force_env%para_env%mepos + 1) THEN
2901 EXIT
2902 ELSE
2903 offset = offset + load_imbalance(work_index(i))
2904 IF (i == force_env%para_env%num_pe) THEN
2905 cumulative_work(i) = load_imbalance(work_index(i))
2906 ELSE
2907 cumulative_work(i) = cumulative_work(i + 1) + load_imbalance(work_index(i))
2908 END IF
2909 END IF
2910 END DO
2911 my_pos = i
2912 j = force_env%para_env%num_pe
2913 nsend_max = load_imbalance(work_index(j))/work_size(work_index(j))
2914 exhausted_work = 0
2915 ! Determine send offset by going through all processors that are more overloaded than my_pos
2916 DO i = 1, no_underloaded
2917 IF (my_pos == force_env%para_env%num_pe) EXIT
2918 nsend = -load_imbalance(work_index(i))/work_size(work_index(j))
2919 IF (nsend < 1) nsend = 1
2920 nsend_max = nsend_max - nsend
2921 IF (nsend_max < 0) nsend = nsend + nsend_max
2922 exhausted_work = exhausted_work + nsend*work_size(work_index(j))
2923 offset = offset - nsend*work_size(work_index(j))
2924 IF (offset < 0) EXIT
2925 IF (exhausted_work == cumulative_work(j)) THEN
2926 j = j - 1
2927 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 IF (i > no_underloaded) THEN
2933 i = no_underloaded
2934 END IF
2935 my_target = i
2936 DEALLOCATE (cumulative_work)
2937 ! Determine how much and who to send slices of my grid points
2938 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 IF (mixed_cdft%is_pencil .OR. mixed_cdft%is_special) THEN
2941 nsend_limit = bo(2, 1) - bo(1, 1) + 1
2942 ELSE
2943 nsend_limit = bo(2, 2) - bo(1, 2) + 1
2944 END IF
2945 IF (.NOT. mixed_cdft%is_special) THEN
2946 ALLOCATE (mixed_cdft%dlb_control%target_list(3, max_targets))
2947 ELSE
2948 ALLOCATE (mixed_cdft%dlb_control%target_list(3 + 2*SIZE(mixed_cdft%dest_list), max_targets))
2949 ALLOCATE (touched(SIZE(mixed_cdft%dest_list)))
2950 touched = .false.
2951 END IF
2952 mixed_cdft%dlb_control%target_list = uninitialized
2953 i = 1
2954 ispecial = 1
2955 offset_special = 0
2956 targets(1, my_pos) = my_target
2957 send_total = 0
2958 ! Main loop. Note, we actually allow my_pos to offload more slices than nsend_max
2959 DO
2960 nsend = -load_imbalance(work_index(my_target))/work_size(force_env%para_env%mepos + 1)
2961 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 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 mixed_cdft%dlb_control%target_list(1, i) = work_index(my_target) - 1 ! This is the actual processor rank
2970 IF (mixed_cdft%is_special) THEN
2971 mixed_cdft%dlb_control%target_list(2, i) = 0
2972 actually_sent = nsend
2973 DO j = ispecial, SIZE(mixed_cdft%dest_list)
2974 mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + 1
2975 touched(j) = .true.
2976 IF (nsend < mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1) THEN
2977 mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2978 mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(1, j) + nsend - 1
2979 mixed_cdft%dest_list_bo(1, j) = mixed_cdft%dest_list_bo(1, j) + nsend
2980 nsend = 0
2981 EXIT
2982 ELSE
2983 mixed_cdft%dlb_control%target_list(3 + 2*j - 1, i) = mixed_cdft%dest_list_bo(1, j)
2984 mixed_cdft%dlb_control%target_list(3 + 2*j, i) = mixed_cdft%dest_list_bo(2, j)
2985 nsend = nsend - (mixed_cdft%dest_list_bo(2, j) - mixed_cdft%dest_list_bo(1, j) + 1)
2986 mixed_cdft%dest_list_bo(1:2, j) = should_deallocate
2987 END IF
2988 IF (nsend <= 0) EXIT
2989 END DO
2990 IF (mixed_cdft%dest_list_bo(1, ispecial) == should_deallocate) ispecial = j + 1
2991 actually_sent = actually_sent - nsend
2992 nsend_max = nsend_max - actually_sent
2993 send_total = send_total + actually_sent
2994 ELSE
2995 mixed_cdft%dlb_control%target_list(2, i) = nsend
2996 nsend_max = nsend_max - nsend
2997 send_total = send_total + nsend
2998 END IF
2999 IF (nsend_max < 0) nsend_max = 0
3000 IF (nsend_max == 0) EXIT
3001 IF (my_target /= no_underloaded) THEN
3002 my_target = my_target + 1
3003 ELSE
3004 ! If multiple processors execute this block load balancing will fail
3005 mixed_cdft%dlb_control%target_list(2, i) = mixed_cdft%dlb_control%target_list(2, i) + nsend_max
3006 nsend_max = 0
3007 EXIT
3008 END IF
3009 i = i + 1
3010 IF (i > max_targets) THEN
3011 CALL cp_abort(__location__, &
3012 "Load balancing error: increase max_targets")
3013 END IF
3014 END DO
3015 IF (.NOT. mixed_cdft%is_special) THEN
3016 CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3, 1, i)
3017 ELSE
3018 CALL reallocate(mixed_cdft%dlb_control%target_list, 1, 3 + 2*SIZE(mixed_cdft%dest_list), 1, i)
3019 END IF
3020 targets(2, my_pos) = my_target
3021 ! Equalize the load on the target processors
3022 IF (.NOT. mixed_cdft%is_special) THEN
3023 IF (send_total > nint(work_factor*nsend_limit)) send_total = nint(work_factor*nsend_limit) - 1
3024 nsend = nint(real(send_total, dp)/real(SIZE(mixed_cdft%dlb_control%target_list, 2), dp))
3025 mixed_cdft%dlb_control%target_list(2, :) = nsend
3026 END IF
3027 ELSE
3028 DO i = 1, no_underloaded
3029 IF (work_index(i) == force_env%para_env%mepos + 1) EXIT
3030 END DO
3031 my_pos = i
3032 END IF
3033 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 IF (force_env%para_env%is_source()) THEN
3046 consistent = .true.
3047 DO i = force_env%para_env%num_pe - 1, force_env%para_env%num_pe - no_overloaded + 1, -1
3048 IF (targets(1, i) > no_underloaded) consistent = .false.
3049 IF (targets(1, i) > targets(2, i + 1)) THEN
3050 cycle
3051 ELSE
3052 consistent = .false.
3053 END IF
3054 END DO
3055 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 " 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 IF (my_pos <= no_underloaded) THEN
3072 DO i = force_env%para_env%num_pe, force_env%para_env%num_pe - no_overloaded + 1, -1
3073 IF (targets(1, i) <= my_pos .AND. targets(2, i) >= my_pos) THEN
3074 mixed_cdft%dlb_control%recv_work = .true.
3075 mixed_cdft%dlb_control%my_source = work_index(i) - 1
3076 EXIT
3077 END IF
3078 END DO
3079 IF (mixed_cdft%dlb_control%recv_work) THEN
3080 IF (.NOT. mixed_cdft%is_special) THEN
3081 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 request=req(1))
3084 CALL req(1)%wait()
3085 mixed_cdft%dlb_control%my_dest_repl = [mixed_cdft%dlb_control%bo(11), mixed_cdft%dlb_control%bo(12)]
3086 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 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 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 mixed_cdft%dlb_control%bo(7):mixed_cdft%dlb_control%bo(8)))
3097 mixed_cdft%dlb_control%gradients = 0.0_dp
3098 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 request=req(1))
3101 CALL req(1)%wait()
3102 DEALLOCATE (mixed_cdft%dlb_control%bo)
3103 ELSE
3104 ALLOCATE (buffsize(1))
3105 CALL force_env%para_env%irecv(msgout=buffsize, source=mixed_cdft%dlb_control%my_source, &
3106 request=req(1))
3107 CALL req(1)%wait()
3108 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 request=req(1))
3111 ALLOCATE (mixed_cdft%dlb_control%sendbuff(buffsize(1)))
3112 ALLOCATE (req_recv(buffsize(1)))
3113 DEALLOCATE (buffsize)
3114 CALL req(1)%wait()
3115 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 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 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 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 mixed_cdft%dlb_control%bo(12*(j - 1) + 8)))
3138 mixed_cdft%dlb_control%sendbuff(j)%weight = 0.0_dp
3139 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 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 mixed_cdft%dlb_control%bo(12*(j - 1) + 12)]
3144 END DO
3145 CALL mp_waitall(req_recv)
3146 DEALLOCATE (req_recv)
3147 END IF
3148 END IF
3149 ELSE
3150 IF (.NOT. mixed_cdft%is_special) THEN
3151 offset = 0
3152 ALLOCATE (sendbuffer(12))
3153 send_total = 0
3154 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 (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets] ! Unique communicator tags
3157 mixed_cdft%dlb_control%target_list(3, i) = tags(1)
3158 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 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 mixed_cdft%dest_list(1), mixed_cdft%dest_list(2)]
3169 END IF
3170 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 request=req(1))
3173 CALL req(1)%wait()
3174 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 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 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 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 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 request=req(1))
3197 CALL req(1)%wait()
3198 offset = offset + mixed_cdft%dlb_control%target_list(2, i)
3199 DEALLOCATE (cavity)
3200 END DO
3201 IF (mixed_cdft%is_pencil) THEN
3202 mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 1)
3203 mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 1) + offset - 1
3204 ELSE
3205 mixed_cdft%dlb_control%distributed(1) = bo_conf(1, 2)
3206 mixed_cdft%dlb_control%distributed(2) = bo_conf(1, 2) + offset - 1
3207 END IF
3208 DEALLOCATE (sendbuffer)
3209 ELSE
3210 ALLOCATE (buffsize(1))
3211 DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3212 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 (i - 1)*3 + 1 + 3*max_targets + force_env%para_env%mepos*6*max_targets]
3216 DO j = 4, SIZE(mixed_cdft%dlb_control%target_list, 1)
3217 IF (mixed_cdft%dlb_control%target_list(j, i) > uninitialized) EXIT
3218 END DO
3219 offset_special = j
3220 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 request=req(1))
3224 CALL req(1)%wait()
3225 ALLOCATE (sendbuffer(12*buffsize(1)))
3226 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 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 request=req(1))
3237 CALL req(1)%wait()
3238 DEALLOCATE (sendbuffer)
3239 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 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 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 request=req(1))
3249 CALL req(1)%wait()
3250 DEALLOCATE (cavity)
3251 END DO
3252 END DO
3253 DEALLOCATE (buffsize)
3254 END IF
3255 END IF
3256 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 IF (mixed_cdft%is_special) THEN
3261 my_special_work = 2
3262 ALLOCATE (mask_send(SIZE(mixed_cdft%dest_list)), mask_recv(SIZE(mixed_cdft%source_list)))
3263 ALLOCATE (nsend_proc(SIZE(mixed_cdft%dest_list)), nrecv(SIZE(mixed_cdft%source_list)))
3264 nrecv = 0
3265 nsend_proc = 0
3266 mask_recv = .false.
3267 mask_send = .false.
3268 ELSE
3269 my_special_work = 1
3270 END IF
3271 ALLOCATE (recvbuffer(SIZE(mixed_cdft%source_list)), sbuff(SIZE(mixed_cdft%dest_list)))
3272 ALLOCATE (req_total(my_special_work*SIZE(mixed_cdft%source_list) + (my_special_work**2)*SIZE(mixed_cdft%dest_list)))
3273 ALLOCATE (mixed_cdft%dlb_control%recv_work_repl(SIZE(mixed_cdft%source_list)))
3274 DO i = 1, SIZE(mixed_cdft%source_list)
3275 NULLIFY (recvbuffer(i)%bv, recvbuffer(i)%iv)
3276 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 request=req_total(i), tag=1)
3280 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 tag=2)
3285 END IF
3286 END DO
3287 DO i = 1, my_special_work
3288 DO j = 1, SIZE(mixed_cdft%dest_list)
3289 IF (i == 1) THEN
3290 NULLIFY (sbuff(j)%iv, sbuff(j)%bv)
3291 ALLOCATE (sbuff(j)%bv(1))
3292 sbuff(j)%bv = mixed_cdft%dlb_control%send_work
3293 IF (mixed_cdft%is_special) THEN
3294 ALLOCATE (sbuff(j)%iv(3))
3295 sbuff(j)%iv(1:2) = mixed_cdft%dest_list_bo(1:2, j)
3296 sbuff(j)%iv(3) = 0
3297 IF (sbuff(j)%iv(1) == should_deallocate) mask_send(j) = .true.
3298 IF (mixed_cdft%dlb_control%send_work) THEN
3299 sbuff(j)%bv = touched(j)
3300 IF (touched(j)) THEN
3301 nsend = 0
3302 DO ispecial = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3303 IF (mixed_cdft%dlb_control%target_list(4 + 2*(j - 1), ispecial) /= uninitialized) THEN
3304 nsend = nsend + 1
3305 END IF
3306 END DO
3307 sbuff(j)%iv(3) = nsend
3308 nsend_proc(j) = nsend
3309 END IF
3310 END IF
3311 END IF
3312 END IF
3313 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 request=req_total(ind), tag=1)
3317 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 request=req_total(ind + 2*SIZE(mixed_cdft%dest_list)), tag=2)
3321 END IF
3322 END DO
3323 END DO
3324 CALL mp_waitall(req_total)
3325 DEALLOCATE (req_total)
3326 DO i = 1, SIZE(mixed_cdft%source_list)
3327 mixed_cdft%dlb_control%recv_work_repl(i) = recvbuffer(i)%bv(1)
3328 IF (mixed_cdft%is_special .AND. mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3329 mixed_cdft%source_list_bo(1:2, i) = recvbuffer(i)%iv(1:2)
3330 nrecv(i) = recvbuffer(i)%iv(3)
3331 IF (recvbuffer(i)%iv(1) == should_deallocate) mask_recv(i) = .true.
3332 END IF
3333 DEALLOCATE (recvbuffer(i)%bv)
3334 IF (ASSOCIATED(recvbuffer(i)%iv)) DEALLOCATE (recvbuffer(i)%iv)
3335 END DO
3336 DO j = 1, SIZE(mixed_cdft%dest_list)
3337 DEALLOCATE (sbuff(j)%bv)
3338 IF (ASSOCIATED(sbuff(j)%iv)) DEALLOCATE (sbuff(j)%iv)
3339 END DO
3340 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 IF (.NOT. mixed_cdft%is_special) THEN
3348 IF (mixed_cdft%dlb_control%send_work) THEN
3349 ALLOCATE (req_total(count(mixed_cdft%dlb_control%recv_work_repl) + 2))
3350 ALLOCATE (sendbuffer(6))
3351 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 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 bo_conf(1, 2), bo_conf(1, 1), bo_conf(2, 1)]
3357 END IF
3358 ELSE IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3359 ALLOCATE (req_total(count(mixed_cdft%dlb_control%recv_work_repl)))
3360 END IF
3361 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3362 ALLOCATE (mixed_cdft%dlb_control%recv_info(2))
3363 NULLIFY (mixed_cdft%dlb_control%recv_info(1)%target_list, mixed_cdft%dlb_control%recv_info(2)%target_list)
3364 ALLOCATE (mixed_cdft%dlb_control%recvbuff(2))
3365 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 IF (mixed_cdft%dlb_control%send_work) THEN
3369 ind = count(mixed_cdft%dlb_control%recv_work_repl) + 1
3370 DO i = 1, 2
3371 CALL force_env%para_env%isend(msgin=sendbuffer, &
3372 dest=mixed_cdft%dest_list(i), &
3373 request=req_total(ind))
3374 ind = ind + 1
3375 END DO
3376 END IF
3377 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3378 ind = 1
3379 DO i = 1, 2
3380 IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3381 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 request=req_total(ind))
3385 ind = ind + 1
3386 END IF
3387 END DO
3388 END IF
3389 IF (ASSOCIATED(req_total)) THEN
3390 CALL mp_waitall(req_total)
3391 END IF
3392 ! Now communicate which processor handles which grid points
3393 IF (mixed_cdft%dlb_control%send_work) THEN
3394 ind = count(mixed_cdft%dlb_control%recv_work_repl) + 1
3395 DO i = 1, 2
3396 IF (i == 2) THEN
3397 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 request=req_total(ind))
3402 ind = ind + 1
3403 END DO
3404 END IF
3405 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3406 ind = 1
3407 DO i = 1, 2
3408 IF (mixed_cdft%dlb_control%recv_work_repl(i)) THEN
3409 ALLOCATE (mixed_cdft%dlb_control%recv_info(i)% &
3410 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 request=req_total(ind))
3414 ind = ind + 1
3415 END IF
3416 END DO
3417 END IF
3418 IF (ASSOCIATED(req_total)) THEN
3419 CALL mp_waitall(req_total)
3420 DEALLOCATE (req_total)
3421 END IF
3422 IF (ASSOCIATED(sendbuffer)) DEALLOCATE (sendbuffer)
3423 ELSE
3424 IF (mixed_cdft%dlb_control%send_work) THEN
3425 ALLOCATE (req_total(count(mixed_cdft%dlb_control%recv_work_repl) + 2*count(touched)))
3426 ELSE IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3427 ALLOCATE (req_total(count(mixed_cdft%dlb_control%recv_work_repl)))
3428 END IF
3429 IF (mixed_cdft%dlb_control%send_work) THEN
3430 ind = count(mixed_cdft%dlb_control%recv_work_repl)
3431 DO j = 1, SIZE(mixed_cdft%dest_list)
3432 IF (touched(j)) THEN
3433 ALLOCATE (sbuff(j)%iv(4 + 3*nsend_proc(j)))
3434 sbuff(j)%iv(1:4) = [bo_conf(1, 2), bo_conf(2, 2), bo_conf(1, 3), bo_conf(2, 3)]
3435 offset = 5
3436 DO i = 1, SIZE(mixed_cdft%dlb_control%target_list, 2)
3437 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 mixed_cdft%dlb_control%target_list(4 + 2*j - 1, i)]
3441 offset = offset + 3
3442 END IF
3443 END DO
3444 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 request=req_total(ind + ispecial))
3448 END DO
3449 ind = ind + my_special_work
3450 END IF
3451 END DO
3452 END IF
3453 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3454 ALLOCATE (mixed_cdft%dlb_control%recv_info(SIZE(mixed_cdft%source_list)))
3455 ALLOCATE (mixed_cdft%dlb_control%recvbuff(SIZE(mixed_cdft%source_list)))
3456 ind = 1
3457 DO j = 1, SIZE(mixed_cdft%source_list)
3458 NULLIFY (mixed_cdft%dlb_control%recv_info(j)%target_list, &
3459 mixed_cdft%dlb_control%recvbuff(j)%buffs)
3460 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3461 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 request=req_total(ind))
3465 ind = ind + 1
3466 END IF
3467 END DO
3468 END IF
3469 IF (ASSOCIATED(req_total)) THEN
3470 CALL mp_waitall(req_total)
3471 DEALLOCATE (req_total)
3472 END IF
3473 IF (any(mask_send)) THEN
3474 ALLOCATE (tmp(SIZE(mixed_cdft%dest_list) - count(mask_send)), &
3475 tmp_bo(2, SIZE(mixed_cdft%dest_list) - count(mask_send)))
3476 i = 1
3477 DO j = 1, SIZE(mixed_cdft%dest_list)
3478 IF (.NOT. mask_send(j)) THEN
3479 tmp(i) = mixed_cdft%dest_list(j)
3480 tmp_bo(1:2, i) = mixed_cdft%dest_list_bo(1:2, j)
3481 i = i + 1
3482 END IF
3483 END DO
3484 DEALLOCATE (mixed_cdft%dest_list, mixed_cdft%dest_list_bo)
3485 ALLOCATE (mixed_cdft%dest_list(SIZE(tmp)), mixed_cdft%dest_list_bo(2, SIZE(tmp)))
3486 mixed_cdft%dest_list = tmp
3487 mixed_cdft%dest_list_bo = tmp_bo
3488 DEALLOCATE (tmp, tmp_bo)
3489 END IF
3490 IF (any(mask_recv)) THEN
3491 ALLOCATE (tmp(SIZE(mixed_cdft%source_list) - count(mask_recv)), &
3492 tmp_bo(4, SIZE(mixed_cdft%source_list) - count(mask_recv)))
3493 i = 1
3494 DO j = 1, SIZE(mixed_cdft%source_list)
3495 IF (.NOT. mask_recv(j)) THEN
3496 tmp(i) = mixed_cdft%source_list(j)
3497 tmp_bo(1:4, i) = mixed_cdft%source_list_bo(1:4, j)
3498 i = i + 1
3499 END IF
3500 END DO
3501 DEALLOCATE (mixed_cdft%source_list, mixed_cdft%source_list_bo)
3502 ALLOCATE (mixed_cdft%source_list(SIZE(tmp)), mixed_cdft%source_list_bo(4, SIZE(tmp)))
3503 mixed_cdft%source_list = tmp
3504 mixed_cdft%source_list_bo = tmp_bo
3505 DEALLOCATE (tmp, tmp_bo)
3506 END IF
3507 DEALLOCATE (mask_recv, mask_send)
3508 DEALLOCATE (nsend_proc, nrecv)
3509 IF (mixed_cdft%dlb_control%send_work) THEN
3510 DO j = 1, SIZE(mixed_cdft%dest_list)
3511 IF (touched(j)) DEALLOCATE (sbuff(j)%iv)
3512 END DO
3513 IF (ASSOCIATED(touched)) DEALLOCATE (touched)
3514 END IF
3515 END IF
3516 DEALLOCATE (sbuff)
3517 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
3518 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3519 CALL timestop(handle)
3520
3521 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 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 INTEGER, DIMENSION(:), POINTER :: work, work_dlb
3563 INTEGER, DIMENSION(:, :), POINTER :: nsent
3564 LOGICAL :: completed_recv, should_communicate
3565 LOGICAL, ALLOCATABLE, DIMENSION(:) :: skip_me
3566 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, distances, ds_dr_i, &
3571 ds_dr_j
3572 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_const_dr, d_sum_pm_dr, &
3573 distance_vecs, dp_i_dri
3574 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 REAL(kind=dp), DIMENSION(:), POINTER :: cutoffs
3579 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: cavity, weight
3580 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 TYPE(mp_request_type), DIMENSION(:), POINTER :: req_recv, req_total
3587 TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_send
3588 TYPE(particle_list_type), POINTER :: particles
3589 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 logger => cp_get_default_logger()
3594 NULLIFY (work, req_recv, req_send, work_dlb, nsent, cutoffs, cavity, &
3595 weight, gradients, cell, subsys_mix, force_env_qs, &
3596 particle_set, particles, auxbas_pw_pool, force_env_section, &
3597 print_section, cdft_control)
3598 CALL timeset(routinen, handle)
3599 nforce_eval = SIZE(force_env%sub_force_env)
3600 CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
3601 print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
3602 iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
3603 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 cell=cell)
3607 CALL cp_subsys_get(subsys=subsys_mix, &
3608 particles=particles, &
3609 particle_set=particle_set)
3610 ELSE
3611 DO iforce_eval = 1, nforce_eval
3612 IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) cycle
3613 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 cell=cell)
3618 CALL cp_subsys_get(subsys=subsys_mix, &
3619 particles=particles, &
3620 particle_set=particle_set)
3621 END IF
3622 natom = SIZE(particles%els)
3623 cdft_control => mixed_cdft%cdft_control
3624 CALL pw_env_get(pw_env=mixed_cdft%pw_env, auxbas_pw_pool=auxbas_pw_pool)
3625 np = auxbas_pw_pool%pw_grid%npts
3626 dr = auxbas_pw_pool%pw_grid%dr
3627 shift = -real(modulo(np, 2), dp)*dr/2.0_dp
3628 ALLOCATE (cell_functions(natom), skip_me(natom))
3629 IF (store_vectors) THEN
3630 ALLOCATE (distances(natom))
3631 ALLOCATE (distance_vecs(3, natom))
3632 END IF
3633 IF (in_memory) THEN
3634 ALLOCATE (ds_dr_j(3))
3635 ALLOCATE (ds_dr_i(3))
3636 ALLOCATE (d_sum_pm_dr(3, natom))
3637 ALLOCATE (d_sum_const_dr(3, natom))
3638 ALLOCATE (dp_i_drj(3, natom, natom))
3639 ALLOCATE (dp_i_dri(3, natom))
3640 th = 1.0e-8_dp
3641 END IF
3642 IF (mixed_cdft%dlb) THEN
3643 ALLOCATE (work(force_env%para_env%num_pe), work_dlb(force_env%para_env%num_pe))
3644 work = 0
3645 work_dlb = 0
3646 END IF
3647 my_work = 1
3648 my_special_work = 1
3649 ! Load balancing: allocate storage for receiving buffers and post recv requests
3650 IF (mixed_cdft%dlb) THEN
3651 IF (mixed_cdft%dlb_control%recv_work) THEN
3652 my_work = 2
3653 IF (.NOT. mixed_cdft%is_special) THEN
3654 ALLOCATE (req_send(2, 3))
3655 ELSE
3656 ALLOCATE (req_send(2, 3*SIZE(mixed_cdft%dlb_control%sendbuff)))
3657 END IF
3658 END IF
3659 IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
3660 IF (.NOT. mixed_cdft%is_special) THEN
3661 offset_repl = 0
3662 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 SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3665 offset_repl = 3*SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2)
3666 ELSE IF (mixed_cdft%dlb_control%recv_work_repl(1)) THEN
3667 ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(1)%target_list, 2))))
3668 ELSE
3669 ALLOCATE (req_recv(3*(SIZE(mixed_cdft%dlb_control%recv_info(2)%target_list, 2))))
3670 END IF
3671 ELSE
3672 nbuffs = 0
3673 offset_repl = 1
3674 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3675 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3676 nbuffs = nbuffs + (SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3
3677 END IF
3678 END DO
3679 ALLOCATE (req_recv(3*nbuffs))
3680 END IF
3681 DO j = 1, SIZE(mixed_cdft%dlb_control%recv_work_repl)
3682 IF (mixed_cdft%dlb_control%recv_work_repl(j)) THEN
3683 IF (.NOT. mixed_cdft%is_special) THEN
3684 offset = 0
3685 index = j + (j/2)
3686 ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)%buffs(SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)))
3687 DO i = 1, SIZE(mixed_cdft%dlb_control%recv_info(j)%target_list, 2)
3688 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 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 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 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 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 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 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 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 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 tag=mixed_cdft%dlb_control%recv_info(j)%target_list(3, i) + 2)
3754 offset = offset + mixed_cdft%dlb_control%recv_info(j)%target_list(2, i)
3755 END DO
3756 DEALLOCATE (mixed_cdft%dlb_control%recv_info(j)%matrix_info)
3757 ELSE
3758 ALLOCATE (mixed_cdft%dlb_control%recvbuff(j)% &
3759 buffs((SIZE(mixed_cdft%dlb_control%recv_info(j)%matrix_info) - 4)/3))
3760 index = 6
3761 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 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 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 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 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 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 request=req_recv(offset_repl + 2), tag=3)
3792 index = index + 3
3793 offset_repl = offset_repl + 3
3794 END DO
3795 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 cutoffs => cdft_control%becke_control%cutoffs
3802 should_communicate = .false.
3803 DO i = 1, 3
3804 cell_v(i) = cell%hmat(i, i)
3805 END DO
3806 DO iwork = my_work, 1, -1
3807 IF (iwork == 2) THEN
3808 IF (.NOT. mixed_cdft%is_special) THEN
3809 cavity => mixed_cdft%dlb_control%cavity
3810 weight => mixed_cdft%dlb_control%weight
3811 gradients => mixed_cdft%dlb_control%gradients
3812 ALLOCATE (completed(2, 3), nsent(2, 3))
3813 ELSE
3814 my_special_work = SIZE(mixed_cdft%dlb_control%sendbuff)
3815 ALLOCATE (completed(2, 3*my_special_work), nsent(2, 3*my_special_work))
3816 END IF
3817 completed = .false.
3818 nsent = 0
3819 ELSE
3820 IF (.NOT. mixed_cdft%is_special) THEN
3821 weight => mixed_cdft%weight
3822 cavity => mixed_cdft%cavity
3823 gradients => cdft_control%group(1)%gradients
3824 ELSE
3825 my_special_work = SIZE(mixed_cdft%dest_list)
3826 END IF
3827 END IF
3828 DO ispecial = 1, my_special_work
3829 nwork = 0
3830 IF (mixed_cdft%is_special) THEN
3831 IF (iwork == 1) THEN
3832 weight => mixed_cdft%sendbuff(ispecial)%weight
3833 cavity => mixed_cdft%sendbuff(ispecial)%cavity
3834 gradients => mixed_cdft%sendbuff(ispecial)%gradients
3835 ELSE
3836 weight => mixed_cdft%dlb_control%sendbuff(ispecial)%weight
3837 cavity => mixed_cdft%dlb_control%sendbuff(ispecial)%cavity
3838 gradients => mixed_cdft%dlb_control%sendbuff(ispecial)%gradients
3839 END IF
3840 END IF
3841 DO k = lbound(weight, 1), ubound(weight, 1)
3842 IF (mixed_cdft%dlb .AND. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3843 IF (mixed_cdft%dlb_control%send_work) THEN
3844 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 DO j = lbound(weight, 2), ubound(weight, 2)
3851 IF (mixed_cdft%dlb .AND. .NOT. mixed_cdft%is_pencil .AND. .NOT. mixed_cdft%is_special) THEN
3852 IF (mixed_cdft%dlb_control%send_work) THEN
3853 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 IF (should_communicate) THEN
3861 DO icomm = 1, SIZE(nsent, 2)
3862 DO jcomm = 1, SIZE(nsent, 1)
3863 IF (nsent(jcomm, icomm) == 1) cycle
3864 completed(jcomm, icomm) = req_send(jcomm, icomm)%test()
3865 IF (completed(jcomm, icomm)) THEN
3866 nsent(jcomm, icomm) = nsent(jcomm, icomm) + 1
3867 nsent_total = nsent_total + 1
3868 IF (nsent_total == SIZE(nsent, 1)*SIZE(nsent, 2)) should_communicate = .false.
3869 END IF
3870 IF (all(completed(:, icomm))) THEN
3871 IF (modulo(icomm, 3) == 1) THEN
3872 IF (.NOT. mixed_cdft%is_special) THEN
3873 DEALLOCATE (mixed_cdft%dlb_control%cavity)
3874 ELSE
3875 DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%cavity)
3876 END IF
3877 ELSE IF (modulo(icomm, 3) == 2) THEN
3878 IF (.NOT. mixed_cdft%is_special) THEN
3879 DEALLOCATE (mixed_cdft%dlb_control%weight)
3880 ELSE
3881 DEALLOCATE (mixed_cdft%dlb_control%sendbuff((icomm - 1)/3 + 1)%weight)
3882 END IF
3883 ELSE
3884 IF (.NOT. mixed_cdft%is_special) THEN
3885 DEALLOCATE (mixed_cdft%dlb_control%gradients)
3886 ELSE
3887 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 IF (ASSOCIATED(req_recv)) THEN
3896 completed_recv = mp_testall(req_recv)
3897 END IF
3898 !
3899 DO i = lbound(weight, 3), ubound(weight, 3)
3900 IF (cdft_control%becke_control%cavity_confine) THEN
3901 IF (cavity(k, j, i) < cdft_control%becke_control%eps_cavity) cycle
3902 END IF
3903 grid_p(1) = k*dr(1) + shift(1)
3904 grid_p(2) = j*dr(2) + shift(2)
3905 grid_p(3) = i*dr(3) + shift(3)
3906 nskipped = 0
3907 cell_functions = 1.0_dp
3908 skip_me = .false.
3909 IF (store_vectors) distances = 0.0_dp
3910 IF (in_memory) THEN
3911 d_sum_pm_dr = 0.0_dp
3912 d_sum_const_dr = 0.0_dp
3913 dp_i_dri = 0.0_dp
3914 END IF
3915 DO iatom = 1, natom
3916 IF (skip_me(iatom)) THEN
3917 cell_functions(iatom) = 0.0_dp
3918 IF (cdft_control%becke_control%should_skip) THEN
3919 IF (is_constraint(iatom)) nskipped = nskipped + 1
3920 IF (nskipped == cdft_control%natoms) THEN
3921 IF (in_memory) THEN
3922 IF (cdft_control%becke_control%cavity_confine) THEN
3923 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 IF (store_vectors) THEN
3932 IF (distances(iatom) == 0.0_dp) THEN
3933 r = position_vecs(:, iatom)
3934 dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
3935 dist1 = norm2(dist_vec)
3936 distance_vecs(:, iatom) = dist_vec
3937 distances(iatom) = dist1
3938 ELSE
3939 dist_vec = distance_vecs(:, iatom)
3940 dist1 = distances(iatom)
3941 END IF
3942 ELSE
3943 r = particle_set(iatom)%r
3944 DO ip = 1, 3
3945 r(ip) = modulo(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3946 END DO
3947 dist_vec = (r - grid_p) - anint((r - grid_p)/cell_v)*cell_v
3948 dist1 = norm2(dist_vec)
3949 END IF
3950 IF (dist1 <= cutoffs(iatom)) THEN
3951 IF (in_memory) THEN
3952 IF (dist1 <= th) dist1 = th
3953 dr_i_dr(:) = dist_vec(:)/dist1
3954 END IF
3955 DO jatom = 1, natom
3956 IF (jatom /= iatom) THEN
3957 IF (jatom < iatom) THEN
3958 IF (.NOT. skip_me(jatom)) cycle
3959 END IF
3960 IF (store_vectors) THEN
3961 IF (distances(jatom) == 0.0_dp) THEN
3962 r1 = position_vecs(:, jatom)
3963 dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
3964 dist2 = norm2(dist_vec)
3965 distance_vecs(:, jatom) = dist_vec
3966 distances(jatom) = dist2
3967 ELSE
3968 dist_vec = distance_vecs(:, jatom)
3969 dist2 = distances(jatom)
3970 END IF
3971 ELSE
3972 r1 = particle_set(jatom)%r
3973 DO ip = 1, 3
3974 r1(ip) = modulo(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp
3975 END DO
3976 dist_vec = (r1 - grid_p) - anint((r1 - grid_p)/cell_v)*cell_v
3977 dist2 = norm2(dist_vec)
3978 END IF
3979 IF (in_memory) THEN
3980 IF (store_vectors) THEN
3981 dr1_r2 = pair_dist_vecs(:, iatom, jatom)
3982 ELSE
3983 dr1_r2 = (r - r1) - anint((r - r1)/cell_v)*cell_v
3984 END IF
3985 IF (dist2 <= th) dist2 = th
3986 tmp_const = (r12(iatom, jatom)**3)
3987 dr_ij_dr(:) = dr1_r2(:)/tmp_const
3988 !derivativ w.r.t. Rj
3989 dr_j_dr = dist_vec(:)/dist2
3990 dmy_dr_j(:) = -(dr_j_dr(:)/r12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:))
3991 !derivativ w.r.t. Ri
3992 dmy_dr_i(:) = dr_i_dr(:)/r12(iatom, jatom) - (dist1 - dist2)*dr_ij_dr(:)
3993 END IF
3994 my1 = (dist1 - dist2)/r12(iatom, jatom)
3995 IF (cdft_control%becke_control%adjust) THEN
3996 my1_homo = my1
3997 my1 = my1 + &
3998 cdft_control%becke_control%aij(iatom, jatom)*(1.0_dp - my1**2)
3999 END IF
4000 myexp = 1.5_dp*my1 - 0.5_dp*my1**3
4001 IF (in_memory) THEN
4002 dmyexp = 1.5_dp - 1.5_dp*my1**2
4003 tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* &
4004 (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2))
4005
4006 ds_dr_i(:) = -0.5_dp*tmp_const*dmy_dr_i(:)
4007 ds_dr_j(:) = -0.5_dp*tmp_const*dmy_dr_j(:)
4008 IF (cdft_control%becke_control%adjust) THEN
4009 tmp_const = 1.0_dp - 2.0_dp*my1_homo*cdft_control%becke_control%aij(iatom, jatom)
4010 ds_dr_i(:) = ds_dr_i(:)*tmp_const
4011 ds_dr_j(:) = ds_dr_j(:)*tmp_const
4012 END IF
4013 END IF
4014 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4015 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3
4016 tmp_const = 0.5_dp*(1.0_dp - myexp)
4017 cell_functions(iatom) = cell_functions(iatom)*tmp_const
4018 IF (in_memory) THEN
4019 IF (abs(tmp_const) <= th) tmp_const = tmp_const + th
4020 dp_i_dri(:, iatom) = dp_i_dri(:, iatom) + ds_dr_i(:)/tmp_const
4021 dp_i_drj(:, iatom, jatom) = ds_dr_j(:)/tmp_const
4022 END IF
4023
4024 IF (dist2 <= cutoffs(jatom)) THEN
4025 tmp_const = 0.5_dp*(1.0_dp + myexp)
4026 cell_functions(jatom) = cell_functions(jatom)*tmp_const
4027 IF (in_memory) THEN
4028 IF (abs(tmp_const) <= th) tmp_const = tmp_const + th
4029 dp_i_drj(:, jatom, iatom) = -ds_dr_i(:)/tmp_const
4030 dp_i_dri(:, jatom) = dp_i_dri(:, jatom) - ds_dr_j(:)/tmp_const
4031 END IF
4032 ELSE
4033 skip_me(jatom) = .true.
4034 END IF
4035 END IF
4036 END DO
4037 IF (in_memory) THEN
4038 dp_i_dri(:, iatom) = cell_functions(iatom)*dp_i_dri(:, iatom)
4039 d_sum_pm_dr(:, iatom) = d_sum_pm_dr(:, iatom) + dp_i_dri(:, iatom)
4040 IF (is_constraint(iatom)) THEN
4041 d_sum_const_dr(:, iatom) = d_sum_const_dr(:, iatom) + dp_i_dri(:, iatom)* &
4042 coefficients(iatom)
4043 END IF
4044 DO jatom = 1, natom
4045 IF (jatom /= iatom) THEN
4046 IF (jatom < iatom) THEN
4047 IF (.NOT. skip_me(jatom)) THEN
4048 dp_i_drj(:, iatom, jatom) = cell_functions(iatom)*dp_i_drj(:, iatom, jatom)
4049 d_sum_pm_dr(:, jatom) = d_sum_pm_dr(:, jatom) + dp_i_drj(:, iatom, jatom)
4050 IF (is_constraint(iatom)) THEN
4051 d_sum_const_dr(:, jatom) = d_sum_const_dr(:, jatom) + &
4052 dp_i_drj(:, iatom, jatom)* &
4053 coefficients(iatom)
4054 END IF
4055 cycle
4056 END IF
4057 END IF
4058 dp_i_drj(:, iatom, jatom) = cell_functions(iatom)*dp_i_drj(:, iatom, jatom)
4059 d_sum_pm_dr(:, jatom) = d_sum_pm_dr(:, jatom) + dp_i_drj(:, iatom, jatom)
4060 IF (is_constraint(iatom)) THEN
4061 d_sum_const_dr(:, jatom) = d_sum_const_dr(:, jatom) + dp_i_drj(:, iatom, jatom)* &
4062 coefficients(iatom)
4063 END IF
4064 END IF
4065 END DO
4066 END IF
4067 ELSE
4068 cell_functions(iatom) = 0.0_dp
4069 skip_me(iatom) = .true.
4070 IF (cdft_control%becke_control%should_skip) THEN
4071 IF (is_constraint(iatom)) nskipped = nskipped + 1
4072 IF (nskipped == cdft_control%natoms) THEN
4073 IF (in_memory) THEN
4074 IF (cdft_control%becke_control%cavity_confine) THEN
4075 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 IF (nskipped == cdft_control%natoms) cycle
4084 sum_cell_f_constr = 0.0_dp
4085 DO ip = 1, cdft_control%natoms
4086 sum_cell_f_constr = sum_cell_f_constr + cell_functions(catom(ip))* &
4087 cdft_control%group(1)%coeff(ip)
4088 END DO
4089 sum_cell_f_all = 0.0_dp
4090 nwork = nwork + 1
4091 DO ip = 1, natom
4092 sum_cell_f_all = sum_cell_f_all + cell_functions(ip)
4093 END DO
4094 IF (in_memory) THEN
4095 DO iatom = 1, natom
4096 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 d_sum_pm_dr(:, iatom)/(sum_cell_f_all**2)
4100 END IF
4101 END DO
4102 END IF
4103 IF (abs(sum_cell_f_all) > eps_sum_cell_f_all) THEN
4104 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 IF (iwork == 2) THEN
4111 IF (.NOT. mixed_cdft%is_special) THEN
4112 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 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 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 tag=mixed_cdft%dlb_control%dest_tags_repl(i) + 2)
4125 END DO
4126 should_communicate = .true.
4127 nsent_total = 0
4128 ELSE
4129 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 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 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 request=req_send(i, 3*(ispecial - 1) + 3), tag=3)
4139 END DO
4140 IF (ispecial == my_special_work) THEN
4141 should_communicate = .true.
4142 nsent_total = 0
4143 END IF
4144 END IF
4145 work(mixed_cdft%dlb_control%my_source + 1) = work(mixed_cdft%dlb_control%my_source + 1) + nwork
4146 work_dlb(force_env%para_env%mepos + 1) = work_dlb(force_env%para_env%mepos + 1) + nwork
4147 ELSE
4148 IF (mixed_cdft%dlb) work(force_env%para_env%mepos + 1) = work(force_env%para_env%mepos + 1) + nwork
4149 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 IF (mixed_cdft%dlb) THEN
4155 IF (mixed_cdft%dlb_control%recv_work .AND. &
4156 any(mixed_cdft%dlb_control%recv_work_repl)) THEN
4157 ALLOCATE (req_total(SIZE(req_recv) + SIZE(req_send, 1)*SIZE(req_send, 2)))
4158 index = SIZE(req_recv)
4159 req_total(1:index) = req_recv
4160 DO i = 1, SIZE(req_send, 2)
4161 DO j = 1, SIZE(req_send, 1)
4162 index = index + 1
4163 req_total(index) = req_send(j, i)
4164 END DO
4165 END DO
4166 CALL mp_waitall(req_total)
4167 DEALLOCATE (req_total)
4168 IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4169 DEALLOCATE (mixed_cdft%dlb_control%cavity)
4170 END IF
4171 IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4172 DEALLOCATE (mixed_cdft%dlb_control%weight)
4173 END IF
4174 IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4175 DEALLOCATE (mixed_cdft%dlb_control%gradients)
4176 END IF
4177 IF (mixed_cdft%is_special) THEN
4178 DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4179 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4180 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4181 END IF
4182 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4183 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4184 END IF
4185 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4186 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4187 END IF
4188 END DO
4189 DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4190 END IF
4191 DEALLOCATE (req_send, req_recv)
4192 ELSE IF (mixed_cdft%dlb_control%recv_work) THEN
4193 IF (should_communicate) THEN
4194 CALL mp_waitall(req_send)
4195 END IF
4196 IF (ASSOCIATED(mixed_cdft%dlb_control%cavity)) THEN
4197 DEALLOCATE (mixed_cdft%dlb_control%cavity)
4198 END IF
4199 IF (ASSOCIATED(mixed_cdft%dlb_control%weight)) THEN
4200 DEALLOCATE (mixed_cdft%dlb_control%weight)
4201 END IF
4202 IF (ASSOCIATED(mixed_cdft%dlb_control%gradients)) THEN
4203 DEALLOCATE (mixed_cdft%dlb_control%gradients)
4204 END IF
4205 IF (mixed_cdft%is_special) THEN
4206 DO j = 1, SIZE(mixed_cdft%dlb_control%sendbuff)
4207 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%cavity)) THEN
4208 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%cavity)
4209 END IF
4210 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%weight)) THEN
4211 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%weight)
4212 END IF
4213 IF (ASSOCIATED(mixed_cdft%dlb_control%sendbuff(j)%gradients)) THEN
4214 DEALLOCATE (mixed_cdft%dlb_control%sendbuff(j)%gradients)
4215 END IF
4216 END DO
4217 DEALLOCATE (mixed_cdft%dlb_control%sendbuff)
4218 END IF
4219 DEALLOCATE (req_send)
4220 ELSE IF (any(mixed_cdft%dlb_control%recv_work_repl)) THEN
4221 CALL mp_waitall(req_recv)
4222 DEALLOCATE (req_recv)
4223 END IF
4224 END IF
4225 IF (mixed_cdft%dlb) THEN
4226 CALL force_env%para_env%sum(work)
4227 CALL force_env%para_env%sum(work_dlb)
4228 IF (.NOT. ASSOCIATED(mixed_cdft%dlb_control%prediction_error)) THEN
4229 ALLOCATE (mixed_cdft%dlb_control%prediction_error(force_env%para_env%num_pe))
4230 END IF
4231 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 DEALLOCATE (work, work_dlb, mixed_cdft%dlb_control%expected_work)
4239 END IF
4240 NULLIFY (gradients, weight, cavity)
4241 IF (ALLOCATED(coefficients)) THEN
4242 DEALLOCATE (coefficients)
4243 END IF
4244 IF (in_memory) THEN
4245 DEALLOCATE (ds_dr_j)
4246 DEALLOCATE (ds_dr_i)
4247 DEALLOCATE (d_sum_pm_dr)
4248 DEALLOCATE (d_sum_const_dr)
4249 DEALLOCATE (dp_i_drj)
4250 DEALLOCATE (dp_i_dri)
4251 NULLIFY (gradients)
4252 IF (store_vectors) THEN
4253 DEALLOCATE (pair_dist_vecs)
4254 END IF
4255 END IF
4256 NULLIFY (cutoffs)
4257 IF (ALLOCATED(is_constraint)) THEN
4258 DEALLOCATE (is_constraint)
4259 END IF
4260 DEALLOCATE (catom)
4261 DEALLOCATE (r12)
4262 DEALLOCATE (cell_functions)
4263 DEALLOCATE (skip_me)
4264 IF (ALLOCATED(completed)) THEN
4265 DEALLOCATE (completed)
4266 END IF
4267 IF (ASSOCIATED(nsent)) THEN
4268 DEALLOCATE (nsent)
4269 END IF
4270 IF (store_vectors) THEN
4271 DEALLOCATE (distances)
4272 DEALLOCATE (distance_vecs)
4273 DEALLOCATE (position_vecs)
4274 END IF
4275 IF (ASSOCIATED(req_send)) THEN
4276 DEALLOCATE (req_send)
4277 END IF
4278 IF (ASSOCIATED(req_recv)) THEN
4279 DEALLOCATE (req_recv)
4280 END IF
4281 CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
4282 "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
4283 CALL timestop(handle)
4284
4285 END SUBROUTINE mixed_becke_constraint_low
4286
4287END MODULE mixed_cdft_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cell_types::pbc
Definition cell_types.F:104
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
message_passing::mp_testall
Definition message_passing.F:890
message_passing::mp_waitall
Definition message_passing.F:886
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_zero
Definition pw_methods.F:83
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius_very_extended
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition ao_util.F:208
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_syevd
subroutine, public cp_dbcsr_syevd(matrix, eigenvectors, eigenvalues, para_env, blacs_env)
...
Definition cp_dbcsr_diag.F:67
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:799
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1690
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1532
cp_fm_basic_linalg::cp_fm_invert
subroutine, public cp_fm_invert(matrix_a, matrix_inverse, det_a, eps_svd, eigval)
Inverts a cp_fm_type matrix, optionally returning the determinant of the input matrix.
Definition cp_fm_basic_linalg.F:1790
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_write_formatted
subroutine, public cp_fm_write_formatted(fm, unit, header, value_format)
Write out a full matrix in plain text.
Definition cp_fm_types.F:2236
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::use_qmmm
integer, parameter, public use_qmmm
Definition force_env_types.F:83
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:83
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition force_env_types.F:83
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition grid_api.F:27
grid_api::collocate_pgf_product
subroutine, public collocate_pgf_product(la_max, zeta, la_min, lb_max, zetb, lb_min, ra, rab, scale, pab, o1, o2, rsgrid, ga_gb_function, radius, use_subpatch, subpatch_pattern)
low level collocation of primitive gaussian functions
Definition grid_api.F:116
hirshfeld_methods
Calculate Hirshfeld charges and related functions.
Definition hirshfeld_methods.F:14
hirshfeld_methods::create_shape_function
subroutine, public create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, radius, radii_list)
creates kind specific shape functions for Hirshfeld charges
Definition hirshfeld_methods.F:179
hirshfeld_types
The types needed for the calculation of Hirshfeld charges and related functions.
Definition hirshfeld_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::mixed_cdft_serial
integer, parameter, public mixed_cdft_serial
Definition input_constants.F:826
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition input_constants.F:835
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:820
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:820
input_constants::becke_cutoff_element
integer, parameter, public becke_cutoff_element
Definition input_constants.F:812
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:820
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:772
input_constants::mixed_cdft_parallel
integer, parameter, public mixed_cdft_parallel
Definition input_constants.F:826
input_constants::becke_cutoff_global
integer, parameter, public becke_cutoff_global
Definition input_constants.F:812
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:820
input_constants::mixed_cdft_parallel_nobuild
integer, parameter, public mixed_cdft_parallel_nobuild
Definition input_constants.F:826
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1017
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mixed_cdft_methods
Methods for mixed CDFT calculations.
Definition mixed_cdft_methods.F:14
mixed_cdft_methods::mixed_cdft_calculate_coupling
subroutine, public mixed_cdft_calculate_coupling(force_env)
Driver routine to calculate the electronic coupling(s) between CDFT states.
Definition mixed_cdft_methods.F:1130
mixed_cdft_methods::mixed_cdft_build_weight
subroutine, public mixed_cdft_build_weight(force_env, calculate_forces, iforce_eval)
Driver routine to handle the build of CDFT weight/gradient in parallel and serial modes.
Definition mixed_cdft_methods.F:420
mixed_cdft_methods::mixed_cdft_init
subroutine, public mixed_cdft_init(force_env, calculate_forces)
Initialize a mixed CDFT calculation.
Definition mixed_cdft_methods.F:154
mixed_cdft_types
Types for mixed CDFT calculations.
Definition mixed_cdft_types.F:14
mixed_cdft_types::mixed_cdft_result_type_set
subroutine, public mixed_cdft_result_type_set(results, lowdin, wfn, nonortho, metric, rotation, h, s, wad, wda, w_diagonal, energy, strength, s_minushalf)
Updates arrays within the mixed CDFT result container.
Definition mixed_cdft_types.F:510
mixed_cdft_types::mixed_cdft_type_create
subroutine, public mixed_cdft_type_create(cdft_control)
inits the given mixed_cdft_type
Definition mixed_cdft_types.F:290
mixed_cdft_types::mixed_cdft_work_type_release
subroutine, public mixed_cdft_work_type_release(matrix)
Releases arrays within the mixed CDFT work matrix container.
Definition mixed_cdft_types.F:631
mixed_cdft_utils
Utility subroutines for mixed CDFT calculations.
Definition mixed_cdft_utils.F:14
mixed_cdft_utils::map_permutation_to_states
subroutine, public map_permutation_to_states(n, ipermutation, i, j)
Given the size of a symmetric matrix and a permutation index, returns indices (i, j) of the off-diago...
Definition mixed_cdft_utils.F:1522
mixed_cdft_utils::mixed_cdft_init_structures
subroutine, public mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
Initialize all the structures needed for a mixed CDFT calculation.
Definition mixed_cdft_utils.F:657
mixed_cdft_utils::mixed_cdft_transfer_settings
subroutine, public mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
Transfer settings to mixed_cdft.
Definition mixed_cdft_utils.F:554
mixed_cdft_utils::mixed_cdft_diagonalize_blocks
subroutine, public mixed_cdft_diagonalize_blocks(blocks, h_block, s_block, eigenvalues)
Diagonalizes each of the matrix blocks.
Definition mixed_cdft_utils.F:1770
mixed_cdft_utils::mixed_cdft_print_couplings
subroutine, public mixed_cdft_print_couplings(force_env)
Routine to print out the electronic coupling(s) between CDFT states.
Definition mixed_cdft_utils.F:1379
mixed_cdft_utils::mixed_cdft_read_block_diag
subroutine, public mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
Read input section related to block diagonalization of the mixed CDFT Hamiltonian matrix.
Definition mixed_cdft_utils.F:1626
mixed_cdft_utils::mixed_cdft_redistribute_arrays
subroutine, public mixed_cdft_redistribute_arrays(force_env)
Redistribute arrays needed for an ET coupling calculation from individual CDFT states to the mixed CD...
Definition mixed_cdft_utils.F:1071
mixed_cdft_utils::mixed_cdft_assemble_block_diag
subroutine, public mixed_cdft_assemble_block_diag(mixed_cdft, blocks, h_block, eigenvalues, n, iounit)
Assembles the new block diagonalized mixed CDFT Hamiltonian and overlap matrices.
Definition mixed_cdft_utils.F:1831
mixed_cdft_utils::hfun_zero
subroutine, public hfun_zero(fun, th, just_zero, bounds, work)
Determine confinement bounds along confinement dir (hardcoded to be z) and determine the number of no...
Definition mixed_cdft_utils.F:1556
mixed_cdft_utils::mixed_cdft_get_blocks
subroutine, public mixed_cdft_get_blocks(mixed_cdft, blocks, h_block, s_block)
Assembles the matrix blocks from the mixed CDFT Hamiltonian.
Definition mixed_cdft_utils.F:1722
mixed_cdft_utils::mixed_cdft_release_work
subroutine, public mixed_cdft_release_work(force_env)
Release storage reserved for mixed CDFT matrices.
Definition mixed_cdft_utils.F:1499
mixed_cdft_utils::mixed_cdft_parse_settings
subroutine, public mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, settings, natom)
Parse settings for mixed cdft calculation and check their consistency.
Definition mixed_cdft_utils.F:141
mixed_environment_types
Definition mixed_environment_types.F:11
mixed_environment_types::get_mixed_env
subroutine, public get_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell, cell_ref, mixed_energy, para_env, sub_para_env, subsys, input, results, cdft_control)
Get the MIXED environment.
Definition mixed_environment_types.F:130
mixed_environment_types::set_mixed_env
subroutine, public set_mixed_env(mixed_env, atomic_kind_set, particle_set, local_particles, local_molecules, molecule_kind_set, molecule_set, cell_ref, mixed_energy, subsys, input, sub_para_env, cdft_control)
Set the MIXED environment.
Definition mixed_environment_types.F:198
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
qs_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_mo_io
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
qs_mo_io::wfn_restart_file_name
subroutine, public wfn_restart_file_name(filename, exist, section, logger, kp, xas, rtp)
...
Definition qs_mo_io.F:460
qs_mo_io::read_mo_set_from_restart
subroutine, public read_mo_set_from_restart(mo_array, qs_kind_set, particle_set, para_env, id_nr, multiplicity, dft_section, natom_mismatch, cdft, out_unit)
...
Definition qs_mo_io.F:526
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_methods::make_basis_simple
subroutine, public make_basis_simple(vmatrix, ncol)
given a set of vectors, return an orthogonal (C^T C == 1) set spanning the same space (notice,...
Definition qs_mo_methods.F:356
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition qs_mo_methods.F:89
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::set_mo_set
subroutine, public set_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, uniform_occupation, kts, mu, flexible_electron_count)
Set the components of a MO set data structure.
Definition qs_mo_types.F:502
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:224
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:394
realspace_grid_types
Definition realspace_grid_types.F:18
realspace_grid_types::transfer_rs2pw
subroutine, public transfer_rs2pw(rs, pw)
...
Definition realspace_grid_types.F:652
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1913
util
All kind of helpful little routines.
Definition util.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_1d_i_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:221
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
hirshfeld_types::hirshfeld_type
quantities needed for a Hirshfeld based partitioning of real space
Definition hirshfeld_types.F:36
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_request_type
Definition message_passing.F:628
mixed_cdft_types::mixed_cdft_settings_type
Container for constraint settings to check consistency of force_evals.
Definition mixed_cdft_types.F:243
mixed_cdft_types::mixed_cdft_type
Main mixed CDFT control type.
Definition mixed_cdft_types.F:204
mixed_environment_types::mixed_environment_type
Definition mixed_environment_types.F:64
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:219
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136