(git:ed6f26b)
Loading...
Searching...
No Matches
mixed_cdft_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods for mixed CDFT calculations
10!> \par History
11!> Separated CDFT routines from mixed_environment_utils
12!> \author Nico Holmberg [01.2017]
13! **************************************************************************************************
14MODULE mixed_cdft_methods
15 USE ao_util, ONLY: exp_radius_very_extended
16 USE atomic_kind_types, ONLY: atomic_kind_type,&
17 get_atomic_kind
18 USE cell_types, ONLY: cell_type,&
19 pbc
20 USE cp_array_utils, ONLY: cp_1d_i_p_type,&
21 cp_1d_r_p_type,&
22 cp_2d_r_p_type
23 USE cp_control_types, ONLY: dft_control_type
24 USE cp_dbcsr_api, ONLY: &
25 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_init_p, dbcsr_p_type, dbcsr_release, &
26 dbcsr_release_p, dbcsr_scale, dbcsr_type
27 USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
28 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
29 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
30 cp_fm_invert,&
31 cp_fm_transpose
32 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 cp_fm_struct_release,&
34 cp_fm_struct_type
35 USE cp_fm_types, ONLY: cp_fm_create,&
36 cp_fm_get_info,&
37 cp_fm_release,&
38 cp_fm_set_all,&
39 cp_fm_to_fm,&
40 cp_fm_type,&
41 cp_fm_write_formatted
42 USE cp_log_handling, ONLY: cp_get_default_logger,&
43 cp_logger_type,&
44 cp_to_string
45 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
46 cp_print_key_unit_nr
47 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
48 USE cp_subsys_types, ONLY: cp_subsys_get,&
49 cp_subsys_type
50 USE cp_units, ONLY: cp_unit_from_cp2k
51 USE force_env_types, ONLY: force_env_get,&
52 force_env_type,&
53 use_qmmm,&
54 use_qmmmx,&
55 use_qs_force
56 USE grid_api, ONLY: grid_func_ab,&
57 collocate_pgf_product
58 USE hirshfeld_methods, ONLY: create_shape_function
59 USE hirshfeld_types, ONLY: hirshfeld_type
60 USE input_constants, ONLY: &
61 becke_cutoff_element, becke_cutoff_global, cdft_alpha_constraint, cdft_beta_constraint, &
62 cdft_charge_constraint, cdft_magnetization_constraint, mix_cdft, mixed_cdft_parallel, &
63 mixed_cdft_parallel_nobuild, mixed_cdft_serial, outer_scf_becke_constraint
64 USE input_section_types, ONLY: section_get_lval,&
65 section_vals_get,&
66 section_vals_get_subs_vals,&
67 section_vals_type,&
68 section_vals_val_get
69 USE kinds, ONLY: default_path_length,&
70 dp
71 USE machine, ONLY: m_walltime
72 USE mathlib, ONLY: diamat_all
73 USE memory_utilities, ONLY: reallocate
74 USE message_passing, ONLY: mp_request_type,&
75 mp_testall,&
76 mp_waitall
77 USE mixed_cdft_types, ONLY: mixed_cdft_result_type_set,&
78 mixed_cdft_settings_type,&
79 mixed_cdft_type,&
80 mixed_cdft_type_create,&
81 mixed_cdft_work_type_release
82 USE mixed_cdft_utils, ONLY: &
83 hfun_zero, map_permutation_to_states, mixed_cdft_assemble_block_diag, &
84 mixed_cdft_diagonalize_blocks, mixed_cdft_get_blocks, mixed_cdft_init_structures, &
85 mixed_cdft_parse_settings, mixed_cdft_print_couplings, mixed_cdft_read_block_diag, &
86 mixed_cdft_redistribute_arrays, mixed_cdft_release_work, mixed_cdft_transfer_settings
87 USE mixed_environment_types, ONLY: get_mixed_env,&
88 mixed_environment_type,&
89 set_mixed_env
90 USE parallel_gemm_api, ONLY: parallel_gemm
91 USE particle_list_types, ONLY: particle_list_type
92 USE particle_types, ONLY: particle_type
93 USE pw_env_types, ONLY: pw_env_get,&
94 pw_env_type
95 USE pw_methods, ONLY: pw_copy,&
96 pw_scale,&
97 pw_zero
98 USE pw_pool_types, ONLY: pw_pool_type
99 USE qs_cdft_types, ONLY: cdft_control_type
100 USE qs_energy_types, ONLY: qs_energy_type
101 USE qs_environment_types, ONLY: get_qs_env,&
102 qs_environment_type
103 USE qs_kind_types, ONLY: qs_kind_type
104 USE qs_mo_io, ONLY: read_mo_set_from_restart,&
105 wfn_restart_file_name
106 USE qs_mo_methods, ONLY: make_basis_simple,&
107 make_basis_sm
108 USE qs_mo_types, ONLY: allocate_mo_set,&
109 deallocate_mo_set,&
110 mo_set_type,&
111 set_mo_set
112 USE realspace_grid_types, ONLY: realspace_grid_type,&
113 rs_grid_zero,&
114 transfer_rs2pw
115 USE util, ONLY: sort
116#include "./base/base_uses.f90"
117
118 IMPLICIT NONE
119
120 PRIVATE
121
122 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_methods'
123 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
124
125 PUBLIC :: mixed_cdft_init, &
126 mixed_cdft_build_weight, &
127 mixed_cdft_calculate_coupling
128
129CONTAINS
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
4210END MODULE mixed_cdft_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:120
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
memory_utilities::reallocate
Definition memory_utilities.F:27
message_passing::mp_testall
Definition message_passing.F:875
message_passing::mp_waitall
Definition message_passing.F:871
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_scale
Definition pw_methods.F:93
pw_methods::pw_zero
Definition pw_methods.F:82
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:220
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:205
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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:552
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:1885
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:1695
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:2019
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
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:1009
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:528
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:164
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:2250
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
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:83
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_types::use_qmmmx
integer, parameter, public use_qmmmx
Definition force_env_types.F:83
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition force_env_types.F:83
grid_api
Fortran API for the grid package, which is written in C.
Definition grid_api.F:12
grid_api::grid_func_ab
integer, parameter, public grid_func_ab
Definition grid_api.F: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:778
input_constants::mix_cdft
integer, parameter, public mix_cdft
Definition input_constants.F:787
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:772
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:772
input_constants::becke_cutoff_element
integer, parameter, public becke_cutoff_element
Definition input_constants.F:764
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:772
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:724
input_constants::mixed_cdft_parallel
integer, parameter, public mixed_cdft_parallel
Definition input_constants.F:778
input_constants::becke_cutoff_global
integer, parameter, public becke_cutoff_global
Definition input_constants.F:764
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:772
input_constants::mixed_cdft_parallel_nobuild
integer, parameter, public mixed_cdft_parallel_nobuild
Definition input_constants.F:778
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:731
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:704
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:1047
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1012
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:147
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:373
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_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_type_create
subroutine, public mixed_cdft_type_create(cdft_control)
inits the given mixed_cdft_type
Definition mixed_cdft_types.F:290
mixed_cdft_types::mixed_cdft_work_type_release
subroutine, public mixed_cdft_work_type_release(matrix)
Releases arrays within the mixed CDFT work matrix container.
Definition mixed_cdft_types.F: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:1484
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:631
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:533
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:1726
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:1342
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:1587
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:1039
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:1787
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:1517
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:1679
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:1461
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:140
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_pp, 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, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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:652
realspace_grid_types::rs_grid_zero
subroutine, public rs_grid_zero(rs)
Initialize grid to zero.
Definition realspace_grid_types.F:1905
util
All kind of helpful little routines.
Definition util.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_array_utils::cp_1d_i_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:221
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
hirshfeld_types::hirshfeld_type
quantities needed for a Hirshfeld based partitioning of real space
Definition hirshfeld_types.F:36
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_request_type
Definition message_passing.F:623
mixed_cdft_types::mixed_cdft_settings_type
Container for constraint settings to check consistency of force_evals.
Definition mixed_cdft_types.F:243
mixed_cdft_types::mixed_cdft_type
Main mixed CDFT control type.
Definition mixed_cdft_types.F:204
mixed_environment_types::mixed_environment_type
Definition mixed_environment_types.F:64
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:219
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:136