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