1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utility subroutines for CDFT calculations 8!> \par History 9!> separated from et_coupling [03.2017] 10!> \author Nico Holmberg [03.2017] 11! ************************************************************************************************** 12MODULE qs_cdft_utils 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind 15 USE bibliography, ONLY: Becke1988b,& 16 Holmberg2017,& 17 Holmberg2018,& 18 cite_reference 19 USE cell_types, ONLY: cell_type,& 20 pbc 21 USE cp_control_types, ONLY: dft_control_type,& 22 qs_control_type 23 USE cp_log_handling, ONLY: cp_get_default_logger,& 24 cp_logger_type,& 25 cp_to_string 26 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 27 cp_print_key_unit_nr 28 USE cp_para_types, ONLY: cp_para_env_type 29 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 30 USE cp_units, ONLY: cp_unit_from_cp2k 31 USE hirshfeld_methods, ONLY: create_shape_function 32 USE hirshfeld_types, ONLY: create_hirshfeld_type,& 33 hirshfeld_type,& 34 set_hirshfeld_info 35 USE input_constants, ONLY: & 36 becke_cutoff_element, becke_cutoff_global, cdft_charge_constraint, & 37 outer_scf_becke_constraint, outer_scf_cdft_constraint, outer_scf_hirshfeld_constraint, & 38 outer_scf_none, radius_user, shape_function_gaussian 39 USE input_section_types, ONLY: section_get_ivals,& 40 section_vals_get,& 41 section_vals_get_subs_vals,& 42 section_vals_type,& 43 section_vals_val_get 44 USE kinds, ONLY: default_path_length,& 45 dp,& 46 int_8 47 USE memory_utilities, ONLY: reallocate 48 USE outer_scf_control_types, ONLY: outer_scf_read_parameters 49 USE particle_list_types, ONLY: particle_list_type 50 USE particle_types, ONLY: particle_type 51 USE pw_env_types, ONLY: pw_env_get,& 52 pw_env_type 53 USE pw_pool_types, ONLY: pw_pool_create_pw,& 54 pw_pool_type 55 USE pw_types, ONLY: REALDATA3D,& 56 REALSPACE 57 USE qs_cdft_types, ONLY: becke_constraint_type,& 58 cdft_control_type,& 59 cdft_group_type,& 60 hirshfeld_constraint_type 61 USE qs_collocate_density, ONLY: collocate_pgf_product_rspace 62 USE qs_environment_types, ONLY: get_qs_env,& 63 qs_environment_type 64 USE qs_kind_types, ONLY: get_qs_kind,& 65 qs_kind_type 66 USE qs_modify_pab_block, ONLY: FUNC_AB 67 USE qs_scf_output, ONLY: qs_scf_cdft_constraint_info 68 USE qs_subsys_types, ONLY: qs_subsys_get,& 69 qs_subsys_type 70 USE realspace_grid_types, ONLY: realspace_grid_type,& 71 rs2pw,& 72 rs_grid_release,& 73 rs_grid_retain,& 74 rs_grid_zero,& 75 rs_pw_transfer 76#include "./base/base_uses.f90" 77 78 IMPLICIT NONE 79 80 PRIVATE 81 82 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_utils' 83 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 84 85! *** Public subroutines *** 86 PUBLIC :: becke_constraint_init, read_becke_section, read_cdft_control_section 87 PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print 88 89CONTAINS 90 91! ************************************************************************************************** 92!> \brief Initializes the Becke constraint environment 93!> \param qs_env the qs_env where to build the constraint 94!> \par History 95!> Created 01.2007 [fschiff] 96!> Extended functionality 12/15-12/16 [Nico Holmberg] 97! ************************************************************************************************** 98 SUBROUTINE becke_constraint_init(qs_env) 99 TYPE(qs_environment_type), POINTER :: qs_env 100 101 CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_init', & 102 routineP = moduleN//':'//routineN 103 104 CHARACTER(len=2) :: element_symbol 105 INTEGER :: atom_a, bounds(2), handle, i, iatom, iex, igroup, ikind, ip, ithread, iw, j, & 106 jatom, katom, natom, nkind, npme, nthread, numexp, unit_nr 107 INTEGER, DIMENSION(2, 3) :: bo 108 INTEGER, DIMENSION(:), POINTER :: atom_list, cores, stride 109 LOGICAL :: build, in_memory, mpi_io 110 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint 111 REAL(KIND=dp) :: alpha, chi, coef, eps_cavity, ircov, & 112 jrcov, uij 113 REAL(KIND=dp), DIMENSION(3) :: cell_v, dist_vec, r, r1, ra 114 REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list 115 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 116 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 117 TYPE(becke_constraint_type), POINTER :: becke_control 118 TYPE(cdft_control_type), POINTER :: cdft_control 119 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 120 TYPE(cell_type), POINTER :: cell 121 TYPE(cp_logger_type), POINTER :: logger 122 TYPE(cp_para_env_type), POINTER :: para_env 123 TYPE(dft_control_type), POINTER :: dft_control 124 TYPE(hirshfeld_type), POINTER :: cavity_env 125 TYPE(particle_list_type), POINTER :: particles 126 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 127 TYPE(pw_env_type), POINTER :: pw_env 128 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 129 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 130 TYPE(qs_subsys_type), POINTER :: subsys 131 TYPE(realspace_grid_type), POINTER :: rs_cavity 132 TYPE(section_vals_type), POINTER :: cdft_constraint_section 133 134 NULLIFY (cores, stride, atom_list, cell, para_env, dft_control, & 135 particle_set, logger, cdft_constraint_section, qs_kind_set, & 136 particles, subsys, pab, pw_env, rs_cavity, cavity_env, & 137 auxbas_pw_pool, atomic_kind_set, group, radii_list, cdft_control) 138 logger => cp_get_default_logger() 139 CALL timeset(routineN, handle) 140 CALL get_qs_env(qs_env, & 141 cell=cell, & 142 particle_set=particle_set, & 143 natom=natom, & 144 dft_control=dft_control, & 145 para_env=para_env) 146 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") 147 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog") 148 cdft_control => dft_control%qs_control%cdft_control 149 becke_control => cdft_control%becke_control 150 group => cdft_control%group 151 in_memory = .FALSE. 152 IF (cdft_control%save_pot) THEN 153 in_memory = becke_control%in_memory 154 END IF 155 IF (becke_control%cavity_confine) THEN 156 ALLOCATE (is_constraint(natom)) 157 is_constraint = .FALSE. 158 DO i = 1, cdft_control%natoms 159 ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges 160 ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct 161 is_constraint(cdft_control%atoms(i)) = .TRUE. 162 END DO 163 END IF 164 eps_cavity = becke_control%eps_cavity 165 ! Setup atomic radii for adjusting cell boundaries 166 IF (becke_control%adjust) THEN 167 IF (.NOT. ASSOCIATED(becke_control%radii)) THEN 168 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set) 169 IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%radii_tmp)) & 170 CALL cp_abort(__LOCATION__, & 171 "Length of keyword BECKE_CONSTRAINT\ATOMIC_RADII does not "// & 172 "match number of atomic kinds in the input coordinate file.") 173 ALLOCATE (becke_control%radii(SIZE(atomic_kind_set))) 174 becke_control%radii(:) = becke_control%radii_tmp(:) 175 DEALLOCATE (becke_control%radii_tmp) 176 END IF 177 END IF 178 ! Setup cutoff scheme 179 IF (.NOT. ASSOCIATED(becke_control%cutoffs)) THEN 180 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set) 181 ALLOCATE (becke_control%cutoffs(natom)) 182 SELECT CASE (becke_control%cutoff_type) 183 CASE (becke_cutoff_global) 184 becke_control%cutoffs(:) = becke_control%rglobal 185 CASE (becke_cutoff_element) 186 IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%cutoffs_tmp)) & 187 CALL cp_abort(__LOCATION__, & 188 "Length of keyword BECKE_CONSTRAINT\ELEMENT_CUTOFFS does not "// & 189 "match number of atomic kinds in the input coordinate file.") 190 DO ikind = 1, SIZE(atomic_kind_set) 191 CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list) 192 DO iatom = 1, katom 193 atom_a = atom_list(iatom) 194 becke_control%cutoffs(atom_a) = becke_control%cutoffs_tmp(ikind) 195 END DO 196 END DO 197 DEALLOCATE (becke_control%cutoffs_tmp) 198 END SELECT 199 END IF 200 ! Zero weight functions 201 DO igroup = 1, SIZE(group) 202 group(igroup)%weight%pw%cr3d = 0.0_dp 203 END DO 204 IF (cdft_control%atomic_charges) THEN 205 DO iatom = 1, cdft_control%natoms 206 cdft_control%charge(iatom)%pw%cr3d = 0.0_dp 207 END DO 208 END IF 209 ! Allocate storage for cell adjustment coefficients and needed distance vectors 210 build = .FALSE. 211 IF (becke_control%adjust .AND. .NOT. ASSOCIATED(becke_control%aij)) THEN 212 ALLOCATE (becke_control%aij(natom, natom)) 213 build = .TRUE. 214 END IF 215 IF (becke_control%vector_buffer%store_vectors) THEN 216 ALLOCATE (becke_control%vector_buffer%distances(natom)) 217 ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom)) 218 IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom)) 219 ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom)) 220 END IF 221 ALLOCATE (becke_control%vector_buffer%R12(natom, natom)) 222 ! Calculate pairwise distances between each atom pair 223 DO i = 1, 3 224 cell_v(i) = cell%hmat(i, i) 225 END DO 226 DO iatom = 1, natom - 1 227 DO jatom = iatom + 1, natom 228 r = particle_set(iatom)%r 229 r1 = particle_set(jatom)%r 230 DO i = 1, 3 231 r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp 232 r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp 233 END DO 234 dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v 235 ! Store pbc corrected position and pairwise distance vectors for later reuse 236 IF (becke_control%vector_buffer%store_vectors) THEN 237 becke_control%vector_buffer%position_vecs(:, iatom) = r(:) 238 IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:) 239 IF (in_memory) THEN 240 becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:) 241 becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:) 242 END IF 243 END IF 244 becke_control%vector_buffer%R12(iatom, jatom) = SQRT(DOT_PRODUCT(dist_vec, dist_vec)) 245 becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom) 246 ! Set up heteronuclear cell partitioning using user defined radii 247 IF (build) THEN 248 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=ikind) 249 ircov = becke_control%radii(ikind) 250 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, kind_number=ikind) 251 jrcov = becke_control%radii(ikind) 252 IF (ircov .NE. jrcov) THEN 253 chi = ircov/jrcov 254 uij = (chi - 1.0_dp)/(chi + 1.0_dp) 255 becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp) 256 IF (becke_control%aij(iatom, jatom) .GT. 0.5_dp) THEN 257 becke_control%aij(iatom, jatom) = 0.5_dp 258 ELSE IF (becke_control%aij(iatom, jatom) .LT. -0.5_dp) THEN 259 becke_control%aij(iatom, jatom) = -0.5_dp 260 END IF 261 ELSE 262 becke_control%aij(iatom, jatom) = 0.0_dp 263 END IF 264 ! Note change of sign 265 becke_control%aij(jatom, iatom) = -becke_control%aij(iatom, jatom) 266 END IF 267 END DO 268 END DO 269 ! Dump some additional information about the calculation 270 IF (cdft_control%first_iteration) THEN 271 IF (iw > 0) THEN 272 WRITE (iw, '(/,T3,A)') & 273 '----------------------- Becke atomic parameters ------------------------' 274 IF (becke_control%adjust) THEN 275 WRITE (iw, '(T3,A)') & 276 'Atom Element Cutoff (angstrom) CDFT Radius (angstrom)' 277 DO iatom = 1, natom 278 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol, & 279 kind_number=ikind) 280 ircov = cp_unit_from_cp2k(becke_control%radii(ikind), "angstrom") 281 WRITE (iw, "(i6,T15,A2,T37,F8.3,T67,F8.3)") & 282 iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(becke_control%cutoffs(iatom), "angstrom"), & 283 ircov 284 END DO 285 ELSE 286 WRITE (iw, '(T3,A)') & 287 'Atom Element Cutoff (angstrom)' 288 DO iatom = 1, natom 289 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol) 290 WRITE (iw, "(i7,T15,A2,T37,F8.3)") & 291 iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(becke_control%cutoffs(iatom), "angstrom") 292 END DO 293 END IF 294 WRITE (iw, '(T3,A)') & 295 '------------------------------------------------------------------------' 296 WRITE (iw, '(/,T3,A,T60)') & 297 '----------------------- Becke group definitions ------------------------' 298 DO igroup = 1, SIZE(group) 299 IF (igroup > 1) WRITE (iw, '(T3,A)') ' ' 300 WRITE (iw, '(T5,A,I5,A,I5)') & 301 'Atomic group', igroup, ' of ', SIZE(group) 302 WRITE (iw, '(T5,A)') 'Atom Element Coefficient' 303 DO ip = 1, SIZE(group(igroup)%atoms) 304 iatom = group(igroup)%atoms(ip) 305 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol) 306 WRITE (iw, '(i8,T16,A2,T23,F8.3)') iatom, ADJUSTR(element_symbol), group(igroup)%coeff(ip) 307 END DO 308 END DO 309 WRITE (iw, '(T3,A)') & 310 '------------------------------------------------------------------------' 311 END IF 312 cdft_control%first_iteration = .FALSE. 313 END IF 314 ! Setup cavity confinement using spherical Gaussians 315 IF (becke_control%cavity_confine) THEN 316 cavity_env => becke_control%cavity_env 317 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, pw_env=pw_env, qs_kind_set=qs_kind_set) 318 CPASSERT(ASSOCIATED(qs_kind_set)) 319 nkind = SIZE(qs_kind_set) 320 ! Setup the Gaussian shape function 321 IF (.NOT. ASSOCIATED(cavity_env%kind_shape_fn)) THEN 322 IF (ASSOCIATED(becke_control%radii)) THEN 323 ALLOCATE (radii_list(SIZE(becke_control%radii))) 324 DO ikind = 1, SIZE(becke_control%radii) 325 IF (cavity_env%use_bohr) THEN 326 radii_list(ikind) = becke_control%radii(ikind) 327 ELSE 328 radii_list(ikind) = cp_unit_from_cp2k(becke_control%radii(ikind), "angstrom") 329 END IF 330 END DO 331 END IF 332 CALL create_shape_function(cavity_env, qs_kind_set, atomic_kind_set, & 333 radius=becke_control%rcavity, & 334 radii_list=radii_list) 335 IF (ASSOCIATED(radii_list)) & 336 DEALLOCATE (radii_list) 337 END IF 338 ! Form cavity by summing isolated Gaussian densities over constraint atoms 339 NULLIFY (rs_cavity) 340 CALL pw_env_get(pw_env, auxbas_rs_grid=rs_cavity, auxbas_pw_pool=auxbas_pw_pool) 341 CALL rs_grid_retain(rs_cavity) 342 CALL rs_grid_zero(rs_cavity) 343 ALLOCATE (pab(1, 1)) 344 nthread = 1 345 ithread = 0 346 DO ikind = 1, SIZE(atomic_kind_set) 347 numexp = cavity_env%kind_shape_fn(ikind)%numexp 348 IF (numexp <= 0) CYCLE 349 CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list) 350 ALLOCATE (cores(katom)) 351 DO iex = 1, numexp 352 alpha = cavity_env%kind_shape_fn(ikind)%zet(iex) 353 coef = cavity_env%kind_shape_fn(ikind)%coef(iex) 354 npme = 0 355 cores = 0 356 DO iatom = 1, katom 357 IF (rs_cavity%desc%parallel .AND. .NOT. rs_cavity%desc%distributed) THEN 358 ! replicated realspace grid, split the atoms up between procs 359 IF (MODULO(iatom, rs_cavity%desc%group_size) == rs_cavity%desc%my_pos) THEN 360 npme = npme + 1 361 cores(npme) = iatom 362 ENDIF 363 ELSE 364 npme = npme + 1 365 cores(npme) = iatom 366 ENDIF 367 END DO 368 DO j = 1, npme 369 iatom = cores(j) 370 atom_a = atom_list(iatom) 371 pab(1, 1) = coef 372 IF (becke_control%vector_buffer%store_vectors) THEN 373 ra(:) = becke_control%vector_buffer%position_vecs(:, atom_a) + cell_v(:)/2._dp 374 ELSE 375 ra(:) = pbc(particle_set(atom_a)%r, cell) 376 END IF 377 IF (is_constraint(atom_a)) & 378 CALL collocate_pgf_product_rspace(0, alpha, 0, 0, 0.0_dp, 0, ra, & 379 (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, 1.0_dp, & 380 pab, 0, 0, rs_cavity, cell, pw_env%cube_info(1), & 381 dft_control%qs_control%eps_rho_rspace, & 382 ga_gb_function=FUNC_AB, ithread=ithread, & 383 use_subpatch=.TRUE., subpatch_pattern=0_int_8, & 384 lmax_global=0) 385 END DO 386 END DO 387 DEALLOCATE (cores) 388 END DO 389 DEALLOCATE (pab) 390 CALL pw_pool_create_pw(auxbas_pw_pool, becke_control%cavity%pw, & 391 use_data=REALDATA3D, in_space=REALSPACE) 392 CALL rs_pw_transfer(rs_cavity, becke_control%cavity%pw, rs2pw) 393 CALL rs_grid_release(rs_cavity) 394 ! Grid points where the Gaussian density falls below eps_cavity are ignored 395 ! We can calculate the smallest/largest values along z-direction outside 396 ! which the cavity is zero at every point (x, y) 397 ! If gradients are needed storage needs to be allocated only for grid points within 398 ! these bounds 399 IF (in_memory .OR. cdft_control%save_pot) THEN 400 CALL hfun_zero(becke_control%cavity%pw%cr3d, eps_cavity, just_bounds=.TRUE., bounds=bounds) 401 ! Save bounds (first nonzero grid point indices) 402 bo = group(1)%weight%pw%pw_grid%bounds_local 403 IF (bounds(2) .LT. bo(2, 3)) THEN 404 bounds(2) = bounds(2) - 1 405 ELSE 406 bounds(2) = bo(2, 3) 407 END IF 408 IF (bounds(1) .GT. bo(1, 3)) THEN 409 ! In the special case bounds(1) == bounds(2) == bo(2, 3), after this check 410 ! bounds(1) > bounds(2) and the subsequent gradient allocation (:, :, :, bounds(1):bounds(2)) 411 ! will correctly allocate a 0-sized array 412 bounds(1) = bounds(1) + 1 413 ELSE 414 bounds(1) = bo(1, 3) 415 END IF 416 becke_control%confine_bounds = bounds 417 END IF 418 ! Optional printing of cavity (meant for testing, so options currently hardcoded...) 419 IF (becke_control%print_cavity) THEN 420 CALL hfun_zero(becke_control%cavity%pw%cr3d, eps_cavity, just_bounds=.FALSE.) 421 ALLOCATE (stride(3)) 422 stride = (/2, 2, 2/) 423 mpi_io = .TRUE. 424 ! Note PROGRAM_RUN_INFO section neeeds to be active! 425 unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", & 426 middle_name="BECKE_CAVITY", & 427 extension=".cube", file_position="REWIND", & 428 log_filename=.FALSE., mpi_io=mpi_io) 429 IF (para_env%ionode .AND. unit_nr .LT. 1) & 430 CALL cp_abort(__LOCATION__, & 431 "Please turn on PROGRAM_RUN_INFO to print cavity") 432 CALL get_qs_env(qs_env, subsys=subsys) 433 CALL qs_subsys_get(subsys, particles=particles) 434 CALL cp_pw_to_cube(becke_control%cavity%pw, unit_nr, "CAVITY", particles=particles, stride=stride, mpi_io=mpi_io) 435 CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io) 436 DEALLOCATE (stride) 437 END IF 438 END IF 439 IF (ALLOCATED(is_constraint)) & 440 DEALLOCATE (is_constraint) 441 CALL timestop(handle) 442 443 END SUBROUTINE becke_constraint_init 444 445! ************************************************************************************************** 446!> \brief reads the input parameters specific to Becke-based CDFT constraints 447!> \param cdft_control the cdft_control which holds the Becke control type 448!> \param becke_section the input section containing Becke constraint information 449!> \par History 450!> Created 01.2007 [fschiff] 451!> Merged Becke into CDFT 09.2018 [Nico Holmberg] 452!> \author Nico Holmberg [09.2018] 453! ************************************************************************************************** 454 SUBROUTINE read_becke_section(cdft_control, becke_section) 455 456 TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control 457 TYPE(section_vals_type), POINTER :: becke_section 458 459 CHARACTER(len=*), PARAMETER :: routineN = 'read_becke_section', & 460 routineP = moduleN//':'//routineN 461 462 INTEGER :: j 463 LOGICAL :: exists 464 REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist 465 TYPE(becke_constraint_type), POINTER :: becke_control 466 467 NULLIFY (rtmplist) 468 becke_control => cdft_control%becke_control 469 CPASSERT(ASSOCIATED(becke_control)) 470 471 ! Atomic size corrections 472 CALL section_vals_val_get(becke_section, "ADJUST_SIZE", l_val=becke_control%adjust) 473 IF (becke_control%adjust) THEN 474 CALL section_vals_val_get(becke_section, "ATOMIC_RADII", explicit=exists) 475 IF (.NOT. exists) CPABORT("Keyword ATOMIC_RADII is missing.") 476 CALL section_vals_val_get(becke_section, "ATOMIC_RADII", r_vals=rtmplist) 477 CPASSERT(SIZE(rtmplist) > 0) 478 ALLOCATE (becke_control%radii_tmp(SIZE(rtmplist))) 479 DO j = 1, SIZE(rtmplist) 480 becke_control%radii_tmp(j) = rtmplist(j) 481 END DO 482 END IF 483 484 ! Cutoff scheme 485 CALL section_vals_val_get(becke_section, "CUTOFF_TYPE", i_val=becke_control%cutoff_type) 486 SELECT CASE (becke_control%cutoff_type) 487 CASE (becke_cutoff_global) 488 CALL section_vals_val_get(becke_section, "GLOBAL_CUTOFF", r_val=becke_control%rglobal) 489 CASE (becke_cutoff_element) 490 CALL section_vals_val_get(becke_section, "ELEMENT_CUTOFF", r_vals=rtmplist) 491 CPASSERT(SIZE(rtmplist) > 0) 492 ALLOCATE (becke_control%cutoffs_tmp(SIZE(rtmplist))) 493 DO j = 1, SIZE(rtmplist) 494 becke_control%cutoffs_tmp(j) = rtmplist(j) 495 END DO 496 END SELECT 497 498 ! Gaussian cavity confinement 499 CALL section_vals_val_get(becke_section, "CAVITY_CONFINE", l_val=becke_control%cavity_confine) 500 CALL section_vals_val_get(becke_section, "SHOULD_SKIP", l_val=becke_control%should_skip) 501 CALL section_vals_val_get(becke_section, "IN_MEMORY", l_val=becke_control%in_memory) 502 IF (cdft_control%becke_control%cavity_confine) THEN 503 CALL section_vals_val_get(becke_section, "CAVITY_SHAPE", i_val=becke_control%cavity_shape) 504 IF (becke_control%cavity_shape == radius_user .AND. .NOT. becke_control%adjust) & 505 CALL cp_abort(__LOCATION__, & 506 "Activate keyword ADJUST_SIZE to use cavity shape USER.") 507 CALL section_vals_val_get(becke_section, "CAVITY_RADIUS", r_val=becke_control%rcavity) 508 CALL section_vals_val_get(becke_section, "EPS_CAVITY", r_val=becke_control%eps_cavity) 509 CALL section_vals_val_get(becke_section, "CAVITY_PRINT", l_val=becke_control%print_cavity) 510 CALL section_vals_val_get(becke_section, "CAVITY_USE_BOHR", l_val=becke_control%use_bohr) 511 IF (.NOT. cdft_control%becke_control%use_bohr) THEN 512 becke_control%rcavity = cp_unit_from_cp2k(becke_control%rcavity, "angstrom") 513 END IF 514 CALL create_hirshfeld_type(becke_control%cavity_env) 515 CALL set_hirshfeld_info(becke_control%cavity_env, & 516 shape_function_type=shape_function_gaussian, iterative=.FALSE., & 517 radius_type=becke_control%cavity_shape, & 518 use_bohr=becke_control%use_bohr) 519 END IF 520 521 CALL cite_reference(Becke1988b) 522 523 END SUBROUTINE read_becke_section 524 525! ************************************************************************************************** 526!> \brief reads the input parameters needed to define CDFT constraints 527!> \param cdft_control the object which holds the CDFT control type 528!> \param cdft_control_section the input section containing CDFT constraint information 529!> \author Nico Holmberg [09.2018] 530! ************************************************************************************************** 531 SUBROUTINE read_constraint_definitions(cdft_control, cdft_control_section) 532 533 TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control 534 TYPE(section_vals_type), INTENT(INOUT), POINTER :: cdft_control_section 535 536 CHARACTER(len=*), PARAMETER :: routineN = 'read_constraint_definitions', & 537 routineP = moduleN//':'//routineN 538 539 INTEGER :: i, j, jj, k, n_rep, natoms, nvar, & 540 tot_natoms 541 INTEGER, DIMENSION(:), POINTER :: atomlist, dummylist, tmplist 542 LOGICAL :: exists, is_duplicate 543 REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist 544 TYPE(section_vals_type), POINTER :: group_section 545 546 NULLIFY (tmplist, rtmplist, atomlist, dummylist, group_section) 547 548 group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP") 549 CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists) 550 IF (.NOT. exists) CPABORT("Section ATOM_GROUP is missing.") 551 ALLOCATE (cdft_control%group(nvar)) 552 tot_natoms = 0 553 ! Parse all ATOM_GROUP sections 554 DO k = 1, nvar 555 ! First determine how much storage is needed 556 natoms = 0 557 CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, n_rep_val=n_rep) 558 DO j = 1, n_rep 559 CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist) 560 IF (SIZE(tmplist) < 1) & 561 CPABORT("Each ATOM_GROUP must contain at least 1 atom.") 562 natoms = natoms + SIZE(tmplist) 563 END DO 564 ALLOCATE (cdft_control%group(k)%atoms(natoms)) 565 ALLOCATE (cdft_control%group(k)%coeff(natoms)) 566 NULLIFY (cdft_control%group(k)%weight%pw) 567 NULLIFY (cdft_control%group(k)%gradients) 568 NULLIFY (cdft_control%group(k)%integrated) 569 tot_natoms = tot_natoms + natoms 570 ! Now parse 571 jj = 0 572 DO j = 1, n_rep 573 CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist) 574 DO i = 1, SIZE(tmplist) 575 jj = jj + 1 576 cdft_control%group(k)%atoms(jj) = tmplist(i) 577 END DO 578 END DO 579 CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, n_rep_val=n_rep) 580 jj = 0 581 DO j = 1, n_rep 582 CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, i_rep_val=j, r_vals=rtmplist) 583 DO i = 1, SIZE(rtmplist) 584 jj = jj + 1 585 IF (jj > natoms) CPABORT("Length of keywords ATOMS and COEFF must match.") 586 IF (ABS(rtmplist(i)) /= 1.0_dp) CPABORT("Keyword COEFF accepts only values +/-1.0") 587 cdft_control%group(k)%coeff(jj) = rtmplist(i) 588 END DO 589 END DO 590 IF (jj < natoms) CPABORT("Length of keywords ATOMS and COEFF must match.") 591 CALL section_vals_val_get(group_section, "CONSTRAINT_TYPE", i_rep_section=k, & 592 i_val=cdft_control%group(k)%constraint_type) 593 CALL section_vals_val_get(group_section, "FRAGMENT_CONSTRAINT", i_rep_section=k, & 594 l_val=cdft_control%group(k)%is_fragment_constraint) 595 IF (cdft_control%group(k)%is_fragment_constraint) cdft_control%fragment_density = .TRUE. 596 END DO 597 ! Create a list containing all constraint atoms 598 ALLOCATE (atomlist(tot_natoms)) 599 atomlist = -1 600 jj = 0 601 DO k = 1, nvar 602 DO j = 1, SIZE(cdft_control%group(k)%atoms) 603 is_duplicate = .FALSE. 604 DO i = 1, jj + 1 605 IF (cdft_control%group(k)%atoms(j) == atomlist(i)) THEN 606 is_duplicate = .TRUE. 607 EXIT 608 END IF 609 END DO 610 IF (.NOT. is_duplicate) THEN 611 jj = jj + 1 612 atomlist(jj) = cdft_control%group(k)%atoms(j) 613 END IF 614 END DO 615 END DO 616 CALL reallocate(atomlist, 1, jj) 617 CALL section_vals_val_get(cdft_control_section, "ATOMIC_CHARGES", & 618 l_val=cdft_control%atomic_charges) 619 ! Parse any dummy atoms (no constraint, just charges) 620 IF (cdft_control%atomic_charges) THEN 621 group_section => section_vals_get_subs_vals(cdft_control_section, "DUMMY_ATOMS") 622 CALL section_vals_get(group_section, explicit=exists) 623 IF (exists) THEN 624 ! First determine how many atoms there are 625 natoms = 0 626 CALL section_vals_val_get(group_section, "ATOMS", n_rep_val=n_rep) 627 DO j = 1, n_rep 628 CALL section_vals_val_get(group_section, "ATOMS", i_rep_val=j, i_vals=tmplist) 629 IF (SIZE(tmplist) < 1) & 630 CPABORT("DUMMY_ATOMS must contain at least 1 atom.") 631 natoms = natoms + SIZE(tmplist) 632 END DO 633 ALLOCATE (dummylist(natoms)) 634 ! Now parse 635 jj = 0 636 DO j = 1, n_rep 637 CALL section_vals_val_get(group_section, "ATOMS", i_rep_val=j, i_vals=tmplist) 638 DO i = 1, SIZE(tmplist) 639 jj = jj + 1 640 dummylist(jj) = tmplist(i) 641 END DO 642 END DO 643 ! Check for duplicates 644 DO j = 1, natoms 645 DO i = j + 1, natoms 646 IF (dummylist(i) == dummylist(j)) & 647 CPABORT("Duplicate atoms defined in section DUMMY_ATOMS.") 648 END DO 649 END DO 650 ! Check that a dummy atom is not included in any ATOM_GROUP 651 DO j = 1, SIZE(atomlist) 652 DO i = 1, SIZE(dummylist) 653 IF (dummylist(i) == atomlist(j)) & 654 CALL cp_abort(__LOCATION__, & 655 "Duplicate atoms defined in sections ATOM_GROUP and DUMMY_ATOMS.") 656 END DO 657 END DO 658 END IF 659 END IF 660 ! Join dummy atoms and constraint atoms into one list 661 IF (ASSOCIATED(dummylist)) THEN 662 cdft_control%natoms = SIZE(atomlist) + SIZE(dummylist) 663 ELSE 664 cdft_control%natoms = SIZE(atomlist) 665 END IF 666 ALLOCATE (cdft_control%atoms(cdft_control%natoms)) 667 ALLOCATE (cdft_control%is_constraint(cdft_control%natoms)) 668 IF (cdft_control%atomic_charges) ALLOCATE (cdft_control%charge(cdft_control%natoms)) 669 cdft_control%atoms(1:SIZE(atomlist)) = atomlist 670 IF (ASSOCIATED(dummylist)) THEN 671 cdft_control%atoms(1 + SIZE(atomlist):) = dummylist 672 DEALLOCATE (dummylist) 673 END IF 674 cdft_control%is_constraint = .FALSE. 675 cdft_control%is_constraint(1:SIZE(atomlist)) = .TRUE. 676 DEALLOCATE (atomlist) 677 ! Get constraint potential definitions from input 678 ALLOCATE (cdft_control%strength(nvar)) 679 ALLOCATE (cdft_control%value(nvar)) 680 ALLOCATE (cdft_control%target(nvar)) 681 CALL section_vals_val_get(cdft_control_section, "STRENGTH", r_vals=rtmplist) 682 IF (SIZE(rtmplist) /= nvar) & 683 CALL cp_abort(__LOCATION__, & 684 "The length of keyword STRENGTH is incorrect. "// & 685 "Expected "//TRIM(ADJUSTL(cp_to_string(nvar)))// & 686 " value(s), got "// & 687 TRIM(ADJUSTL(cp_to_string(SIZE(rtmplist))))//" value(s).") 688 DO j = 1, nvar 689 cdft_control%strength(j) = rtmplist(j) 690 END DO 691 CALL section_vals_val_get(cdft_control_section, "TARGET", r_vals=rtmplist) 692 IF (SIZE(rtmplist) /= nvar) & 693 CALL cp_abort(__LOCATION__, & 694 "The length of keyword TARGET is incorrect. "// & 695 "Expected "//TRIM(ADJUSTL(cp_to_string(nvar)))// & 696 " value(s), got "// & 697 TRIM(ADJUSTL(cp_to_string(SIZE(rtmplist))))//" value(s).") 698 DO j = 1, nvar 699 cdft_control%target(j) = rtmplist(j) 700 END DO 701 ! Read fragment constraint definitions 702 IF (cdft_control%fragment_density) THEN 703 CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_FILE_NAME", & 704 c_val=cdft_control%fragment_a_fname) 705 CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_FILE_NAME", & 706 c_val=cdft_control%fragment_b_fname) 707 CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_SPIN_FILE", & 708 c_val=cdft_control%fragment_a_spin_fname) 709 CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_SPIN_FILE", & 710 c_val=cdft_control%fragment_b_spin_fname) 711 CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_A", & 712 l_val=cdft_control%flip_fragment(1)) 713 CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_B", & 714 l_val=cdft_control%flip_fragment(2)) 715 END IF 716 717 END SUBROUTINE read_constraint_definitions 718 719! ************************************************************************************************** 720!> \brief reads the input parameters needed for CDFT with OT 721!> \param qs_control the qs_control which holds the CDFT control type 722!> \param cdft_control_section the input section for CDFT 723!> \author Nico Holmberg [12.2015] 724! ************************************************************************************************** 725 SUBROUTINE read_cdft_control_section(qs_control, cdft_control_section) 726 TYPE(qs_control_type), INTENT(INOUT) :: qs_control 727 TYPE(section_vals_type), POINTER :: cdft_control_section 728 729 CHARACTER(len=*), PARAMETER :: routineN = 'read_cdft_control_section', & 730 routineP = moduleN//':'//routineN 731 732 LOGICAL :: exists 733 TYPE(cdft_control_type), POINTER :: cdft_control 734 TYPE(section_vals_type), POINTER :: becke_constraint_section, & 735 hirshfeld_constraint_section, & 736 outer_scf_section, print_section 737 738 NULLIFY (outer_scf_section, hirshfeld_constraint_section, becke_constraint_section, & 739 print_section) 740 cdft_control => qs_control%cdft_control 741 CPASSERT(ASSOCIATED(cdft_control)) 742 743 CALL section_vals_val_get(cdft_control_section, "TYPE_OF_CONSTRAINT", & 744 i_val=qs_control%cdft_control%type) 745 746 IF (cdft_control%type /= outer_scf_none) THEN 747 CALL section_vals_val_get(cdft_control_section, "REUSE_PRECOND", & 748 l_val=cdft_control%reuse_precond) 749 CALL section_vals_val_get(cdft_control_section, "PRECOND_FREQ", & 750 i_val=cdft_control%precond_freq) 751 CALL section_vals_val_get(cdft_control_section, "MAX_REUSE", & 752 i_val=cdft_control%max_reuse) 753 CALL section_vals_val_get(cdft_control_section, "PURGE_HISTORY", & 754 l_val=cdft_control%purge_history) 755 CALL section_vals_val_get(cdft_control_section, "PURGE_FREQ", & 756 i_val=cdft_control%purge_freq) 757 CALL section_vals_val_get(cdft_control_section, "PURGE_OFFSET", & 758 i_val=cdft_control%purge_offset) 759 CALL section_vals_val_get(cdft_control_section, "COUNTER", & 760 i_val=cdft_control%ienergy) 761 print_section => section_vals_get_subs_vals(cdft_control_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION") 762 CALL section_vals_get(print_section, explicit=cdft_control%print_weight) 763 764 outer_scf_section => section_vals_get_subs_vals(cdft_control_section, "OUTER_SCF") 765 CALL outer_scf_read_parameters(cdft_control%constraint_control, outer_scf_section) 766 IF (cdft_control%constraint_control%have_scf) THEN 767 IF (cdft_control%constraint_control%type /= outer_scf_cdft_constraint) & 768 CPABORT("Unsupported CDFT constraint.") 769 ! Constraint definitions 770 CALL read_constraint_definitions(cdft_control, cdft_control_section) 771 ! Constraint-specific initializations 772 SELECT CASE (cdft_control%type) 773 CASE (outer_scf_becke_constraint) 774 becke_constraint_section => section_vals_get_subs_vals(cdft_control_section, "BECKE_CONSTRAINT") 775 CALL section_vals_get(becke_constraint_section, explicit=exists) 776 IF (.NOT. exists) CPABORT("BECKE_CONSTRAINT section is missing.") 777 CALL read_becke_section(cdft_control, becke_constraint_section) 778 CASE (outer_scf_hirshfeld_constraint) 779 hirshfeld_constraint_section => section_vals_get_subs_vals(cdft_control_section, "HIRSHFELD_CONSTRAINT") 780 CALL section_vals_get(hirshfeld_constraint_section, explicit=exists) 781 IF (.NOT. exists) CPABORT("HIRSHFELD_CONSTRAINT section is missing.") 782 CALL read_hirshfeld_constraint_section(cdft_control, hirshfeld_constraint_section) 783 CASE DEFAULT 784 CPABORT("Unknown constraint type.") 785 END SELECT 786 787 CALL cite_reference(Holmberg2017) 788 CALL cite_reference(Holmberg2018) 789 ELSE 790 qs_control%cdft = .FALSE. 791 END IF 792 ELSE 793 qs_control%cdft = .FALSE. 794 END IF 795 796 END SUBROUTINE read_cdft_control_section 797 798! ************************************************************************************************** 799!> \brief reads the input parameters needed for Hirshfeld constraint 800!> \param cdft_control the cdft_control which holds the Hirshfeld constraint 801!> \param hirshfeld_section the input section for a Hirshfeld constraint 802!> \author Nico Holmberg [09.2018] 803! ************************************************************************************************** 804 SUBROUTINE read_hirshfeld_constraint_section(cdft_control, hirshfeld_section) 805 TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control 806 TYPE(section_vals_type), POINTER :: hirshfeld_section 807 808 CHARACTER(len=*), PARAMETER :: routineN = 'read_hirshfeld_constraint_section', & 809 routineP = moduleN//':'//routineN 810 811 LOGICAL :: exists 812 REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist 813 TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control 814 815 NULLIFY (rtmplist) 816 hirshfeld_control => cdft_control%hirshfeld_control 817 CPASSERT(ASSOCIATED(hirshfeld_control)) 818 819 CALL section_vals_val_get(hirshfeld_section, "SHAPE_FUNCTION", i_val=hirshfeld_control%shape_function) 820 CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_SHAPE", i_val=hirshfeld_control%gaussian_shape) 821 CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_RADIUS", r_val=hirshfeld_control%radius) 822 CALL section_vals_val_get(hirshfeld_section, "USE_BOHR", l_val=hirshfeld_control%use_bohr) 823 IF (.NOT. hirshfeld_control%use_bohr) THEN 824 hirshfeld_control%radius = cp_unit_from_cp2k(hirshfeld_control%radius, "angstrom") 825 END IF 826 827 IF (hirshfeld_control%shape_function == shape_function_gaussian .AND. & 828 hirshfeld_control%gaussian_shape == radius_user) THEN 829 CALL section_vals_val_get(hirshfeld_section, "ATOMIC_RADII", explicit=exists) 830 IF (.NOT. exists) CPABORT("Keyword ATOMIC_RADII is missing.") 831 CALL section_vals_val_get(hirshfeld_section, "ATOMIC_RADII", r_vals=rtmplist) 832 CPASSERT(SIZE(rtmplist) > 0) 833 ALLOCATE (hirshfeld_control%radii(SIZE(rtmplist))) 834 hirshfeld_control%radii(:) = rtmplist 835 END IF 836 837 CALL create_hirshfeld_type(hirshfeld_control%hirshfeld_env) 838 CALL set_hirshfeld_info(hirshfeld_control%hirshfeld_env, & 839 shape_function_type=hirshfeld_control%shape_function, & 840 iterative=.FALSE., & 841 radius_type=hirshfeld_control%gaussian_shape, & 842 use_bohr=hirshfeld_control%use_bohr) 843 844 END SUBROUTINE read_hirshfeld_constraint_section 845 846! ************************************************************************************************** 847!> \brief Calculate fout = fun1/fun2 or fout = fun1*fun2 848!> \param fout the output 3D potential 849!> \param fun1 the first input 3D potential 850!> \param fun2 the second input 3D potential 851!> \param divide logical that decides whether to divide or multiply the input potentials 852! ************************************************************************************************** 853 SUBROUTINE hfun_scale(fout, fun1, fun2, divide) 854 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout 855 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2 856 LOGICAL, INTENT(IN) :: divide 857 858 CHARACTER(len=*), PARAMETER :: routineN = 'hfun_scale', routineP = moduleN//':'//routineN 859 REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp 860 861 INTEGER :: i1, i2, i3, n1, n2, n3 862 863 n1 = SIZE(fout, 1) 864 n2 = SIZE(fout, 2) 865 n3 = SIZE(fout, 3) 866 CPASSERT(n1 == SIZE(fun1, 1)) 867 CPASSERT(n2 == SIZE(fun1, 2)) 868 CPASSERT(n3 == SIZE(fun1, 3)) 869 CPASSERT(n1 == SIZE(fun2, 1)) 870 CPASSERT(n2 == SIZE(fun2, 2)) 871 CPASSERT(n3 == SIZE(fun2, 3)) 872 873 IF (divide) THEN 874 DO i3 = 1, n3 875 DO i2 = 1, n2 876 DO i1 = 1, n1 877 IF (fun2(i1, i2, i3) > small) THEN 878 fout(i1, i2, i3) = fun1(i1, i2, i3)/fun2(i1, i2, i3) 879 ELSE 880 fout(i1, i2, i3) = 0.0_dp 881 END IF 882 END DO 883 END DO 884 END DO 885 ELSE 886 DO i3 = 1, n3 887 DO i2 = 1, n2 888 DO i1 = 1, n1 889 fout(i1, i2, i3) = fun1(i1, i2, i3)*fun2(i1, i2, i3) 890 END DO 891 END DO 892 END DO 893 END IF 894 895 END SUBROUTINE hfun_scale 896 897! ************************************************************************************************** 898!> \brief Determine confinement bounds along confinement dir (hardcoded to be z) 899!> and optionally zero entries below a given threshold 900!> \param fun input 3D potential (real space) 901!> \param th threshold for screening values 902!> \param just_bounds if the bounds should be computed without zeroing values 903!> \param bounds the confinement bounds: fun is nonzero only between these values along 3rd dimension 904! ************************************************************************************************** 905 SUBROUTINE hfun_zero(fun, th, just_bounds, bounds) 906 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: fun 907 REAL(KIND=dp), INTENT(IN) :: th 908 LOGICAL :: just_bounds 909 INTEGER, OPTIONAL :: bounds(2) 910 911 CHARACTER(len=*), PARAMETER :: routineN = 'hfun_zero', routineP = moduleN//':'//routineN 912 913 INTEGER :: i1, i2, i3, lb, n1, n2, n3, nzeroed, & 914 nzeroed_inner, ub 915 LOGICAL :: lb_final, ub_final 916 917 n1 = SIZE(fun, 1) 918 n2 = SIZE(fun, 2) 919 n3 = SIZE(fun, 3) 920 IF (just_bounds) THEN 921 CPASSERT(PRESENT(bounds)) 922 lb = 1 923 lb_final = .FALSE. 924 ub_final = .FALSE. 925 END IF 926 927 DO i3 = 1, n3 928 IF (just_bounds) nzeroed = 0 929 DO i2 = 1, n2 930 IF (just_bounds) nzeroed_inner = 0 931 DO i1 = 1, n1 932 IF (fun(i1, i2, i3) < th) THEN 933 IF (just_bounds) THEN 934 nzeroed_inner = nzeroed_inner + 1 935 ELSE 936 fun(i1, i2, i3) = 0.0_dp 937 END IF 938 ELSE 939 IF (just_bounds) EXIT 940 END IF 941 END DO 942 IF (just_bounds) THEN 943 IF (nzeroed_inner < n1) EXIT 944 nzeroed = nzeroed + nzeroed_inner 945 END IF 946 END DO 947 IF (just_bounds) THEN 948 IF (nzeroed == (n2*n1)) THEN 949 IF (.NOT. lb_final) THEN 950 lb = i3 951 ELSE IF (.NOT. ub_final) THEN 952 ub = i3 953 ub_final = .TRUE. 954 END IF 955 ELSE 956 IF (.NOT. lb_final) lb_final = .TRUE. 957 IF (ub_final) ub_final = .FALSE. ! Safeguard against "holes" 958 END IF 959 END IF 960 END DO 961 IF (just_bounds) THEN 962 IF (.NOT. ub_final) ub = n3 963 bounds(1) = lb 964 bounds(2) = ub 965 bounds = bounds - (n3/2) - 1 966 END IF 967 968 END SUBROUTINE hfun_zero 969 970! ************************************************************************************************** 971!> \brief Initializes Gaussian Hirshfeld constraints 972!> \param qs_env the qs_env where to build the constraint 973!> \author Nico Holmberg (09.2018) 974! ************************************************************************************************** 975 SUBROUTINE hirshfeld_constraint_init(qs_env) 976 TYPE(qs_environment_type), POINTER :: qs_env 977 978 CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_init', & 979 routineP = moduleN//':'//routineN 980 981 CHARACTER(len=2) :: element_symbol 982 INTEGER :: handle, iat, iatom, igroup, ikind, ip, & 983 iw, natom, nkind 984 INTEGER, DIMENSION(:), POINTER :: atom_list 985 REAL(KIND=dp) :: zeff 986 REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list 987 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 988 TYPE(atomic_kind_type), POINTER :: atomic_kind 989 TYPE(cdft_control_type), POINTER :: cdft_control 990 TYPE(cdft_group_type), DIMENSION(:), POINTER :: group 991 TYPE(cp_logger_type), POINTER :: logger 992 TYPE(dft_control_type), POINTER :: dft_control 993 TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control 994 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 995 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 996 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 997 TYPE(section_vals_type), POINTER :: print_section 998 999 NULLIFY (cdft_control, hirshfeld_control, hirshfeld_env, qs_kind_set, atomic_kind_set, & 1000 radii_list, dft_control, group, atomic_kind, atom_list) 1001 CALL timeset(routineN, handle) 1002 1003 logger => cp_get_default_logger() 1004 print_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") 1005 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".cdftLog") 1006 1007 CALL get_qs_env(qs_env, dft_control=dft_control) 1008 cdft_control => dft_control%qs_control%cdft_control 1009 hirshfeld_control => cdft_control%hirshfeld_control 1010 hirshfeld_env => hirshfeld_control%hirshfeld_env 1011 1012 ! Setup the Hirshfeld shape function 1013 IF (.NOT. ASSOCIATED(hirshfeld_env%kind_shape_fn)) THEN 1014 hirshfeld_env => hirshfeld_control%hirshfeld_env 1015 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) 1016 CPASSERT(ASSOCIATED(qs_kind_set)) 1017 nkind = SIZE(qs_kind_set) 1018 ! Parse atomic radii for setting up Gaussian shape function 1019 IF (ASSOCIATED(hirshfeld_control%radii)) THEN 1020 IF (.NOT. SIZE(atomic_kind_set) == SIZE(hirshfeld_control%radii)) & 1021 CALL cp_abort(__LOCATION__, & 1022 "Length of keyword HIRSHFELD_CONSTRAINT\ATOMIC_RADII does not "// & 1023 "match number of atomic kinds in the input coordinate file.") 1024 1025 ALLOCATE (radii_list(SIZE(hirshfeld_control%radii))) 1026 DO ikind = 1, SIZE(hirshfeld_control%radii) 1027 IF (hirshfeld_control%use_bohr) THEN 1028 radii_list(ikind) = hirshfeld_control%radii(ikind) 1029 ELSE 1030 radii_list(ikind) = cp_unit_from_cp2k(hirshfeld_control%radii(ikind), "angstrom") 1031 END IF 1032 END DO 1033 END IF 1034 ! radius/radii_list parameters are optional for shape_function_density 1035 CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, & 1036 radius=hirshfeld_control%radius, & 1037 radii_list=radii_list) 1038 IF (ASSOCIATED(radii_list)) DEALLOCATE (radii_list) 1039 END IF 1040 1041 ! Atomic reference charges (Mulliken not supported) 1042 IF (.NOT. ASSOCIATED(hirshfeld_env%charges)) THEN 1043 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & 1044 nkind=nkind, natom=natom) 1045 ALLOCATE (hirshfeld_env%charges(natom)) 1046 DO ikind = 1, nkind 1047 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 1048 atomic_kind => atomic_kind_set(ikind) 1049 CALL get_atomic_kind(atomic_kind, atom_list=atom_list) 1050 DO iat = 1, SIZE(atom_list) 1051 iatom = atom_list(iat) 1052 hirshfeld_env%charges(iatom) = zeff 1053 END DO 1054 END DO 1055 END IF 1056 1057 ! Print some additional information about the calculation on the first iteration 1058 IF (cdft_control%first_iteration) THEN 1059 IF (iw > 0) THEN 1060 group => cdft_control%group 1061 CALL get_qs_env(qs_env, particle_set=particle_set) 1062 IF (ASSOCIATED(hirshfeld_control%radii)) THEN 1063 WRITE (iw, '(T3,A)') & 1064 'Atom Element Gaussian radius (angstrom)' 1065 DO iatom = 1, natom 1066 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol) 1067 WRITE (iw, "(i7,T15,A2,T37,F8.3)") & 1068 iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(hirshfeld_control%radii(iatom), "angstrom") 1069 END DO 1070 WRITE (iw, '(T3,A)') & 1071 '------------------------------------------------------------------------' 1072 END IF 1073 WRITE (iw, '(/,T3,A,T60)') & 1074 '----------------------- CDFT group definitions -------------------------' 1075 DO igroup = 1, SIZE(group) 1076 IF (igroup > 1) WRITE (iw, '(T3,A)') ' ' 1077 WRITE (iw, '(T5,A,I5,A,I5)') & 1078 'Atomic group', igroup, ' of ', SIZE(group) 1079 WRITE (iw, '(T5,A)') 'Atom Element Coefficient' 1080 DO ip = 1, SIZE(group(igroup)%atoms) 1081 iatom = group(igroup)%atoms(ip) 1082 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol) 1083 WRITE (iw, '(i8,T16,A2,T23,F8.3)') iatom, ADJUSTR(element_symbol), group(igroup)%coeff(ip) 1084 END DO 1085 END DO 1086 WRITE (iw, '(T3,A)') & 1087 '------------------------------------------------------------------------' 1088 END IF 1089 cdft_control%first_iteration = .FALSE. 1090 END IF 1091 1092 ! Radii no longer needed 1093 IF (ASSOCIATED(hirshfeld_control%radii)) DEALLOCATE (hirshfeld_control%radii) 1094 CALL timestop(handle) 1095 1096 END SUBROUTINE hirshfeld_constraint_init 1097 1098! ************************************************************************************************** 1099!> \brief Prints information about CDFT constraints 1100!> \param qs_env the qs_env where to build the constraint 1101!> \param electronic_charge the CDFT charges 1102!> \par History 1103!> Created 9.2018 [Nico Holmberg] 1104! ************************************************************************************************** 1105 SUBROUTINE cdft_constraint_print(qs_env, electronic_charge) 1106 TYPE(qs_environment_type), POINTER :: qs_env 1107 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge 1108 1109 CHARACTER(len=*), PARAMETER :: routineN = 'cdft_constraint_print', & 1110 routineP = moduleN//':'//routineN 1111 1112 CHARACTER(len=2) :: element_symbol 1113 INTEGER :: iatom, ikind, iw, jatom 1114 REAL(kind=dp) :: tc(2), zeff 1115 TYPE(cdft_control_type), POINTER :: cdft_control 1116 TYPE(cp_logger_type), POINTER :: logger 1117 TYPE(dft_control_type), POINTER :: dft_control 1118 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1119 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1120 TYPE(section_vals_type), POINTER :: cdft_constraint_section 1121 1122 NULLIFY (cdft_constraint_section, logger, particle_set, dft_control, qs_kind_set) 1123 logger => cp_get_default_logger() 1124 1125 CALL get_qs_env(qs_env, & 1126 particle_set=particle_set, & 1127 dft_control=dft_control, & 1128 qs_kind_set=qs_kind_set) 1129 CPASSERT(ASSOCIATED(qs_kind_set)) 1130 1131 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") 1132 iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog") 1133 cdft_control => dft_control%qs_control%cdft_control 1134 1135 ! Print constraint information 1136 CALL qs_scf_cdft_constraint_info(iw, cdft_control) 1137 1138 ! Print weight function(s) to cube file(s) whenever weight is (re)built 1139 IF (cdft_control%print_weight .AND. cdft_control%need_pot) & 1140 CALL cdft_print_weight_function(qs_env) 1141 1142 ! Print atomic CDFT charges 1143 IF (iw > 0 .AND. cdft_control%atomic_charges) THEN 1144 IF (.NOT. cdft_control%fragment_density) THEN 1145 IF (dft_control%nspins == 1) THEN 1146 WRITE (iw, '(/,T3,A)') & 1147 '-------------------------------- CDFT atomic charges --------------------------------' 1148 WRITE (iw, '(T3,A,A)') & 1149 '#Atom Element Is_constraint', ' Core charge Population (total)'// & 1150 ' Net charge' 1151 tc = 0.0_dp 1152 DO iatom = 1, cdft_control%natoms 1153 jatom = cdft_control%atoms(iatom) 1154 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, & 1155 element_symbol=element_symbol, & 1156 kind_number=ikind) 1157 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 1158 WRITE (iw, "(i7,T15,A2,T23,L10,T39,F8.3,T61,F8.3,T81,F8.3)") & 1159 jatom, ADJUSTR(element_symbol), cdft_control%is_constraint(iatom), zeff, electronic_charge(iatom, 1), & 1160 (zeff - electronic_charge(iatom, 1)) 1161 tc(1) = tc(1) + (zeff - electronic_charge(iatom, 1)) 1162 END DO 1163 WRITE (iw, '(/,T3,A,T81,F8.3,/)') "Total Charge: ", tc(1) 1164 ELSE 1165 WRITE (iw, '(/,T3,A)') & 1166 '------------------------------------------ CDFT atomic charges -------------------------------------------' 1167 WRITE (iw, '(T3,A,A)') & 1168 '#Atom Element Is_constraint', ' Core charge Population (alpha, beta)'// & 1169 ' Net charge Spin population' 1170 tc = 0.0_dp 1171 DO iatom = 1, cdft_control%natoms 1172 jatom = cdft_control%atoms(iatom) 1173 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, & 1174 element_symbol=element_symbol, & 1175 kind_number=ikind) 1176 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 1177 WRITE (iw, "(i7,T15,A2,T23,L10,T39,F8.3,T53,F8.3,T67,F8.3,T81,F8.3,T102,F8.3)") & 1178 jatom, ADJUSTR(element_symbol), & 1179 cdft_control%is_constraint(iatom), & 1180 zeff, electronic_charge(iatom, 1), electronic_charge(iatom, 2), & 1181 (zeff - electronic_charge(iatom, 1) - electronic_charge(iatom, 2)), & 1182 electronic_charge(iatom, 1) - electronic_charge(iatom, 2) 1183 tc(1) = tc(1) + (zeff - electronic_charge(iatom, 1) - electronic_charge(iatom, 2)) 1184 tc(2) = tc(2) + (electronic_charge(iatom, 1) - electronic_charge(iatom, 2)) 1185 END DO 1186 WRITE (iw, '(/,T3,A,T81,F8.3,T102,F8.3/)') "Total Charge and Spin Moment: ", tc(1), tc(2) 1187 END IF 1188 ELSE 1189 IF (ALL(cdft_control%group(:)%constraint_type == cdft_charge_constraint)) THEN 1190 WRITE (iw, '(/,T3,A)') & 1191 '-------------------------------- CDFT atomic charges --------------------------------' 1192 IF (dft_control%nspins == 1) THEN 1193 WRITE (iw, '(T3,A,A)') & 1194 '#Atom Element Is_constraint', ' Fragment charge Population (total)'// & 1195 ' Net charge' 1196 ELSE 1197 WRITE (iw, '(T3,A,A)') & 1198 '#Atom Element Is_constraint', ' Fragment charge Population (alpha, beta)'// & 1199 ' Net charge' 1200 END IF 1201 tc = 0.0_dp 1202 DO iatom = 1, cdft_control%natoms 1203 jatom = cdft_control%atoms(iatom) 1204 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, & 1205 element_symbol=element_symbol, & 1206 kind_number=ikind) 1207 IF (dft_control%nspins == 1) THEN 1208 WRITE (iw, "(i7,T15,A2,T23,L10,T43,F8.3,T65,F8.3,T81,F8.3)") & 1209 jatom, ADJUSTR(element_symbol), & 1210 cdft_control%is_constraint(iatom), & 1211 cdft_control%charges_fragment(iatom, 1), & 1212 electronic_charge(iatom, 1), & 1213 (electronic_charge(iatom, 1) - & 1214 cdft_control%charges_fragment(iatom, 1)) 1215 tc(1) = tc(1) + (electronic_charge(iatom, 1) - & 1216 cdft_control%charges_fragment(iatom, 1)) 1217 ELSE 1218 WRITE (iw, "(i7,T15,A2,T23,L10,T43,F8.3,T57,F8.3,T69,F8.3,T81,F8.3)") & 1219 jatom, ADJUSTR(element_symbol), & 1220 cdft_control%is_constraint(iatom), & 1221 cdft_control%charges_fragment(iatom, 1), & 1222 electronic_charge(iatom, 1), electronic_charge(iatom, 2), & 1223 (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - & 1224 cdft_control%charges_fragment(iatom, 1)) 1225 tc(1) = tc(1) + (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - & 1226 cdft_control%charges_fragment(iatom, 1)) 1227 END IF 1228 END DO 1229 WRITE (iw, '(/,T3,A,T81,F8.3,/)') "Total Charge: ", tc(1) 1230 ELSE 1231 WRITE (iw, '(/,T3,A)') & 1232 '------------------------------------------ CDFT atomic charges -------------------------------------------' 1233 WRITE (iw, '(T3,A,A)') & 1234 '#Atom Element Is_constraint', ' Fragment charge/spin moment'// & 1235 ' Population (alpha, beta) Net charge/spin moment' 1236 tc = 0.0_dp 1237 DO iatom = 1, cdft_control%natoms 1238 jatom = cdft_control%atoms(iatom) 1239 CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, & 1240 element_symbol=element_symbol, & 1241 kind_number=ikind) 1242 WRITE (iw, "(i7,T15,A2,T22,L10,T40,F8.3,T52,F8.3,T66,F8.3,T78,F8.3,T90,F8.3,T102,F8.3)") & 1243 jatom, ADJUSTR(element_symbol), & 1244 cdft_control%is_constraint(iatom), & 1245 cdft_control%charges_fragment(iatom, 1), & 1246 cdft_control%charges_fragment(iatom, 2), & 1247 electronic_charge(iatom, 1), electronic_charge(iatom, 2), & 1248 (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - & 1249 cdft_control%charges_fragment(iatom, 1)), & 1250 (electronic_charge(iatom, 1) - electronic_charge(iatom, 2) - & 1251 cdft_control%charges_fragment(iatom, 2)) 1252 tc(1) = tc(1) + (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - & 1253 cdft_control%charges_fragment(iatom, 1)) 1254 tc(2) = tc(2) + (electronic_charge(iatom, 1) - electronic_charge(iatom, 2) - & 1255 cdft_control%charges_fragment(iatom, 2)) 1256 END DO 1257 WRITE (iw, '(/,T3,A,T90,F8.3,T102,F8.3/)') "Total Charge and Spin Moment: ", tc(1), tc(2) 1258 END IF 1259 END IF 1260 END IF 1261 1262 END SUBROUTINE cdft_constraint_print 1263 1264! ************************************************************************************************** 1265!> \brief Prints CDFT weight functions to cube files 1266!> \param qs_env the qs_env where to build the constraint 1267!> \par History 1268!> Created 9.2018 [Nico Holmberg] 1269! ************************************************************************************************** 1270 SUBROUTINE cdft_print_weight_function(qs_env) 1271 TYPE(qs_environment_type), POINTER :: qs_env 1272 1273 CHARACTER(len=*), PARAMETER :: routineN = 'cdft_print_weight_function', & 1274 routineP = moduleN//':'//routineN 1275 1276 CHARACTER(LEN=default_path_length) :: middle_name 1277 INTEGER :: igroup, unit_nr 1278 LOGICAL :: mpi_io 1279 TYPE(cdft_control_type), POINTER :: cdft_control 1280 TYPE(cp_logger_type), POINTER :: logger 1281 TYPE(cp_para_env_type), POINTER :: para_env 1282 TYPE(dft_control_type), POINTER :: dft_control 1283 TYPE(particle_list_type), POINTER :: particles 1284 TYPE(qs_subsys_type), POINTER :: subsys 1285 TYPE(section_vals_type), POINTER :: cdft_constraint_section 1286 1287 NULLIFY (cdft_constraint_section, logger, particles, dft_control, & 1288 para_env, subsys, cdft_control) 1289 logger => cp_get_default_logger() 1290 1291 CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control) 1292 CALL qs_subsys_get(subsys, particles=particles) 1293 cdft_control => dft_control%qs_control%cdft_control 1294 cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT") 1295 1296 DO igroup = 1, SIZE(cdft_control%group) 1297 mpi_io = .TRUE. 1298 middle_name = "cdft_weight_"//TRIM(ADJUSTL(cp_to_string(igroup))) 1299 unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", & 1300 middle_name=middle_name, & 1301 extension=".cube", file_position="REWIND", & 1302 log_filename=.FALSE., mpi_io=mpi_io) 1303 ! Note PROGRAM_RUN_INFO section neeeds to be active! 1304 IF (para_env%ionode .AND. unit_nr .LT. 1) & 1305 CALL cp_abort(__LOCATION__, & 1306 "Please turn on PROGRAM_RUN_INFO to print CDFT weight function.") 1307 1308 CALL cp_pw_to_cube(cdft_control%group(igroup)%weight%pw, & 1309 unit_nr, & 1310 "CDFT Weight Function", & 1311 particles=particles, & 1312 stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"), & 1313 mpi_io=mpi_io) 1314 CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io) 1315 END DO 1316 1317 END SUBROUTINE cdft_print_weight_function 1318 1319END MODULE qs_cdft_utils 1320