1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines to handle an external electrostatic field 8!> The external field can be generic and is provided by user input 9! ************************************************************************************************** 10MODULE qs_external_potential 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind 13 USE cell_types, ONLY: cell_type,& 14 pbc 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_log_handling, ONLY: cp_to_string 17 USE cp_realspace_grid_cube, ONLY: cp_cube_to_pw 18 USE force_fields_util, ONLY: get_generic_info 19 USE fparser, ONLY: evalf,& 20 evalfd,& 21 finalizef,& 22 initf,& 23 parsef 24 USE input_section_types, ONLY: section_vals_get_subs_vals,& 25 section_vals_type,& 26 section_vals_val_get 27 USE kinds, ONLY: default_path_length,& 28 default_string_length,& 29 dp 30 USE message_passing, ONLY: mp_bcast 31 USE particle_types, ONLY: particle_type 32 USE pw_grid_types, ONLY: PW_MODE_LOCAL 33 USE pw_types, ONLY: pw_p_type,& 34 pw_type 35 USE qs_energy_types, ONLY: qs_energy_type 36 USE qs_environment_types, ONLY: get_qs_env,& 37 qs_environment_type 38 USE qs_force_types, ONLY: qs_force_type 39 USE qs_kind_types, ONLY: get_qs_kind,& 40 qs_kind_type 41 USE string_utilities, ONLY: compress 42#include "./base/base_uses.f90" 43 44 IMPLICIT NONE 45 46 PRIVATE 47 48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_external_potential' 49 50! *** Public subroutines *** 51 PUBLIC :: external_e_potential, & 52 external_c_potential 53 54CONTAINS 55 56! ************************************************************************************************** 57!> \brief Computes the external potential on the grid 58!> \param qs_env ... 59!> \date 12.2009 60!> \author Teodoro Laino [tlaino] 61! ************************************************************************************************** 62 SUBROUTINE external_e_potential(qs_env) 63 64 TYPE(qs_environment_type), POINTER :: qs_env 65 66 CHARACTER(len=*), PARAMETER :: routineN = 'external_e_potential', & 67 routineP = moduleN//':'//routineN 68 69 INTEGER :: handle, i, j, k 70 INTEGER, DIMENSION(2, 3) :: bo_global, bo_local 71 LOGICAL :: read_from_cube, static_potential 72 REAL(kind=dp) :: dvol, efunc, scaling_factor 73 REAL(kind=dp), DIMENSION(3) :: dr, grid_p 74 TYPE(dft_control_type), POINTER :: dft_control 75 TYPE(pw_p_type), POINTER :: v_ee 76 TYPE(section_vals_type), POINTER :: ext_pot_section, input 77 78 CALL timeset(routineN, handle) 79 NULLIFY (v_ee, input, ext_pot_section, dft_control) 80 CALL get_qs_env(qs_env, & 81 vee=v_ee, & 82 input=input, & 83 dft_control=dft_control) 84 IF (dft_control%apply_external_potential) THEN 85 ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL") 86 CALL section_vals_val_get(ext_pot_section, "STATIC", l_val=static_potential) 87 CALL section_vals_val_get(ext_pot_section, "READ_FROM_CUBE", l_val=read_from_cube) 88 89 IF ((.NOT. static_potential) .OR. dft_control%eval_external_potential) THEN 90 IF (read_from_cube) THEN 91 CALL section_vals_val_get(ext_pot_section, "SCALING_FACTOR", r_val=scaling_factor) 92 CALL cp_cube_to_pw(v_ee%pw, 'pot.cube', scaling_factor) 93 dft_control%eval_external_potential = .FALSE. 94 ELSE 95 dr = v_ee%pw%pw_grid%dr 96 dvol = v_ee%pw%pw_grid%dvol 97 v_ee%pw%cr3d = 0.0_dp 98 99 bo_local = v_ee%pw%pw_grid%bounds_local 100 bo_global = v_ee%pw%pw_grid%bounds 101 102 DO k = bo_local(1, 3), bo_local(2, 3) 103 DO j = bo_local(1, 2), bo_local(2, 2) 104 DO i = bo_local(1, 1), bo_local(2, 1) 105 grid_p(1) = (i - bo_global(1, 1))*dr(1) 106 grid_p(2) = (j - bo_global(1, 2))*dr(2) 107 grid_p(3) = (k - bo_global(1, 3))*dr(3) 108 CALL get_external_potential(grid_p, ext_pot_section, func=efunc) 109 v_ee%pw%cr3d(i, j, k) = v_ee%pw%cr3d(i, j, k) + efunc 110 END DO 111 END DO 112 END DO 113 dft_control%eval_external_potential = .FALSE. 114 END IF 115 END IF 116 END IF 117 CALL timestop(handle) 118 END SUBROUTINE external_e_potential 119 120! ************************************************************************************************** 121!> \brief Computes the force and the energy due to the external potential on the cores 122!> \param qs_env ... 123!> \param calculate_forces ... 124!> \date 12.2009 125!> \author Teodoro Laino [tlaino] 126! ************************************************************************************************** 127 SUBROUTINE external_c_potential(qs_env, calculate_forces) 128 129 TYPE(qs_environment_type), POINTER :: qs_env 130 LOGICAL, OPTIONAL :: calculate_forces 131 132 CHARACTER(len=*), PARAMETER :: routineN = 'external_c_potential', & 133 routineP = moduleN//':'//routineN 134 135 INTEGER :: atom_a, handle, iatom, ikind, natom, & 136 nkind 137 INTEGER, DIMENSION(:), POINTER :: list 138 LOGICAL :: my_force, pot_on_grid 139 REAL(KIND=dp) :: ee_core_ener, efunc, zeff 140 REAL(KIND=dp), DIMENSION(3) :: dfunc, r 141 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 142 TYPE(cell_type), POINTER :: cell 143 TYPE(dft_control_type), POINTER :: dft_control 144 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 145 TYPE(pw_p_type), POINTER :: v_ee 146 TYPE(qs_energy_type), POINTER :: energy 147 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 148 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 149 TYPE(section_vals_type), POINTER :: ext_pot_section, input 150 151 CALL timeset(routineN, handle) 152 NULLIFY (dft_control) 153 154 CALL get_qs_env(qs_env=qs_env, & 155 atomic_kind_set=atomic_kind_set, & 156 qs_kind_set=qs_kind_set, & 157 energy=energy, & 158 particle_set=particle_set, & 159 input=input, & 160 cell=cell, & 161 dft_control=dft_control) 162 163 IF (dft_control%apply_external_potential) THEN 164 IF (dft_control%eval_external_potential) THEN !ensure that external potential is loaded to grid 165 CALL external_e_potential(qs_env) 166 END IF 167 my_force = .FALSE. 168 IF (PRESENT(calculate_forces)) my_force = calculate_forces 169 ext_pot_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_POTENTIAL") 170 ee_core_ener = 0.0_dp 171 nkind = SIZE(atomic_kind_set) 172 173 !check if external potential on grid has been loaded from a file instead of giving a function 174 CALL section_vals_val_get(ext_pot_section, "READ_FROM_CUBE", l_val=pot_on_grid) 175 IF (pot_on_grid) CALL get_qs_env(qs_env, vee=v_ee, input=input) 176 177 DO ikind = 1, SIZE(atomic_kind_set) 178 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom) 179 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 180 181 natom = SIZE(list) 182 DO iatom = 1, natom 183 atom_a = list(iatom) 184 !pbc returns r(i) in range [-cell%hmat(i,i)/2, cell%hmat(i,i)/2] 185 !for periodic dimensions (assuming the cell is orthorombic). 186 !This is not consistent with the potential on grid, where r(i) is 187 !in range [0, cell%hmat(i,i)] 188 !Use new pbc function with switch positive_range=.TRUE. 189 r(:) = pbc(particle_set(atom_a)%r(:), cell, positive_range=.TRUE.) 190 !if potential is on grid, interpolate the value at r, 191 !otherwise evaluate the given function 192 IF (pot_on_grid) THEN 193 CALL interpolate_external_potential(r, v_ee%pw, func=efunc, & 194 dfunc=dfunc, calc_derivatives=my_force) 195 ELSE 196 CALL get_external_potential(r, ext_pot_section, func=efunc, & 197 dfunc=dfunc, calc_derivatives=my_force) 198 END IF 199 ee_core_ener = ee_core_ener + zeff*efunc 200 IF (my_force) THEN 201 CALL get_qs_env(qs_env=qs_env, force=force) 202 force(ikind)%eev(:, iatom) = dfunc*zeff 203 END IF 204 END DO 205 END DO 206 energy%ee_core = ee_core_ener 207 END IF 208 CALL timestop(handle) 209 END SUBROUTINE external_c_potential 210 211! ************************************************************************************************** 212!> \brief Low level function for computing the potential and the derivatives 213!> \param r position in realspace 214!> \param ext_pot_section ... 215!> \param func external potential at r 216!> \param dfunc derivative of the external potential at r 217!> \param calc_derivatives Whether to calculate dfunc 218!> \date 12.2009 219!> \par History 220!> 12.2009 created [tlaino] 221!> 11.2014 reading external cube file added [Juha Ritala & Matt Watkins] 222!> \author Teodoro Laino [tlaino] 223! ************************************************************************************************** 224 SUBROUTINE get_external_potential(r, ext_pot_section, func, dfunc, calc_derivatives) 225 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 226 TYPE(section_vals_type), POINTER :: ext_pot_section 227 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: func, dfunc(3) 228 LOGICAL, INTENT(IN), OPTIONAL :: calc_derivatives 229 230 CHARACTER(len=*), PARAMETER :: routineN = 'get_external_potential', & 231 routineP = moduleN//':'//routineN 232 233 CHARACTER(LEN=default_path_length) :: coupling_function 234 CHARACTER(LEN=default_string_length) :: def_error, this_error 235 CHARACTER(LEN=default_string_length), & 236 DIMENSION(:), POINTER :: my_par 237 INTEGER :: handle, j 238 LOGICAL :: check, my_force 239 REAL(KIND=dp) :: dedf, dx, err, lerr 240 REAL(KIND=dp), DIMENSION(:), POINTER :: my_val 241 242 CALL timeset(routineN, handle) 243 NULLIFY (my_par, my_val) 244 my_force = .FALSE. 245 IF (PRESENT(calc_derivatives)) my_force = calc_derivatives 246 check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives) 247 CPASSERT(check) 248 CALL section_vals_val_get(ext_pot_section, "DX", r_val=dx) 249 CALL section_vals_val_get(ext_pot_section, "ERROR_LIMIT", r_val=lerr) 250 CALL get_generic_info(ext_pot_section, "FUNCTION", coupling_function, my_par, my_val, & 251 input_variables=(/"X", "Y", "Z"/), i_rep_sec=1) 252 CALL initf(1) 253 CALL parsef(1, TRIM(coupling_function), my_par) 254 255 my_val(1) = r(1) 256 my_val(2) = r(2) 257 my_val(3) = r(3) 258 259 IF (PRESENT(func)) func = evalf(1, my_val) 260 IF (my_force) THEN 261 DO j = 1, 3 262 dedf = evalfd(1, j, my_val, dx, err) 263 IF (ABS(err) > lerr) THEN 264 WRITE (this_error, "(A,G12.6,A)") "(", err, ")" 265 WRITE (def_error, "(A,G12.6,A)") "(", lerr, ")" 266 CALL compress(this_error, .TRUE.) 267 CALL compress(def_error, .TRUE.) 268 CALL cp_warn(__LOCATION__, & 269 'ASSERTION (cond) failed at line '//cp_to_string(__LINE__)// & 270 ' Error '//TRIM(this_error)//' in computing numerical derivatives larger then'// & 271 TRIM(def_error)//' .') 272 END IF 273 dfunc(j) = dedf 274 END DO 275 END IF 276 DEALLOCATE (my_par) 277 DEALLOCATE (my_val) 278 CALL finalizef() 279 CALL timestop(handle) 280 END SUBROUTINE get_external_potential 281 282! ************************************************************************************************** 283!> \brief subroutine that interpolates the value of the external 284!> potential at position r based on the values on the realspace grid 285!> \param r ... 286!> \param grid external potential pw grid, vee 287!> \param func value of vee at r 288!> \param dfunc derivatives of vee at r 289!> \param calc_derivatives calc dfunc 290! ************************************************************************************************** 291 SUBROUTINE interpolate_external_potential(r, grid, func, dfunc, calc_derivatives) 292 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 293 TYPE(pw_type), POINTER :: grid 294 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: func, dfunc(3) 295 LOGICAL, INTENT(IN), OPTIONAL :: calc_derivatives 296 297 CHARACTER(len=*), PARAMETER :: routineN = 'interpolate_external_potential', & 298 routineP = moduleN//':'//routineN 299 300 INTEGER :: buffer_i, buffer_j, buffer_k, & 301 data_source, fd_extra_point, gid, & 302 handle, i, i_pbc, ip, j, j_pbc, k, & 303 k_pbc, my_rank, num_pe, tag 304 INTEGER, DIMENSION(3) :: lbounds, lbounds_local, lower_inds, & 305 ubounds, ubounds_local, upper_inds 306 LOGICAL :: check, my_force 307 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: bcast_buffer 308 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: grid_buffer 309 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dgrid 310 REAL(KIND=dp), DIMENSION(3) :: dr, subgrid_origin 311 312 CALL timeset(routineN, handle) 313 my_force = .FALSE. 314 IF (PRESENT(calc_derivatives)) my_force = calc_derivatives 315 check = PRESENT(dfunc) .EQV. PRESENT(calc_derivatives) 316 CPASSERT(check) 317 318 IF (my_force) THEN 319 ALLOCATE (grid_buffer(0:3, 0:3, 0:3)) 320 ALLOCATE (bcast_buffer(0:3)) 321 ALLOCATE (dgrid(1:2, 1:2, 1:2, 3)) 322 fd_extra_point = 1 323 ELSE 324 ALLOCATE (grid_buffer(1:2, 1:2, 1:2)) 325 ALLOCATE (bcast_buffer(1:2)) 326 fd_extra_point = 0 327 END IF 328 329 ! The values of external potential on grid are distributed among the 330 ! processes, so first we have to gather them up 331 gid = grid%pw_grid%para%group 332 my_rank = grid%pw_grid%para%my_pos 333 num_pe = grid%pw_grid%para%group_size 334 tag = 1 335 336 dr = grid%pw_grid%dr 337 lbounds = grid%pw_grid%bounds(1, :) 338 ubounds = grid%pw_grid%bounds(2, :) 339 lbounds_local = grid%pw_grid%bounds_local(1, :) 340 ubounds_local = grid%pw_grid%bounds_local(2, :) 341 342 ! Determine the indices of grid points that are needed 343 lower_inds = lbounds + FLOOR(r/dr) - fd_extra_point 344 upper_inds = lower_inds + 1 + 2*fd_extra_point 345 346 DO i = lower_inds(1), upper_inds(1) 347 ! If index is out of global bounds, assume periodic boundary conditions 348 i_pbc = pbc_index(i, lbounds(1), ubounds(1)) 349 buffer_i = i - lower_inds(1) + 1 - fd_extra_point 350 DO j = lower_inds(2), upper_inds(2) 351 j_pbc = pbc_index(j, lbounds(2), ubounds(2)) 352 buffer_j = j - lower_inds(2) + 1 - fd_extra_point 353 354 ! Find the process that has the data for indices i_pbc and j_pbc 355 ! and store the data to bcast_buffer. Assuming that each process has full z data 356 IF (grid%pw_grid%para%mode .NE. PW_MODE_LOCAL) THEN 357 DO ip = 0, num_pe - 1 358 IF (grid%pw_grid%para%bo(1, 1, ip, 1) <= i_pbc - lbounds(1) + 1 .AND. & 359 grid%pw_grid%para%bo(2, 1, ip, 1) >= i_pbc - lbounds(1) + 1 .AND. & 360 grid%pw_grid%para%bo(1, 2, ip, 1) <= j_pbc - lbounds(2) + 1 .AND. & 361 grid%pw_grid%para%bo(2, 2, ip, 1) >= j_pbc - lbounds(2) + 1) THEN 362 data_source = ip 363 EXIT 364 END IF 365 END DO 366 IF (my_rank == data_source) THEN 367 IF (lower_inds(3) >= lbounds(3) .AND. upper_inds(3) <= ubounds(3)) THEN 368 bcast_buffer(:) = & 369 grid%cr3d(i_pbc, j_pbc, lower_inds(3):upper_inds(3)) 370 ELSE 371 DO k = lower_inds(3), upper_inds(3) 372 k_pbc = pbc_index(k, lbounds(3), ubounds(3)) 373 buffer_k = k - lower_inds(3) + 1 - fd_extra_point 374 bcast_buffer(buffer_k) = & 375 grid%cr3d(i_pbc, j_pbc, k_pbc) 376 END DO 377 END IF 378 END IF 379 ! data_source sends data to everyone 380 CALL mp_bcast(bcast_buffer, data_source, gid) 381 grid_buffer(buffer_i, buffer_j, :) = bcast_buffer 382 ELSE 383 grid_buffer(buffer_i, buffer_j, :) = grid%cr3d(i_pbc, j_pbc, lower_inds(3):upper_inds(3)) 384 END IF 385 END DO 386 END DO 387 388 ! Now that all the processes have local external potential data around r, 389 ! interpolate the value at r 390 subgrid_origin = (lower_inds - lbounds + fd_extra_point)*dr 391 func = trilinear_interpolation(r, grid_buffer(1:2, 1:2, 1:2), subgrid_origin, dr) 392 393 ! If the derivative of the potential is needed, approximate the derivative at grid 394 ! points using finite differences, and then interpolate the value at r 395 IF (my_force) THEN 396 CALL d_finite_difference(grid_buffer, dr, dgrid) 397 DO i = 1, 3 398 dfunc(i) = trilinear_interpolation(r, dgrid(:, :, :, i), subgrid_origin, dr) 399 END DO 400 DEALLOCATE (dgrid) 401 END IF 402 403 DEALLOCATE (grid_buffer) 404 CALL timestop(handle) 405 END SUBROUTINE interpolate_external_potential 406 407! ************************************************************************************************** 408!> \brief subroutine that uses finite differences to approximate the partial 409!> derivatives of the potential based on the given values on grid 410!> \param grid tiny bit of external potential vee 411!> \param dr step size for finite difference 412!> \param dgrid derivatives of grid 413! ************************************************************************************************** 414 SUBROUTINE d_finite_difference(grid, dr, dgrid) 415 REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(IN) :: grid 416 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: dr 417 REAL(KIND=dp), DIMENSION(1:, 1:, 1:, :), & 418 INTENT(OUT) :: dgrid 419 420 INTEGER :: i, j, k 421 422 DO i = 1, SIZE(dgrid, 1) 423 DO j = 1, SIZE(dgrid, 2) 424 DO k = 1, SIZE(dgrid, 3) 425 dgrid(i, j, k, 1) = 0.5*(grid(i + 1, j, k) - grid(i - 1, j, k))/dr(1) 426 dgrid(i, j, k, 2) = 0.5*(grid(i, j + 1, k) - grid(i, j - 1, k))/dr(2) 427 dgrid(i, j, k, 3) = 0.5*(grid(i, j, k + 1) - grid(i, j, k - 1))/dr(3) 428 END DO 429 END DO 430 END DO 431 END SUBROUTINE d_finite_difference 432 433! ************************************************************************************************** 434!> \brief trilinear interpolation function that interpolates value at r based 435!> on 2x2x2 grid points around r in subgrid 436!> \param r where to interpolate to 437!> \param subgrid part of external potential on a grid 438!> \param origin center of grid 439!> \param dr step size 440!> \return interpolated value of external potential 441! ************************************************************************************************** 442 FUNCTION trilinear_interpolation(r, subgrid, origin, dr) RESULT(value_at_r) 443 REAL(KIND=dp), DIMENSION(3) :: r 444 REAL(KIND=dp), DIMENSION(:, :, :) :: subgrid 445 REAL(KIND=dp), DIMENSION(3) :: origin, dr 446 REAL(KIND=dp) :: value_at_r 447 448 REAL(KIND=dp), DIMENSION(3) :: norm_r, norm_r_rev 449 450 norm_r = (r - origin)/dr 451 norm_r_rev = 1 - norm_r 452 value_at_r = subgrid(1, 1, 1)*PRODUCT(norm_r_rev) + & 453 subgrid(2, 1, 1)*norm_r(1)*norm_r_rev(2)*norm_r_rev(3) + & 454 subgrid(1, 2, 1)*norm_r_rev(1)*norm_r(2)*norm_r_rev(3) + & 455 subgrid(1, 1, 2)*norm_r_rev(1)*norm_r_rev(2)*norm_r(3) + & 456 subgrid(1, 2, 2)*norm_r_rev(1)*norm_r(2)*norm_r(3) + & 457 subgrid(2, 1, 2)*norm_r(1)*norm_r_rev(2)*norm_r(3) + & 458 subgrid(2, 2, 1)*norm_r(1)*norm_r(2)*norm_r_rev(3) + & 459 subgrid(2, 2, 2)*PRODUCT(norm_r) 460 END FUNCTION trilinear_interpolation 461 462! ************************************************************************************************** 463!> \brief get a correct value for possible out of bounds index using periodic 464!> boundary conditions 465!> \param i ... 466!> \param lowbound ... 467!> \param upbound ... 468!> \return ... 469! ************************************************************************************************** 470 FUNCTION pbc_index(i, lowbound, upbound) 471 INTEGER :: i, lowbound, upbound, pbc_index 472 473 IF (i < lowbound) THEN 474 pbc_index = upbound + i - lowbound + 1 475 ELSE IF (i > upbound) THEN 476 pbc_index = lowbound + i - upbound - 1 477 ELSE 478 pbc_index = i 479 END IF 480 END FUNCTION pbc_index 481 482END MODULE qs_external_potential 483