1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Subroutines for building CDFT constraints 8!> \par History 9!> separated from et_coupling [03.2017] 10!> \author Nico Holmberg [03.2017] 11! ************************************************************************************************** 12MODULE qs_cdft_methods 13 USE ao_util, ONLY: exp_radius_very_extended 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind,& 16 get_atomic_kind_set 17 USE cell_types, ONLY: cell_type,& 18 pbc 19 USE cp_control_types, ONLY: dft_control_type 20 USE cp_log_handling, ONLY: cp_get_default_logger,& 21 cp_logger_type 22 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 23 cp_print_key_unit_nr 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw 26 USE cube_utils, ONLY: cube_info_type 27 USE grid_api, ONLY: GRID_FUNC_AB,& 28 collocate_pgf_product 29 USE hirshfeld_types, ONLY: get_hirshfeld_info,& 30 hirshfeld_type,& 31 set_hirshfeld_info 32 USE input_constants, ONLY: cdft_alpha_constraint,& 33 cdft_beta_constraint,& 34 cdft_charge_constraint,& 35 cdft_magnetization_constraint,& 36 outer_scf_becke_constraint,& 37 outer_scf_hirshfeld_constraint 38 USE input_section_types, ONLY: section_vals_get_subs_vals,& 39 section_vals_type 40 USE kahan_sum, ONLY: accurate_dot_product,& 41 accurate_sum 42 USE kinds, ONLY: dp,& 43 int_8 44 USE message_passing, ONLY: mp_sum 45 USE particle_types, ONLY: particle_type 46 USE pw_env_types, ONLY: pw_env_get,& 47 pw_env_type 48 USE pw_methods, ONLY: pw_axpy,& 49 pw_copy,& 50 pw_integrate_function 51 USE pw_pool_types, ONLY: pw_pool_create_pw,& 52 pw_pool_give_back_pw,& 53 pw_pool_type 54 USE pw_types, ONLY: REALDATA3D,& 55 REALSPACE,& 56 pw_p_type 57 USE qs_cdft_types, ONLY: becke_constraint_type,& 58 cdft_control_type,& 59 cdft_group_type,& 60 hirshfeld_constraint_type 61 USE qs_cdft_utils, ONLY: becke_constraint_init,& 62 cdft_constraint_print,& 63 hfun_scale,& 64 hirshfeld_constraint_init 65 USE qs_energy_types, ONLY: qs_energy_type 66 USE qs_environment_types, ONLY: get_qs_env,& 67 qs_environment_type 68 USE qs_force_types, ONLY: qs_force_type 69 USE qs_kind_types, ONLY: get_qs_kind,& 70 qs_kind_type 71 USE qs_rho0_types, ONLY: get_rho0_mpole,& 72 mpole_rho_atom,& 73 rho0_mpole_type 74 USE qs_rho_types, ONLY: qs_rho_get,& 75 qs_rho_type 76 USE qs_subsys_types, ONLY: qs_subsys_get,& 77 qs_subsys_type 78 USE realspace_grid_types, ONLY: & 79 realspace_grid_desc_type, realspace_grid_p_type, realspace_grid_type, rs2pw, & 80 rs_grid_create, rs_grid_release, rs_grid_retain, rs_grid_zero, rs_pw_transfer 81#include "./base/base_uses.f90" 82 83 IMPLICIT NONE 84 85 PRIVATE 86 87 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_methods' 88 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 89 90! *** Public subroutines *** 91 92 PUBLIC :: becke_constraint, hirshfeld_constraint 93 94CONTAINS 95 96! ************************************************************************************************** 97!> \brief Driver routine for calculating a Becke constraint 98!> \param qs_env the qs_env where to build the constraint 99!> \param calc_pot if the potential needs to be recalculated or just integrated 100!> \param calculate_forces logical if potential has to be calculated or only_energy 101!> \par History 102!> Created 01.2007 [fschiff] 103!> Extended functionality 12/15-12/16 [Nico Holmberg] 104! ************************************************************************************************** 105 SUBROUTINE becke_constraint(qs_env, calc_pot, calculate_forces) 106 TYPE(qs_environment_type), POINTER :: qs_env 107 LOGICAL :: calc_pot, calculate_forces 108 109 CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint', & 110 routineP = moduleN//':'//routineN 111 112 INTEGER :: handle 113 TYPE(cdft_control_type), POINTER :: cdft_control 114 TYPE(dft_control_type), POINTER :: dft_control 115 116 CALL timeset(routineN, handle) 117 CALL get_qs_env(qs_env, dft_control=dft_control) 118 cdft_control => dft_control%qs_control%cdft_control 119 IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_becke_constraint) THEN 120 IF (calc_pot) THEN 121 ! Initialize the Becke constraint environment 122 CALL becke_constraint_init(qs_env) 123 ! Calculate the Becke weight function and possibly the gradients 124 CALL becke_constraint_low(qs_env) 125 END IF 126 ! Integrate the Becke constraint 127 CALL cdft_constraint_integrate(qs_env) 128 ! Calculate forces 129 IF (calculate_forces) CALL becke_constraint_force(qs_env) 130 END IF 131 CALL timestop(handle) 132 133 END SUBROUTINE becke_constraint 134 135! ************************************************************************************************** 136!> \brief Low level routine to build a Becke weight function and its gradients 137!> \param qs_env the qs_env where to build the constraint 138!> \param just_gradients optional logical which determines if only the gradients should be calculated 139!> \par History 140!> Created 03.2017 [Nico Holmberg] 141! ************************************************************************************************** 142 SUBROUTINE becke_constraint_low(qs_env, just_gradients) 143 TYPE(qs_environment_type), POINTER :: qs_env 144 LOGICAL, OPTIONAL :: just_gradients 145 146 CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_low', & 147 routineP = moduleN//':'//routineN 148 149 INTEGER :: handle, i, iatom, igroup, ind(3), ip, j, & 150 jatom, jp, k, natom, np(3), nskipped 151 INTEGER, ALLOCATABLE, DIMENSION(:) :: catom 152 INTEGER, DIMENSION(2, 3) :: bo, bo_conf 153 LOGICAL :: in_memory, my_just_gradients 154 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint, skip_me 155 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: atom_in_group 156 REAL(kind=dp) :: dist1, dist2, dmyexp, dvol, eps_cavity, & 157 my1, my1_homo, myexp, sum_cell_f_all, & 158 th, tmp_const 159 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cell_functions, ds_dR_i, ds_dR_j, & 160 sum_cell_f_group 161 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_sum_Pm_dR, dP_i_dRi 162 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dP_i_dRj 163 REAL(kind=dp), DIMENSION(3) :: cell_v, dist_vec, dmy_dR_i, dmy_dR_j, & 164 dr, dr1_r2, dr_i_dR, dr_ij_dR, & 165 dr_j_dR, grid_p, r, r1, shift 166 REAL(KIND=dp), DIMENSION(:), POINTER :: cutoffs 167 TYPE(becke_constraint_type), POINTER :: becke_control 168 TYPE(cdft_control_type), POINTER :: cdft_control 169 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 170 TYPE(cell_type), POINTER :: cell 171 TYPE(dft_control_type), POINTER :: dft_control 172 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 173 TYPE(pw_p_type), DIMENSION(:), POINTER :: charge 174 175 NULLIFY (cutoffs, cell, dft_control, particle_set, group, charge, cdft_control) 176 CALL timeset(routineN, handle) 177 ! Get simulation environment 178 CALL get_qs_env(qs_env, & 179 cell=cell, & 180 particle_set=particle_set, & 181 natom=natom, & 182 dft_control=dft_control) 183 cdft_control => dft_control%qs_control%cdft_control 184 becke_control => cdft_control%becke_control 185 group => cdft_control%group 186 cutoffs => becke_control%cutoffs 187 IF (cdft_control%atomic_charges) & 188 charge => cdft_control%charge 189 in_memory = .FALSE. 190 IF (cdft_control%save_pot) THEN 191 in_memory = becke_control%in_memory 192 END IF 193 eps_cavity = becke_control%eps_cavity 194 ! Decide if only gradients need to be calculated 195 my_just_gradients = .FALSE. 196 IF (PRESENT(just_gradients)) my_just_gradients = just_gradients 197 IF (my_just_gradients) THEN 198 in_memory = .TRUE. 199 ! Pairwise distances need to be recalculated 200 IF (becke_control%vector_buffer%store_vectors) THEN 201 ALLOCATE (becke_control%vector_buffer%distances(natom)) 202 ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom)) 203 IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom)) 204 ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom)) 205 END IF 206 ALLOCATE (becke_control%vector_buffer%R12(natom, natom)) 207 DO i = 1, 3 208 cell_v(i) = cell%hmat(i, i) 209 END DO 210 DO iatom = 1, natom - 1 211 DO jatom = iatom + 1, natom 212 r = particle_set(iatom)%r 213 r1 = particle_set(jatom)%r 214 DO i = 1, 3 215 r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp 216 r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp 217 END DO 218 dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v 219 IF (becke_control%vector_buffer%store_vectors) THEN 220 becke_control%vector_buffer%position_vecs(:, iatom) = r(:) 221 IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:) 222 IF (in_memory) THEN 223 becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:) 224 becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:) 225 END IF 226 END IF 227 becke_control%vector_buffer%R12(iatom, jatom) = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 228 becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom) 229 END DO 230 END DO 231 END IF 232 ALLOCATE (catom(cdft_control%natoms)) 233 IF (cdft_control%save_pot .OR. & 234 becke_control%cavity_confine .OR. & 235 becke_control%should_skip) THEN 236 ALLOCATE (is_constraint(natom)) 237 is_constraint = .FALSE. 238 END IF 239 ! This boolean is needed to prevent calculation of atom pairs ji when the pair ij has 240 ! already been calculated (data for pair ji is set using symmetry) 241 ! With gradient precomputation, symmetry exploited for both weight function and gradients 242 ALLOCATE (skip_me(natom)) 243 DO i = 1, cdft_control%natoms 244 catom(i) = cdft_control%atoms(i) 245 ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges 246 ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct 247 IF (cdft_control%save_pot .OR. & 248 becke_control%cavity_confine .OR. & 249 becke_control%should_skip) & 250 is_constraint(catom(i)) = .TRUE. 251 END DO 252 bo = group(1)%weight%pw%pw_grid%bounds_local 253 dvol = group(1)%weight%pw%pw_grid%dvol 254 dr = group(1)%weight%pw%pw_grid%dr 255 np = group(1)%weight%pw%pw_grid%npts 256 shift = -REAL(MODULO(np, 2), dp)*dr/2.0_dp 257 DO i = 1, 3 258 cell_v(i) = cell%hmat(i, i) 259 END DO 260 ! If requested, allocate storage for gradients 261 IF (in_memory) THEN 262 bo_conf = bo 263 ! With confinement active, we dont need to store gradients outside 264 ! the confinement bounds since they vanish for all particles 265 IF (becke_control%cavity_confine) THEN 266 bo_conf(1, 3) = becke_control%confine_bounds(1) 267 bo_conf(2, 3) = becke_control%confine_bounds(2) 268 END IF 269 ALLOCATE (atom_in_group(SIZE(group), natom)) 270 atom_in_group = .FALSE. 271 DO igroup = 1, SIZE(group) 272 ALLOCATE (group(igroup)%gradients(3*natom, bo_conf(1, 1):bo_conf(2, 1), & 273 bo_conf(1, 2):bo_conf(2, 2), & 274 bo_conf(1, 3):bo_conf(2, 3))) 275 group(igroup)%gradients = 0.0_dp 276 ALLOCATE (group(igroup)%d_sum_const_dR(3, natom)) 277 group(igroup)%d_sum_const_dR = 0.0_dp 278 DO ip = 1, SIZE(group(igroup)%atoms) 279 atom_in_group(igroup, group(igroup)%atoms(ip)) = .TRUE. 280 END DO 281 END DO 282 END IF 283 ! Allocate remaining work 284 ALLOCATE (sum_cell_f_group(SIZE(group))) 285 ALLOCATE (cell_functions(natom)) 286 IF (in_memory) THEN 287 ALLOCATE (ds_dR_j(3)) 288 ALLOCATE (ds_dR_i(3)) 289 ALLOCATE (d_sum_Pm_dR(3, natom)) 290 ALLOCATE (dP_i_dRj(3, natom, natom)) 291 ALLOCATE (dP_i_dRi(3, natom)) 292 th = 1.0e-8_dp 293 END IF 294 ! Build constraint 295 DO k = bo(1, 1), bo(2, 1) 296 DO j = bo(1, 2), bo(2, 2) 297 DO i = bo(1, 3), bo(2, 3) 298 ! If the grid point is too far from all constraint atoms and cavity confinement is active, 299 ! we can skip this grid point as it does not contribute to the weight or gradients 300 IF (becke_control%cavity_confine) THEN 301 IF (becke_control%cavity%pw%cr3d(k, j, i) < eps_cavity) CYCLE 302 END IF 303 ind = (/k, j, i/) 304 grid_p(1) = k*dr(1) + shift(1) 305 grid_p(2) = j*dr(2) + shift(2) 306 grid_p(3) = i*dr(3) + shift(3) 307 nskipped = 0 308 cell_functions = 1.0_dp 309 skip_me = .FALSE. 310 IF (becke_control%vector_buffer%store_vectors) becke_control%vector_buffer%distances = 0.0_dp 311 IF (in_memory) THEN 312 d_sum_Pm_dR = 0.0_dp 313 DO igroup = 1, SIZE(group) 314 group(igroup)%d_sum_const_dR = 0.0_dp 315 END DO 316 dP_i_dRi = 0.0_dp 317 END IF 318 ! Iterate over all atoms in the system 319 DO iatom = 1, natom 320 IF (skip_me(iatom)) THEN 321 cell_functions(iatom) = 0.0_dp 322 IF (becke_control%should_skip) THEN 323 IF (is_constraint(iatom)) nskipped = nskipped + 1 324 IF (nskipped == cdft_control%natoms) THEN 325 IF (in_memory) THEN 326 IF (becke_control%cavity_confine) THEN 327 becke_control%cavity%pw%cr3d(k, j, i) = 0.0_dp 328 END IF 329 END IF 330 EXIT 331 END IF 332 END IF 333 CYCLE 334 END IF 335 IF (becke_control%vector_buffer%store_vectors) THEN 336 IF (becke_control%vector_buffer%distances(iatom) .EQ. 0.0_dp) THEN 337 r = becke_control%vector_buffer%position_vecs(:, iatom) 338 dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v 339 dist1 = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 340 becke_control%vector_buffer%distance_vecs(:, iatom) = dist_vec 341 becke_control%vector_buffer%distances(iatom) = dist1 342 ELSE 343 dist_vec = becke_control%vector_buffer%distance_vecs(:, iatom) 344 dist1 = becke_control%vector_buffer%distances(iatom) 345 END IF 346 ELSE 347 r = particle_set(iatom)%r 348 DO ip = 1, 3 349 r(ip) = MODULO(r(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp 350 END DO 351 dist_vec = (r - grid_p) - ANINT((r - grid_p)/cell_v)*cell_v 352 dist1 = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 353 END IF 354 IF (dist1 .LE. cutoffs(iatom)) THEN 355 IF (in_memory) THEN 356 IF (dist1 .LE. th) dist1 = th 357 dr_i_dR(:) = dist_vec(:)/dist1 358 END IF 359 DO jatom = 1, natom 360 IF (jatom .NE. iatom) THEN 361 ! Using pairwise symmetry, execute block only for such j<i 362 ! that have previously not been looped over 363 ! Note that if skip_me(jatom) = .TRUE., this means that the outer 364 ! loop over iatom skipped this index when iatom=jatom, but we still 365 ! need to compute the pair for iatom>jatom 366 IF (jatom < iatom) THEN 367 IF (.NOT. skip_me(jatom)) CYCLE 368 END IF 369 IF (becke_control%vector_buffer%store_vectors) THEN 370 IF (becke_control%vector_buffer%distances(jatom) .EQ. 0.0_dp) THEN 371 r1 = becke_control%vector_buffer%position_vecs(:, jatom) 372 dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v 373 dist2 = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 374 becke_control%vector_buffer%distance_vecs(:, jatom) = dist_vec 375 becke_control%vector_buffer%distances(jatom) = dist2 376 ELSE 377 dist_vec = becke_control%vector_buffer%distance_vecs(:, jatom) 378 dist2 = becke_control%vector_buffer%distances(jatom) 379 END IF 380 ELSE 381 r1 = particle_set(jatom)%r 382 DO ip = 1, 3 383 r1(ip) = MODULO(r1(ip), cell%hmat(ip, ip)) - cell%hmat(ip, ip)/2._dp 384 END DO 385 dist_vec = (r1 - grid_p) - ANINT((r1 - grid_p)/cell_v)*cell_v 386 dist2 = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 387 END IF 388 IF (in_memory) THEN 389 IF (becke_control%vector_buffer%store_vectors) THEN 390 dr1_r2 = becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) 391 ELSE 392 dr1_r2 = (r - r1) - ANINT((r - r1)/cell_v)*cell_v 393 END IF 394 IF (dist2 .LE. th) dist2 = th 395 tmp_const = (becke_control%vector_buffer%R12(iatom, jatom)**3) 396 dr_ij_dR(:) = dr1_r2(:)/tmp_const 397 !derivative w.r.t. Rj 398 dr_j_dR = dist_vec(:)/dist2 399 dmy_dR_j(:) = -(dr_j_dR(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:)) 400 !derivative w.r.t. Ri 401 dmy_dR_i(:) = dr_i_dR(:)/becke_control%vector_buffer%R12(iatom, jatom) - (dist1 - dist2)*dr_ij_dR(:) 402 END IF 403 ! myij 404 my1 = (dist1 - dist2)/becke_control%vector_buffer%R12(iatom, jatom) 405 IF (becke_control%adjust) THEN 406 my1_homo = my1 ! Homonuclear quantity needed for gradient 407 my1 = my1 + becke_control%aij(iatom, jatom)*(1.0_dp - my1**2) 408 END IF 409 ! f(myij) 410 myexp = 1.5_dp*my1 - 0.5_dp*my1**3 411 IF (in_memory) THEN 412 dmyexp = 1.5_dp - 1.5_dp*my1**2 413 tmp_const = (1.5_dp**2)*dmyexp*(1 - myexp**2)* & 414 (1.0_dp - ((1.5_dp*myexp - 0.5_dp*(myexp**3))**2)) 415 ! d s(myij)/d R_i 416 ds_dR_i(:) = -0.5_dp*tmp_const*dmy_dR_i(:) 417 ! d s(myij)/d R_j 418 ds_dR_j(:) = -0.5_dp*tmp_const*dmy_dR_j(:) 419 IF (becke_control%adjust) THEN 420 tmp_const = 1.0_dp - 2.0_dp*my1_homo* & 421 becke_control%aij(iatom, jatom) 422 ds_dR_i(:) = ds_dR_i(:)*tmp_const 423 ! tmp_const is same for both since aij=-aji and myij=-myji 424 ds_dR_j(:) = ds_dR_j(:)*tmp_const 425 END IF 426 END IF 427 ! s(myij) = f[f(f{myij})] 428 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3 429 myexp = 1.5_dp*myexp - 0.5_dp*myexp**3 430 tmp_const = 0.5_dp*(1.0_dp - myexp) 431 cell_functions(iatom) = cell_functions(iatom)*tmp_const 432 IF (in_memory) THEN 433 IF (ABS(tmp_const) .LE. th) tmp_const = tmp_const + th 434 ! P_i independent part of dP_i/dR_i 435 dP_i_dRi(:, iatom) = dP_i_dRi(:, iatom) + ds_dR_i(:)/tmp_const 436 ! P_i independent part of dP_i/dR_j 437 dP_i_dRj(:, iatom, jatom) = ds_dR_j(:)/tmp_const 438 END IF 439 440 IF (dist2 .LE. cutoffs(jatom)) THEN 441 tmp_const = 0.5_dp*(1.0_dp + myexp) ! s(myji) 442 cell_functions(jatom) = cell_functions(jatom)*tmp_const 443 IF (in_memory) THEN 444 IF (ABS(tmp_const) .LE. th) tmp_const = tmp_const + th 445 ! P_j independent part of dP_j/dR_i 446 ! d s(myji)/d R_i = -d s(myij)/d R_i 447 dP_i_dRj(:, jatom, iatom) = -ds_dR_i(:)/tmp_const 448 ! P_j independent part of dP_j/dR_j 449 ! d s(myji)/d R_j = -d s(myij)/d R_j 450 dP_i_dRi(:, jatom) = dP_i_dRi(:, jatom) - ds_dR_j(:)/tmp_const 451 END IF 452 ELSE 453 skip_me(jatom) = .TRUE. 454 END IF 455 END IF 456 END DO ! jatom 457 IF (in_memory) THEN 458 ! Final value of dP_i_dRi 459 dP_i_dRi(:, iatom) = cell_functions(iatom)*dP_i_dRi(:, iatom) 460 ! Update relevant sums with value 461 d_sum_Pm_dR(:, iatom) = d_sum_Pm_dR(:, iatom) + dP_i_dRi(:, iatom) 462 IF (is_constraint(iatom)) THEN 463 DO igroup = 1, SIZE(group) 464 IF (.NOT. atom_in_group(igroup, iatom)) CYCLE 465 DO jp = 1, SIZE(group(igroup)%atoms) 466 IF (iatom == group(igroup)%atoms(jp)) THEN 467 ip = jp 468 EXIT 469 END IF 470 END DO 471 group(igroup)%d_sum_const_dR(1:3, iatom) = group(igroup)%d_sum_const_dR(1:3, iatom) + & 472 group(igroup)%coeff(ip)*dP_i_dRi(:, iatom) 473 END DO 474 END IF 475 DO jatom = 1, natom 476 IF (jatom .NE. iatom) THEN 477 ! Final value of dP_i_dRj 478 dP_i_dRj(:, iatom, jatom) = cell_functions(iatom)*dP_i_dRj(:, iatom, jatom) 479 ! Update where needed 480 d_sum_Pm_dR(:, jatom) = d_sum_Pm_dR(:, jatom) + dP_i_dRj(:, iatom, jatom) 481 IF (is_constraint(iatom)) THEN 482 DO igroup = 1, SIZE(group) 483 IF (.NOT. atom_in_group(igroup, iatom)) CYCLE 484 ip = -1 485 DO jp = 1, SIZE(group(igroup)%atoms) 486 IF (iatom == group(igroup)%atoms(jp)) THEN 487 ip = jp 488 EXIT 489 END IF 490 END DO 491 group(igroup)%d_sum_const_dR(1:3, jatom) = group(igroup)%d_sum_const_dR(1:3, jatom) + & 492 group(igroup)%coeff(ip)* & 493 dP_i_dRj(:, iatom, jatom) 494 END DO 495 END IF 496 END IF 497 END DO 498 END IF 499 ELSE 500 cell_functions(iatom) = 0.0_dp 501 skip_me(iatom) = .TRUE. 502 IF (becke_control%should_skip) THEN 503 IF (is_constraint(iatom)) nskipped = nskipped + 1 504 IF (nskipped == cdft_control%natoms) THEN 505 IF (in_memory) THEN 506 IF (becke_control%cavity_confine) THEN 507 becke_control%cavity%pw%cr3d(k, j, i) = 0.0_dp 508 END IF 509 END IF 510 EXIT 511 END IF 512 END IF 513 END IF 514 END DO !iatom 515 IF (nskipped == cdft_control%natoms) CYCLE 516 ! Sum up cell functions 517 sum_cell_f_group = 0.0_dp 518 DO igroup = 1, SIZE(group) 519 DO ip = 1, SIZE(group(igroup)%atoms) 520 sum_cell_f_group(igroup) = sum_cell_f_group(igroup) + group(igroup)%coeff(ip)* & 521 cell_functions(group(igroup)%atoms(ip)) 522 END DO 523 END DO 524 sum_cell_f_all = 0.0_dp 525 DO ip = 1, natom 526 sum_cell_f_all = sum_cell_f_all + cell_functions(ip) 527 END DO 528 ! Gradients at (k,j,i) 529 IF (in_memory .AND. ABS(sum_cell_f_all) .GT. 0.0_dp) THEN 530 DO igroup = 1, SIZE(group) 531 DO iatom = 1, natom 532 group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i) = & 533 group(igroup)%d_sum_const_dR(1:3, iatom)/sum_cell_f_all - sum_cell_f_group(igroup)* & 534 d_sum_Pm_dR(1:3, iatom)/(sum_cell_f_all**2) 535 END DO 536 END DO 537 END IF 538 ! Weight function(s) at (k,j,i) 539 IF (.NOT. my_just_gradients .AND. ABS(sum_cell_f_all) .GT. 0.000001) THEN 540 DO igroup = 1, SIZE(group) 541 group(igroup)%weight%pw%cr3d(k, j, i) = sum_cell_f_group(igroup)/sum_cell_f_all 542 END DO 543 IF (cdft_control%atomic_charges) THEN 544 DO iatom = 1, cdft_control%natoms 545 charge(iatom)%pw%cr3d(k, j, i) = cell_functions(catom(iatom))/sum_cell_f_all 546 END DO 547 END IF 548 END IF 549 END DO 550 END DO 551 END DO 552 ! Release storage 553 IF (in_memory) THEN 554 DEALLOCATE (ds_dR_j) 555 DEALLOCATE (ds_dR_i) 556 DEALLOCATE (d_sum_Pm_dR) 557 DEALLOCATE (dP_i_dRj) 558 DEALLOCATE (dP_i_dRi) 559 DO igroup = 1, SIZE(group) 560 DEALLOCATE (group(igroup)%d_sum_const_dR) 561 END DO 562 DEALLOCATE (atom_in_group) 563 IF (becke_control%vector_buffer%store_vectors) THEN 564 DEALLOCATE (becke_control%vector_buffer%pair_dist_vecs) 565 END IF 566 END IF 567 NULLIFY (cutoffs) 568 IF (ALLOCATED(is_constraint)) & 569 DEALLOCATE (is_constraint) 570 DEALLOCATE (catom) 571 DEALLOCATE (cell_functions) 572 DEALLOCATE (skip_me) 573 DEALLOCATE (sum_cell_f_group) 574 DEALLOCATE (becke_control%vector_buffer%R12) 575 IF (becke_control%vector_buffer%store_vectors) THEN 576 DEALLOCATE (becke_control%vector_buffer%distances) 577 DEALLOCATE (becke_control%vector_buffer%distance_vecs) 578 DEALLOCATE (becke_control%vector_buffer%position_vecs) 579 END IF 580 CALL timestop(handle) 581 582 END SUBROUTINE becke_constraint_low 583 584! ************************************************************************************************** 585!> \brief Calculates atomic forces due to a Becke constraint 586!> \param qs_env the qs_env where to build the gradients 587!> \par History 588!> Created 01.2007 [fschiff] 589!> Extended functionality 12/15-12/16 [Nico Holmberg] 590! ************************************************************************************************** 591 SUBROUTINE becke_constraint_force(qs_env) 592 TYPE(qs_environment_type), POINTER :: qs_env 593 594 CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_force', & 595 routineP = moduleN//':'//routineN 596 597 INTEGER :: handle, i, iatom, igroup, ikind, ind(3), & 598 ispin, j, k, natom, nvar 599 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 600 INTEGER, DIMENSION(2, 3) :: bo 601 REAL(kind=dp) :: dvol, eps_cavity, sign, th 602 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: strength 603 REAL(KIND=dp), DIMENSION(:), POINTER :: cutoffs 604 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 605 TYPE(becke_constraint_type), POINTER :: becke_control 606 TYPE(cdft_control_type), POINTER :: cdft_control 607 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 608 TYPE(cell_type), POINTER :: cell 609 TYPE(cp_para_env_type), POINTER :: para_env 610 TYPE(dft_control_type), POINTER :: dft_control 611 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 612 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 613 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 614 TYPE(qs_rho_type), POINTER :: rho 615 616 CALL timeset(routineN, handle) 617 NULLIFY (atomic_kind_set, cell, para_env, dft_control, particle_set, & 618 rho, rho_r, force, cutoffs, becke_control, group) 619 620 CALL get_qs_env(qs_env, & 621 atomic_kind_set=atomic_kind_set, & 622 natom=natom, & 623 particle_set=particle_set, & 624 cell=cell, & 625 rho=rho, & 626 force=force, & 627 dft_control=dft_control, & 628 para_env=para_env) 629 CALL qs_rho_get(rho, rho_r=rho_r) 630 631 th = 1.0e-8_dp 632 cdft_control => dft_control%qs_control%cdft_control 633 becke_control => cdft_control%becke_control 634 group => cdft_control%group 635 nvar = SIZE(cdft_control%target) 636 ALLOCATE (strength(nvar)) 637 strength(:) = cdft_control%strength(:) 638 cutoffs => becke_control%cutoffs 639 eps_cavity = becke_control%eps_cavity 640 ALLOCATE (atom_of_kind(natom)) 641 ALLOCATE (kind_of(natom)) 642 643 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 644 atom_of_kind=atom_of_kind, & 645 kind_of=kind_of) 646 DO igroup = 1, SIZE(group) 647 ALLOCATE (group(igroup)%integrated(3, natom)) 648 group(igroup)%integrated = 0.0_dp 649 END DO 650 bo = group(1)%weight%pw%pw_grid%bounds_local 651 dvol = group(1)%weight%pw%pw_grid%dvol 652 653 IF (.NOT. becke_control%in_memory) & 654 ! Gradients were not precomputed => calculate them now 655 CALL becke_constraint_low(qs_env, just_gradients=.TRUE.) 656 ! Integrate gradients 657 IF (.NOT. ASSOCIATED(becke_control%cavity_mat)) THEN 658 ! No external control 659 DO k = bo(1, 1), bo(2, 1) 660 DO j = bo(1, 2), bo(2, 2) 661 DO i = bo(1, 3), bo(2, 3) 662 ! First check if this grid point should be skipped 663 IF (becke_control%cavity_confine) THEN 664 IF (becke_control%cavity%pw%cr3d(k, j, i) < eps_cavity) CYCLE 665 END IF 666 ind = (/k, j, i/) 667 DO igroup = 1, SIZE(group) 668 DO iatom = 1, natom 669 DO ispin = 1, dft_control%nspins 670 SELECT CASE (group(igroup)%constraint_type) 671 CASE (cdft_charge_constraint) 672 sign = 1.0_dp 673 CASE (cdft_magnetization_constraint) 674 IF (ispin == 1) THEN 675 sign = 1.0_dp 676 ELSE 677 sign = -1.0_dp 678 END IF 679 CASE (cdft_alpha_constraint) 680 sign = 1.0_dp 681 IF (ispin == 2) CYCLE 682 CASE (cdft_beta_constraint) 683 sign = 1.0_dp 684 IF (ispin == 1) CYCLE 685 CASE DEFAULT 686 CPABORT("Unknown constraint type.") 687 END SELECT 688 group(igroup)%integrated(:, iatom) = group(igroup)%integrated(:, iatom) + sign* & 689 group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i)* & 690 rho_r(ispin)%pw%cr3d(k, j, i)*dvol 691 END DO 692 END DO 693 END DO 694 END DO 695 END DO 696 END DO 697 ELSE 698 DO k = LBOUND(becke_control%cavity_mat, 1), UBOUND(becke_control%cavity_mat, 1) 699 DO j = LBOUND(becke_control%cavity_mat, 2), UBOUND(becke_control%cavity_mat, 2) 700 DO i = LBOUND(becke_control%cavity_mat, 3), UBOUND(becke_control%cavity_mat, 3) 701 ! First check if this grid point should be skipped 702 IF (becke_control%cavity_mat(k, j, i) < eps_cavity) CYCLE 703 DO igroup = 1, SIZE(group) 704 DO iatom = 1, natom 705 DO ispin = 1, dft_control%nspins 706 SELECT CASE (group(igroup)%constraint_type) 707 CASE (cdft_charge_constraint) 708 sign = 1.0_dp 709 CASE (cdft_magnetization_constraint) 710 IF (ispin == 1) THEN 711 sign = 1.0_dp 712 ELSE 713 sign = -1.0_dp 714 END IF 715 CASE (cdft_alpha_constraint) 716 sign = 1.0_dp 717 IF (ispin == 2) CYCLE 718 CASE (cdft_beta_constraint) 719 sign = 1.0_dp 720 IF (ispin == 1) CYCLE 721 CASE DEFAULT 722 CPABORT("Unknown constraint type.") 723 END SELECT 724 group(igroup)%integrated(:, iatom) = group(igroup)%integrated(:, iatom) + sign* & 725 group(igroup)%gradients(3*(iatom - 1) + 1:3*(iatom - 1) + 3, k, j, i)* & 726 rho_r(ispin)%pw%cr3d(k, j, i)*dvol 727 END DO 728 END DO 729 END DO 730 END DO 731 END DO 732 END DO 733 END IF 734 IF (.NOT. cdft_control%transfer_pot) THEN 735 DO igroup = 1, SIZE(group) 736 DEALLOCATE (group(igroup)%gradients) 737 END DO 738 END IF 739 DO igroup = 1, SIZE(group) 740 CALL mp_sum(group(igroup)%integrated, para_env%group) 741 END DO 742 ! Update force only on master process. Otherwise force due to constraint becomes multiplied 743 ! by the number of processes when the final force%rho_elec is constructed in qs_force 744 ! by mp_summing [the final integrated(:,:) is distributed on all processors] 745 IF (para_env%ionode) THEN 746 DO igroup = 1, SIZE(group) 747 DO iatom = 1, natom 748 ikind = kind_of(iatom) 749 i = atom_of_kind(iatom) 750 force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + group(igroup)%integrated(:, iatom)*strength(igroup) 751 END DO 752 END DO 753 END IF 754 DEALLOCATE (strength) 755 DEALLOCATE (atom_of_kind) 756 DEALLOCATE (kind_of) 757 DO igroup = 1, SIZE(group) 758 DEALLOCATE (group(igroup)%integrated) 759 END DO 760 NULLIFY (group) 761 762 CALL timestop(handle) 763 764 END SUBROUTINE becke_constraint_force 765 766! ************************************************************************************************** 767!> \brief Calculates the value of a CDFT constraint by integrating the product of the CDFT 768!> weight function and the realspace electron density 769!> \param qs_env the qs_env where to build the constraint 770!> \par History 771!> Created 3.2017 [Nico Holmberg] 772!> Generalized to any CDFT weight function 9.2018 [Nico Holmberg] 773! ************************************************************************************************** 774 SUBROUTINE cdft_constraint_integrate(qs_env) 775 TYPE(qs_environment_type), POINTER :: qs_env 776 777 CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_integrate', & 778 routineP = moduleN//':'//routineN 779 780 INTEGER :: handle, i, iatom, igroup, ikind, ivar, & 781 iw, jatom, natom, nvar 782 LOGICAL :: is_becke, paw_atom 783 REAL(kind=dp) :: dvol, eps_cavity, sign 784 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dE, strength, target_val 785 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge, gapw_offset 786 TYPE(becke_constraint_type), POINTER :: becke_control 787 TYPE(cdft_control_type), POINTER :: cdft_control 788 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 789 TYPE(cp_logger_type), POINTER :: logger 790 TYPE(cp_para_env_type), POINTER :: para_env 791 TYPE(dft_control_type), POINTER :: dft_control 792 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 793 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 794 TYPE(pw_p_type), DIMENSION(:), POINTER :: charge, rho_r 795 TYPE(qs_energy_type), POINTER :: energy 796 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 797 TYPE(qs_rho_type), POINTER :: rho 798 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 799 TYPE(section_vals_type), POINTER :: cdft_constraint_section 800 801 NULLIFY (para_env, dft_control, particle_set, rho_r, energy, rho, & 802 logger, cdft_constraint_section, qs_kind_set, mp_rho, & 803 rho0_mpole, group, charge) 804 CALL timeset(routineN, handle) 805 logger => cp_get_default_logger() 806 CALL get_qs_env(qs_env, & 807 particle_set=particle_set, & 808 rho=rho, & 809 natom=natom, & 810 dft_control=dft_control, & 811 para_env=para_env, & 812 qs_kind_set=qs_kind_set) 813 CALL qs_rho_get(rho, rho_r=rho_r) 814 CPASSERT(ASSOCIATED(qs_kind_set)) 815 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") 816 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog") 817 cdft_control => dft_control%qs_control%cdft_control 818 is_becke = (cdft_control%type == outer_scf_becke_constraint) 819 becke_control => cdft_control%becke_control 820 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) & 821 CPABORT("Becke control has not been allocated.") 822 group => cdft_control%group 823 ! Initialize 824 nvar = SIZE(cdft_control%target) 825 ALLOCATE (strength(nvar)) 826 ALLOCATE (target_val(nvar)) 827 ALLOCATE (dE(nvar)) 828 strength(:) = cdft_control%strength(:) 829 target_val(:) = cdft_control%target(:) 830 dE = 0.0_dp 831 dvol = group(1)%weight%pw%pw_grid%dvol 832 IF (cdft_control%atomic_charges) THEN 833 charge => cdft_control%charge 834 ALLOCATE (electronic_charge(cdft_control%natoms, dft_control%nspins)) 835 electronic_charge = 0.0_dp 836 END IF 837 ! Calculate value of constraint i.e. int ( rho(r) w(r) dr) 838 DO i = 1, dft_control%nspins 839 DO igroup = 1, SIZE(group) 840 SELECT CASE (group(igroup)%constraint_type) 841 CASE (cdft_charge_constraint) 842 sign = 1.0_dp 843 CASE (cdft_magnetization_constraint) 844 IF (i == 1) THEN 845 sign = 1.0_dp 846 ELSE 847 sign = -1.0_dp 848 END IF 849 CASE (cdft_alpha_constraint) 850 sign = 1.0_dp 851 IF (i == 2) CYCLE 852 CASE (cdft_beta_constraint) 853 sign = 1.0_dp 854 IF (i == 1) CYCLE 855 CASE DEFAULT 856 CPABORT("Unknown constraint type.") 857 END SELECT 858 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN 859 ! With external control, we can use cavity_mat as a mask to kahan sum 860 eps_cavity = becke_control%eps_cavity 861 IF (igroup /= 1) & 862 CALL cp_abort(__LOCATION__, & 863 "Multiple constraints not yet supported by parallel mixed calculations.") 864 dE(igroup) = dE(igroup) + sign*accurate_dot_product(group(igroup)%weight%pw%cr3d, rho_r(i)%pw%cr3d, & 865 becke_control%cavity_mat, eps_cavity)*dvol 866 ELSE 867 dE(igroup) = dE(igroup) + sign*accurate_sum(group(igroup)%weight%pw%cr3d*rho_r(i)%pw%cr3d)*dvol 868 END IF 869 END DO 870 IF (cdft_control%atomic_charges) THEN 871 DO iatom = 1, cdft_control%natoms 872 electronic_charge(iatom, i) = accurate_sum(charge(iatom)%pw%cr3d*rho_r(i)%pw%cr3d)*dvol 873 END DO 874 END IF 875 END DO 876 CALL get_qs_env(qs_env, energy=energy) 877 CALL mp_sum(dE, para_env%group) 878 IF (cdft_control%atomic_charges) THEN 879 CALL mp_sum(electronic_charge, para_env%group) 880 END IF 881 ! Use fragment densities as reference value (= Becke deformation density) 882 IF (cdft_control%fragment_density .AND. .NOT. cdft_control%fragments_integrated) THEN 883 CALL prepare_fragment_constraint(qs_env) 884 END IF 885 IF (dft_control%qs_control%gapw) THEN 886 ! GAPW: add core charges (rho_hard - rho_soft) 887 IF (cdft_control%fragment_density) & 888 CALL cp_abort(__LOCATION__, & 889 "Fragment constraints not yet compatible with GAPW.") 890 ALLOCATE (gapw_offset(nvar, dft_control%nspins)) 891 gapw_offset = 0.0_dp 892 CALL get_qs_env(qs_env, rho0_mpole=rho0_mpole) 893 CALL get_rho0_mpole(rho0_mpole, mp_rho=mp_rho) 894 DO i = 1, dft_control%nspins 895 DO igroup = 1, SIZE(group) 896 DO iatom = 1, SIZE(group(igroup)%atoms) 897 SELECT CASE (group(igroup)%constraint_type) 898 CASE (cdft_charge_constraint) 899 sign = 1.0_dp 900 CASE (cdft_magnetization_constraint) 901 IF (i == 1) THEN 902 sign = 1.0_dp 903 ELSE 904 sign = -1.0_dp 905 END IF 906 CASE (cdft_alpha_constraint) 907 sign = 1.0_dp 908 IF (i == 2) CYCLE 909 CASE (cdft_beta_constraint) 910 sign = 1.0_dp 911 IF (i == 1) CYCLE 912 CASE DEFAULT 913 CPABORT("Unknown constraint type.") 914 END SELECT 915 jatom = group(igroup)%atoms(iatom) 916 CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind) 917 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom) 918 IF (paw_atom) THEN 919 gapw_offset(igroup, i) = gapw_offset(igroup, i) + sign*group(igroup)%coeff(iatom)*mp_rho(jatom)%q0(i) 920 END IF 921 END DO 922 END DO 923 END DO 924 IF (cdft_control%atomic_charges) THEN 925 DO iatom = 1, cdft_control%natoms 926 jatom = cdft_control%atoms(iatom) 927 CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=ikind) 928 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom) 929 IF (paw_atom) THEN 930 DO i = 1, dft_control%nspins 931 electronic_charge(iatom, i) = electronic_charge(iatom, i) + mp_rho(jatom)%q0(i) 932 END DO 933 END IF 934 END DO 935 END IF 936 DO i = 1, dft_control%nspins 937 DO ivar = 1, nvar 938 dE(ivar) = dE(ivar) + gapw_offset(ivar, i) 939 END DO 940 END DO 941 DEALLOCATE (gapw_offset) 942 END IF 943 ! Update constraint value and energy 944 cdft_control%value(:) = dE(:) 945 energy%cdft = 0.0_dp 946 DO ivar = 1, nvar 947 energy%cdft = energy%cdft + (dE(ivar) - target_val(ivar))*strength(ivar) 948 END DO 949 ! Print constraint info and atomic CDFT charges 950 CALL cdft_constraint_print(qs_env, electronic_charge) 951 ! Deallocate tmp storage 952 DEALLOCATE (dE, strength, target_val) 953 IF (cdft_control%atomic_charges) DEALLOCATE (electronic_charge) 954 CALL cp_print_key_finished_output(iw, logger, cdft_constraint_section, "PROGRAM_RUN_INFO") 955 CALL timestop(handle) 956 957 END SUBROUTINE cdft_constraint_integrate 958 959! ************************************************************************************************** 960!> \brief Prepare CDFT fragment constraints. Fragment densities are read from cube files, multiplied 961!> by the CDFT weight functions and integrated over the realspace grid. 962!> \param qs_env the qs_env where to build the constraint 963!> \par History 964!> Created 9.2018 [Nico Holmberg] 965! ************************************************************************************************** 966 SUBROUTINE prepare_fragment_constraint(qs_env) 967 TYPE(qs_environment_type), POINTER :: qs_env 968 969 CHARACTER(len=*), PARAMETER :: routineN = 'prepare_fragment_constraint', & 970 routineP = moduleN//':'//routineN 971 972 INTEGER :: handle, i, iatom, igroup, natom, & 973 nelectron_total, nfrag_spins 974 LOGICAL :: is_becke, needs_spin_density 975 REAL(kind=dp) :: dvol, multiplier(2), nelectron_frag 976 TYPE(becke_constraint_type), POINTER :: becke_control 977 TYPE(cdft_control_type), POINTER :: cdft_control 978 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 979 TYPE(cp_logger_type), POINTER :: logger 980 TYPE(cp_para_env_type), POINTER :: para_env 981 TYPE(dft_control_type), POINTER :: dft_control 982 TYPE(pw_env_type), POINTER :: pw_env 983 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_frag 984 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 985 TYPE(qs_subsys_type), POINTER :: subsys 986 987 NULLIFY (para_env, dft_control, logger, subsys, pw_env, auxbas_pw_pool, group, rho_frag) 988 CALL timeset(routineN, handle) 989 logger => cp_get_default_logger() 990 CALL get_qs_env(qs_env, & 991 natom=natom, & 992 dft_control=dft_control, & 993 para_env=para_env) 994 995 cdft_control => dft_control%qs_control%cdft_control 996 is_becke = (cdft_control%type == outer_scf_becke_constraint) 997 becke_control => cdft_control%becke_control 998 IF (is_becke .AND. .NOT. ASSOCIATED(becke_control)) & 999 CPABORT("Becke control has not been allocated.") 1000 group => cdft_control%group 1001 dvol = group(1)%weight%pw%pw_grid%dvol 1002 ! Fragment densities are meaningful only for some calculation types 1003 IF (.NOT. qs_env%single_point_run) & 1004 CALL cp_abort(__LOCATION__, & 1005 "CDFT fragment constraints are only compatible with single "// & 1006 "point calculations (run_type ENERGY or ENERGY_FORCE).") 1007 IF (dft_control%qs_control%gapw) & 1008 CALL cp_abort(__LOCATION__, & 1009 "CDFT fragment constraint not compatible with GAPW.") 1010 needs_spin_density = .FALSE. 1011 multiplier = 1.0_dp 1012 nfrag_spins = 1 1013 DO igroup = 1, SIZE(group) 1014 SELECT CASE (group(igroup)%constraint_type) 1015 CASE (cdft_charge_constraint) 1016 ! Do nothing 1017 CASE (cdft_magnetization_constraint) 1018 needs_spin_density = .TRUE. 1019 CASE (cdft_alpha_constraint, cdft_beta_constraint) 1020 CALL cp_abort(__LOCATION__, & 1021 "CDFT fragment constraint not yet compatible with "// & 1022 "spin specific constraints.") 1023 CASE DEFAULT 1024 CPABORT("Unknown constraint type.") 1025 END SELECT 1026 END DO 1027 IF (needs_spin_density) THEN 1028 nfrag_spins = 2 1029 DO i = 1, 2 1030 IF (cdft_control%flip_fragment(i)) multiplier(i) = -1.0_dp 1031 END DO 1032 END IF 1033 ! Read fragment reference densities 1034 ALLOCATE (cdft_control%fragments(nfrag_spins, 2)) 1035 ALLOCATE (rho_frag(nfrag_spins)) 1036 CALL get_qs_env(qs_env, pw_env=pw_env) 1037 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 1038 ! Total density (rho_alpha + rho_beta) 1039 CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%fragments(1, 1)%pw, & 1040 use_data=REALDATA3D, in_space=REALSPACE) 1041 CALL cp_cube_to_pw(cdft_control%fragments(1, 1)%pw, & 1042 cdft_control%fragment_a_fname, 1.0_dp) 1043 CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%fragments(1, 2)%pw, & 1044 use_data=REALDATA3D, in_space=REALSPACE) 1045 CALL cp_cube_to_pw(cdft_control%fragments(1, 2)%pw, & 1046 cdft_control%fragment_b_fname, 1.0_dp) 1047 ! Spin difference density (rho_alpha - rho_beta) if needed 1048 IF (needs_spin_density) THEN 1049 CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%fragments(2, 1)%pw, & 1050 use_data=REALDATA3D, in_space=REALSPACE) 1051 CALL cp_cube_to_pw(cdft_control%fragments(2, 1)%pw, & 1052 cdft_control%fragment_a_spin_fname, multiplier(1)) 1053 CALL pw_pool_create_pw(auxbas_pw_pool, cdft_control%fragments(2, 2)%pw, & 1054 use_data=REALDATA3D, in_space=REALSPACE) 1055 CALL cp_cube_to_pw(cdft_control%fragments(2, 2)%pw, & 1056 cdft_control%fragment_b_spin_fname, multiplier(2)) 1057 END IF 1058 ! Sum up fragments 1059 DO i = 1, nfrag_spins 1060 CALL pw_pool_create_pw(auxbas_pw_pool, rho_frag(i)%pw, use_data=REALDATA3D, & 1061 in_space=REALSPACE) 1062 CALL pw_copy(cdft_control%fragments(i, 1)%pw, rho_frag(i)%pw) 1063 CALL pw_axpy(cdft_control%fragments(i, 2)%pw, rho_frag(i)%pw, 1.0_dp) 1064 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%fragments(i, 1)%pw) 1065 CALL pw_pool_give_back_pw(auxbas_pw_pool, cdft_control%fragments(i, 2)%pw) 1066 END DO 1067 DEALLOCATE (cdft_control%fragments) 1068 ! Check that the number of electrons is consistent 1069 CALL get_qs_env(qs_env, subsys=subsys) 1070 CALL qs_subsys_get(subsys, nelectron_total=nelectron_total) 1071 nelectron_frag = pw_integrate_function(rho_frag(1)%pw) 1072 IF (NINT(nelectron_frag) /= nelectron_total) & 1073 CALL cp_abort(__LOCATION__, & 1074 "The number of electrons in the reference and interacting "// & 1075 "configurations does not match. Check your fragment cube files.") 1076 ! Update constraint target value i.e. perform integration w_i*rho_frag_{tot/spin}*dr 1077 cdft_control%target = 0.0_dp 1078 DO igroup = 1, SIZE(group) 1079 IF (group(igroup)%constraint_type == cdft_charge_constraint) THEN 1080 i = 1 1081 ELSE 1082 i = 2 1083 END IF 1084 IF (is_becke .AND. (cdft_control%external_control .AND. becke_control%cavity_confine)) THEN 1085 cdft_control%target(igroup) = cdft_control%target(igroup) + & 1086 accurate_dot_product(group(igroup)%weight%pw%cr3d, rho_frag(i)%pw%cr3d, & 1087 becke_control%cavity_mat, becke_control%eps_cavity)*dvol 1088 ELSE 1089 cdft_control%target(igroup) = cdft_control%target(igroup) + & 1090 accurate_sum(group(igroup)%weight%pw%cr3d*rho_frag(i)%pw%cr3d)*dvol 1091 END IF 1092 END DO 1093 CALL mp_sum(cdft_control%target, para_env%group) 1094 ! Calculate reference atomic charges int( w_i * rho_frag * dr ) 1095 IF (cdft_control%atomic_charges) THEN 1096 ALLOCATE (cdft_control%charges_fragment(cdft_control%natoms, nfrag_spins)) 1097 DO i = 1, nfrag_spins 1098 DO iatom = 1, cdft_control%natoms 1099 cdft_control%charges_fragment(iatom, i) = accurate_sum(cdft_control%charge(iatom)%pw%cr3d*rho_frag(i)%pw%cr3d)*dvol 1100 END DO 1101 END DO 1102 CALL mp_sum(cdft_control%charges_fragment, para_env%group) 1103 END IF 1104 DO i = 1, nfrag_spins 1105 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_frag(i)%pw) 1106 END DO 1107 DEALLOCATE (rho_frag) 1108 cdft_control%fragments_integrated = .TRUE. 1109 1110 CALL timestop(handle) 1111 1112 END SUBROUTINE prepare_fragment_constraint 1113 1114! ************************************************************************************************** 1115!> \brief Driver routine for calculating a Gaussian Hirshfeld constraint 1116!> \param qs_env the qs_env where to build the constraint 1117!> \param calc_pot if the potential needs to be recalculated or just integrated 1118!> \param calculate_forces logical if potential has to be calculated or only_energy 1119!> \par History 1120!> Created 09.2018 [Nico Holmberg] 1121! ************************************************************************************************** 1122 SUBROUTINE hirshfeld_constraint(qs_env, calc_pot, calculate_forces) 1123 TYPE(qs_environment_type), POINTER :: qs_env 1124 LOGICAL :: calc_pot, calculate_forces 1125 1126 CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint', & 1127 routineP = moduleN//':'//routineN 1128 1129 INTEGER :: handle 1130 TYPE(cdft_control_type), POINTER :: cdft_control 1131 TYPE(dft_control_type), POINTER :: dft_control 1132 1133 CALL timeset(routineN, handle) 1134 CALL get_qs_env(qs_env, dft_control=dft_control) 1135 cdft_control => dft_control%qs_control%cdft_control 1136 IF (dft_control%qs_control%cdft .AND. cdft_control%type == outer_scf_hirshfeld_constraint) THEN 1137 IF (calc_pot) THEN 1138 ! Initialize the Hirshfeld constraint environment 1139 CALL hirshfeld_constraint_init(qs_env) 1140 ! Calculate the Hirshfeld weight function and possibly the gradients 1141 CALL hirshfeld_constraint_low(qs_env) 1142 END IF 1143 ! Integrate the Hirshfeld constraint 1144 CALL cdft_constraint_integrate(qs_env) 1145 ! Calculate forces 1146 IF (calculate_forces) CPABORT("Hirshfeld force NYI.") 1147 !CALL hirshfeld_constraint_force(qs_env) 1148 END IF 1149 CALL timestop(handle) 1150 1151 END SUBROUTINE hirshfeld_constraint 1152 1153! ************************************************************************************************** 1154!> \brief Calculates Gaussian Hirshfeld constraints 1155!> \param qs_env the qs_env where to build the constraint 1156!> \author Nico Holmberg (09.2018) 1157! ************************************************************************************************** 1158 SUBROUTINE hirshfeld_constraint_low(qs_env) 1159 TYPE(qs_environment_type), POINTER :: qs_env 1160 1161 CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_low', & 1162 routineP = moduleN//':'//routineN 1163 1164 INTEGER :: atom_a, handle, i, iatom, iex, igroup, & 1165 ikind, ithread, j, natom, npme, & 1166 nthread, numexp 1167 INTEGER(KIND=int_8) :: subpatch_pattern 1168 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 1169 LOGICAL, ALLOCATABLE, DIMENSION(:) :: compute_charge, is_constraint 1170 REAL(kind=dp) :: alpha, coef, dvol, eps_rho_rspace, & 1171 radius, radius_constr 1172 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coefficients 1173 REAL(kind=dp), DIMENSION(3) :: ra 1174 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 1175 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1176 TYPE(cdft_control_type), POINTER :: cdft_control 1177 TYPE(cell_type), POINTER :: cell 1178 TYPE(cube_info_type) :: cube_info 1179 TYPE(dft_control_type), POINTER :: dft_control 1180 TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control 1181 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 1182 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1183 TYPE(pw_env_type), POINTER :: pw_env 1184 TYPE(pw_p_type) :: tmp 1185 TYPE(pw_p_type), POINTER :: fnorm 1186 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1187 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 1188 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_single 1189 TYPE(realspace_grid_type), POINTER :: rs_rho_all, rs_rho_constr 1190 1191 NULLIFY (atom_list, cores, pab, atomic_kind_set, cell, dft_control, & 1192 hirshfeld_env, particle_set, pw_env, fnorm, auxbas_pw_pool, & 1193 auxbas_rs_desc, rs_rho_all, rs_rho_constr, cdft_control, & 1194 hirshfeld_control, rs_single) 1195 1196 CALL timeset(routineN, handle) 1197 CALL get_qs_env(qs_env, & 1198 atomic_kind_set=atomic_kind_set, & 1199 cell=cell, & 1200 particle_set=particle_set, & 1201 natom=natom, & 1202 dft_control=dft_control, & 1203 pw_env=pw_env) 1204 1205 cdft_control => dft_control%qs_control%cdft_control 1206 hirshfeld_control => cdft_control%hirshfeld_control 1207 hirshfeld_env => hirshfeld_control%hirshfeld_env 1208 1209 ! Build weight function 1210 ALLOCATE (is_constraint(natom)) 1211 ALLOCATE (coefficients(natom)) 1212 IF (cdft_control%atomic_charges) THEN 1213 ALLOCATE (compute_charge(natom)) 1214 compute_charge = .FALSE. 1215 END IF 1216 1217 DO igroup = 1, SIZE(cdft_control%group) 1218 dvol = cdft_control%group(igroup)%weight%pw%pw_grid%dvol 1219 cdft_control%group(igroup)%weight%pw%cr3d = 0.0_dp 1220 ! Keep track of those atoms that are active in this constraint group 1221 coefficients(:) = 0.0_dp 1222 is_constraint = .FALSE. 1223 DO i = 1, SIZE(cdft_control%group(igroup)%atoms) 1224 coefficients(cdft_control%group(igroup)%atoms(i)) = cdft_control%group(igroup)%coeff(i) 1225 is_constraint(cdft_control%group(igroup)%atoms(i)) = .TRUE. 1226 END DO 1227 ! Calculate atomic Hirshfeld weight functions: rs_rho_constr / rs_rho_all 1228 ! rs_rho_all: Superposition of isolated Gaussian densities over all atoms in system 1229 ! rs_rho_constr: Sum of isolated Gaussian densities over constraint atoms in this constraint group 1230 IF (.NOT. ASSOCIATED(rs_rho_all)) THEN 1231 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho_all, & 1232 auxbas_pw_pool=auxbas_pw_pool) 1233 cube_info = pw_env%cube_info(1) 1234 CALL rs_grid_retain(rs_rho_all) 1235 CALL rs_grid_zero(rs_rho_all) 1236 END IF 1237 CALL rs_grid_create(rs_rho_constr, auxbas_rs_desc) 1238 CALL rs_grid_zero(rs_rho_constr) 1239 ! Compute Gaussian density over single atoms (rs_single) when atomic charges are requested 1240 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN 1241 ALLOCATE (rs_single(cdft_control%natoms)) 1242 DO i = 1, cdft_control%natoms 1243 NULLIFY (rs_single(i)%rs_grid) 1244 CALL rs_grid_create(rs_single(i)%rs_grid, auxbas_rs_desc) 1245 CALL rs_grid_zero(rs_single(i)%rs_grid) 1246 compute_charge(cdft_control%atoms(i)) = .TRUE. 1247 END DO 1248 END IF 1249 1250 ! Collocate Gaussians 1251 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 1252 ALLOCATE (pab(1, 1)) 1253 nthread = 1 1254 ithread = 0 1255 1256 DO ikind = 1, SIZE(atomic_kind_set) 1257 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp 1258 IF (numexp <= 0) CYCLE 1259 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list) 1260 ALLOCATE (cores(natom)) 1261 1262 DO iex = 1, numexp 1263 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) 1264 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) 1265 npme = 0 1266 cores = 0 1267 DO iatom = 1, natom 1268 atom_a = atom_list(iatom) 1269 ra(:) = pbc(particle_set(atom_a)%r, cell) 1270 IF (rs_rho_all%desc%parallel .AND. .NOT. rs_rho_all%desc%distributed) THEN 1271 ! replicated realspace grid, split the atoms up between procs 1272 IF (MODULO(iatom, rs_rho_all%desc%group_size) == rs_rho_all%desc%my_pos) THEN 1273 npme = npme + 1 1274 cores(npme) = iatom 1275 ENDIF 1276 ELSE 1277 npme = npme + 1 1278 cores(npme) = iatom 1279 ENDIF 1280 END DO 1281 DO j = 1, npme 1282 iatom = cores(j) 1283 atom_a = atom_list(iatom) 1284 pab(1, 1) = coef*hirshfeld_env%charges(atom_a) 1285 ra(:) = pbc(particle_set(atom_a)%r, cell) 1286 subpatch_pattern = 0 1287 radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & 1288 ra=ra, rb=ra, rp=ra, & 1289 zetp=alpha, eps=eps_rho_rspace, & 1290 pab=pab, o1=0, o2=0, & ! without map_consistent 1291 prefactor=1.0_dp, cutoff=0.0_dp) 1292 1293 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & 1294 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & 1295 rs_rho_all, cell, cube_info, radius=radius, & 1296 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & 1297 subpatch_pattern=subpatch_pattern) 1298 IF (is_constraint(atom_a)) THEN 1299 radius_constr = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, & 1300 ra=ra, rb=ra, rp=ra, & 1301 zetp=alpha, eps=eps_rho_rspace, & 1302 pab=pab, o1=0, o2=0, & ! without map_consistent 1303 prefactor=coefficients(atom_a), cutoff=0.0_dp) 1304 1305 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & 1306 (/0.0_dp, 0.0_dp, 0.0_dp/), coefficients(atom_a), & 1307 pab, 0, 0, rs_rho_constr, cell, cube_info, & 1308 radius=radius_constr, & 1309 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & 1310 subpatch_pattern=subpatch_pattern) 1311 END IF 1312 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN 1313 IF (compute_charge(atom_a)) THEN 1314 DO iatom = 1, cdft_control%natoms 1315 IF (atom_a == cdft_control%atoms(iatom)) EXIT 1316 END DO 1317 CPASSERT(iatom <= cdft_control%natoms) 1318 CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, & 1319 (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, & 1320 rs_single(iatom)%rs_grid, cell, cube_info, radius=radius, & 1321 ga_gb_function=GRID_FUNC_AB, use_subpatch=.TRUE., & 1322 subpatch_pattern=subpatch_pattern) 1323 END IF 1324 END IF 1325 END DO 1326 END DO 1327 DEALLOCATE (cores) 1328 END DO 1329 DEALLOCATE (pab) 1330 1331 ! Transfer rs_rho_all to the correct grid and save it 1332 CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm) 1333 IF (igroup == 1) THEN 1334 IF (ASSOCIATED(fnorm)) THEN 1335 CALL pw_pool_give_back_pw(auxbas_pw_pool, fnorm%pw) 1336 END IF 1337 ALLOCATE (fnorm) 1338 CALL pw_pool_create_pw(auxbas_pw_pool, fnorm%pw, use_data=REALDATA3D, & 1339 in_space=REALSPACE) 1340 CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm) 1341 CALL rs_pw_transfer(rs_rho_all, fnorm%pw, rs2pw) 1342 CALL rs_grid_release(rs_rho_all) 1343 END IF 1344 ! Compute CDFT weight function 1345 CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, & 1346 in_space=REALSPACE) 1347 CALL rs_pw_transfer(rs_rho_constr, tmp%pw, rs2pw) 1348 CALL rs_grid_release(rs_rho_constr) 1349 CALL hfun_scale(cdft_control%group(igroup)%weight%pw%cr3d, tmp%pw%cr3d, & 1350 fnorm%pw%cr3d, divide=.TRUE.) 1351 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw) 1352 ! Compute atomic weight functions if charges are needed 1353 IF (cdft_control%atomic_charges .AND. igroup == 1) THEN 1354 CALL pw_pool_create_pw(auxbas_pw_pool, tmp%pw, use_data=REALDATA3D, & 1355 in_space=REALSPACE) 1356 DO i = 1, cdft_control%natoms 1357 CALL rs_pw_transfer(rs_single(i)%rs_grid, tmp%pw, rs2pw) 1358 CALL rs_grid_release(rs_single(i)%rs_grid) 1359 CALL hfun_scale(cdft_control%charge(i)%pw%cr3d, tmp%pw%cr3d, & 1360 fnorm%pw%cr3d, divide=.TRUE.) 1361 END DO 1362 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp%pw) 1363 DEALLOCATE (rs_single) 1364 DEALLOCATE (compute_charge) 1365 END IF 1366 END DO 1367 1368 DEALLOCATE (is_constraint) 1369 DEALLOCATE (coefficients) 1370 CALL timestop(handle) 1371 1372 END SUBROUTINE hirshfeld_constraint_low 1373 1374END MODULE qs_cdft_methods 1375