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