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