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