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