46 #include "./base/base_uses.f90"
52 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_external_potential'
68 TYPE(qs_environment_type),
POINTER :: qs_env
70 CHARACTER(len=*),
PARAMETER :: routinen =
'external_e_potential'
72 INTEGER :: handle, i, j, k
73 INTEGER(kind=int_8) :: npoints
74 INTEGER,
DIMENSION(2, 3) :: bo_global, bo_local
75 REAL(kind=
dp) :: dvol, scaling_factor
76 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: efunc, grid_p_i, grid_p_j, grid_p_k
77 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: grid_p
78 REAL(kind=
dp),
DIMENSION(3) :: dr
79 TYPE(dft_control_type),
POINTER :: dft_control
80 TYPE(pw_r3d_rs_type),
POINTER :: v_ee
81 TYPE(section_vals_type),
POINTER :: ext_pot_section, input
83 CALL timeset(routinen, handle)
84 CALL get_qs_env(qs_env, dft_control=dft_control)
85 IF (dft_control%apply_external_potential)
THEN
86 IF (dft_control%eval_external_potential)
THEN
88 IF (dft_control%expot_control%maxwell_solver)
THEN
89 scaling_factor = dft_control%expot_control%scaling_factor
91 qs_env%sim_step, qs_env%sim_time, &
93 dft_control%eval_external_potential = .false.
94 ELSEIF (dft_control%expot_control%read_from_cube)
THEN
95 scaling_factor = dft_control%expot_control%scaling_factor
97 dft_control%eval_external_potential = .false.
103 dvol = v_ee%pw_grid%dvol
106 bo_local = v_ee%pw_grid%bounds_local
107 bo_global = v_ee%pw_grid%bounds
109 npoints = int(bo_local(2, 1) - bo_local(1, 1) + 1, kind=
int_8)* &
110 int(bo_local(2, 2) - bo_local(1, 2) + 1, kind=
int_8)* &
111 int(bo_local(2, 3) - bo_local(1, 3) + 1, kind=
int_8)
112 ALLOCATE (efunc(npoints))
113 ALLOCATE (grid_p(3, npoints))
114 ALLOCATE (grid_p_i(bo_local(1, 1):bo_local(2, 1)))
115 ALLOCATE (grid_p_j(bo_local(1, 2):bo_local(2, 2)))
116 ALLOCATE (grid_p_k(bo_local(1, 3):bo_local(2, 3)))
118 DO i = bo_local(1, 1), bo_local(2, 1)
119 grid_p_i(i) = (i - bo_global(1, 1))*dr(1)
121 DO j = bo_local(1, 2), bo_local(2, 2)
122 grid_p_j(j) = (j - bo_global(1, 2))*dr(2)
124 DO k = bo_local(1, 3), bo_local(2, 3)
125 grid_p_k(k) = (k - bo_global(1, 3))*dr(3)
129 DO k = bo_local(1, 3), bo_local(2, 3)
130 DO j = bo_local(1, 2), bo_local(2, 2)
131 DO i = bo_local(1, 1), bo_local(2, 1)
132 npoints = npoints + 1
133 grid_p(1, npoints) = grid_p_i(i)
134 grid_p(2, npoints) = grid_p_j(j)
135 grid_p(3, npoints) = grid_p_k(k)
140 DEALLOCATE (grid_p_i, grid_p_j, grid_p_k)
142 CALL get_external_potential(grid_p, ext_pot_section, func=efunc)
145 DO k = bo_local(1, 3), bo_local(2, 3)
146 DO j = bo_local(1, 2), bo_local(2, 2)
147 DO i = bo_local(1, 1), bo_local(2, 1)
148 npoints = npoints + 1
149 v_ee%array(i, j, k) = v_ee%array(i, j, k) + efunc(npoints)
154 DEALLOCATE (grid_p, efunc)
156 dft_control%eval_external_potential = .false.
160 CALL timestop(handle)
172 TYPE(qs_environment_type),
POINTER :: qs_env
173 LOGICAL,
OPTIONAL :: calculate_forces
175 CHARACTER(len=*),
PARAMETER :: routinen =
'external_c_potential'
177 INTEGER :: atom_a, handle, iatom, ikind, natom, &
179 INTEGER,
DIMENSION(:),
POINTER ::
list
180 LOGICAL :: my_force, pot_on_grid
181 REAL(kind=
dp) :: ee_core_ener, zeff
182 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: efunc
183 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dfunc, r
184 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
185 TYPE(cell_type),
POINTER :: cell
186 TYPE(dft_control_type),
POINTER :: dft_control
187 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
188 TYPE(pw_r3d_rs_type),
POINTER :: v_ee
189 TYPE(qs_energy_type),
POINTER :: energy
190 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
191 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
192 TYPE(section_vals_type),
POINTER :: ext_pot_section, input
194 CALL timeset(routinen, handle)
195 NULLIFY (dft_control)
198 atomic_kind_set=atomic_kind_set, &
199 qs_kind_set=qs_kind_set, &
201 particle_set=particle_set, &
204 dft_control=dft_control)
206 IF (dft_control%apply_external_potential)
THEN
208 IF (dft_control%eval_external_potential)
THEN
212 IF (
PRESENT(calculate_forces)) my_force = calculate_forces
213 ee_core_ener = 0.0_dp
214 nkind =
SIZE(atomic_kind_set)
217 IF (dft_control%expot_control%read_from_cube .OR. &
218 dft_control%expot_control%maxwell_solver)
THEN
222 pot_on_grid = .false.
227 DO ikind = 1,
SIZE(atomic_kind_set)
229 nparticles = nparticles + max(natom, 0)
232 ALLOCATE (efunc(nparticles))
233 ALLOCATE (dfunc(3, nparticles), r(3, nparticles))
236 DO ikind = 1,
SIZE(atomic_kind_set)
241 nparticles = nparticles + 1
247 r(:, nparticles) =
pbc(particle_set(atom_a)%r(:), cell, positive_range=.true.)
253 IF (pot_on_grid)
THEN
254 DO iatom = 1, nparticles
255 CALL interpolate_external_potential(r(:, iatom), v_ee, func=efunc(iatom), &
256 dfunc=dfunc(:, iatom), calc_derivatives=my_force)
259 CALL get_external_potential(r, ext_pot_section, func=efunc, dfunc=dfunc, calc_derivatives=my_force)
267 DO ikind = 1,
SIZE(atomic_kind_set)
272 nparticles = nparticles + 1
274 ee_core_ener = ee_core_ener + zeff*efunc(nparticles)
276 force(ikind)%eev(1:3, iatom) = dfunc(1:3, nparticles)*zeff
280 energy%ee_core = ee_core_ener
282 DEALLOCATE (dfunc, r)
285 CALL timestop(handle)
301 SUBROUTINE get_external_potential(r, ext_pot_section, func, dfunc, calc_derivatives)
302 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: r
303 TYPE(section_vals_type),
POINTER :: ext_pot_section
304 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT),
OPTIONAL :: func
305 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT), &
307 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_derivatives
309 CHARACTER(len=*),
PARAMETER :: routinen =
'get_external_potential'
311 CHARACTER(LEN=default_path_length) :: coupling_function
312 CHARACTER(LEN=default_string_length) :: def_error, this_error
313 CHARACTER(LEN=default_string_length), &
314 DIMENSION(:),
POINTER :: my_par
316 INTEGER(kind=int_8) :: ipoint, npoints
317 LOGICAL :: check, my_force
318 REAL(kind=
dp) :: dedf, dx, err, lerr
319 REAL(kind=
dp),
DIMENSION(:),
POINTER :: my_val
321 CALL timeset(routinen, handle)
322 NULLIFY (my_par, my_val)
324 IF (
PRESENT(calc_derivatives)) my_force = calc_derivatives
325 check =
PRESENT(dfunc) .EQV.
PRESENT(calc_derivatives)
329 CALL get_generic_info(ext_pot_section,
"FUNCTION", coupling_function, my_par, my_val, &
330 input_variables=(/
"X",
"Y",
"Z"/), i_rep_sec=1)
332 CALL parsef(1, trim(coupling_function), my_par)
334 npoints =
SIZE(r, 2, kind=
int_8)
336 DO ipoint = 1, npoints
337 my_val(1) = r(1, ipoint)
338 my_val(2) = r(2, ipoint)
339 my_val(3) = r(3, ipoint)
341 IF (
PRESENT(func)) func(ipoint) =
evalf(1, my_val)
344 dedf =
evalfd(1, j, my_val, dx, err)
345 IF (abs(err) > lerr)
THEN
346 WRITE (this_error,
"(A,G12.6,A)")
"(", err,
")"
347 WRITE (def_error,
"(A,G12.6,A)")
"(", lerr,
")"
350 CALL cp_warn(__location__, &
351 'ASSERTION (cond) failed at line '//cp_to_string(__line__)// &
352 ' Error '//trim(this_error)//
' in computing numerical derivatives larger then'// &
353 trim(def_error)//
' .')
355 dfunc(j, ipoint) = dedf
362 CALL timestop(handle)
363 END SUBROUTINE get_external_potential
374 SUBROUTINE interpolate_external_potential(r, grid, func, dfunc, calc_derivatives)
375 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: r
376 TYPE(pw_r3d_rs_type),
POINTER :: grid
377 REAL(kind=
dp),
INTENT(OUT),
OPTIONAL :: func, dfunc(3)
378 LOGICAL,
INTENT(IN),
OPTIONAL :: calc_derivatives
380 CHARACTER(len=*),
PARAMETER :: routinen =
'interpolate_external_potential'
382 INTEGER :: buffer_i, buffer_j, buffer_k, &
383 data_source, fd_extra_point, handle, &
384 i, i_pbc, ip, j, j_pbc, k, k_pbc, &
386 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, lower_inds, &
387 ubounds, ubounds_local, upper_inds
388 LOGICAL :: check, my_force
389 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: bcast_buffer
390 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: grid_buffer
391 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: dgrid
392 REAL(kind=
dp),
DIMENSION(3) :: dr, subgrid_origin
393 TYPE(mp_comm_type) :: gid
395 CALL timeset(routinen, handle)
397 IF (
PRESENT(calc_derivatives)) my_force = calc_derivatives
398 check =
PRESENT(dfunc) .EQV.
PRESENT(calc_derivatives)
402 ALLOCATE (grid_buffer(0:3, 0:3, 0:3))
403 ALLOCATE (bcast_buffer(0:3))
404 ALLOCATE (dgrid(1:2, 1:2, 1:2, 3))
407 ALLOCATE (grid_buffer(1:2, 1:2, 1:2))
408 ALLOCATE (bcast_buffer(1:2))
414 gid = grid%pw_grid%para%group
415 my_rank = grid%pw_grid%para%my_pos
416 num_pe = grid%pw_grid%para%group_size
420 lbounds = grid%pw_grid%bounds(1, :)
421 ubounds = grid%pw_grid%bounds(2, :)
422 lbounds_local = grid%pw_grid%bounds_local(1, :)
423 ubounds_local = grid%pw_grid%bounds_local(2, :)
426 lower_inds = lbounds + floor(r/dr) - fd_extra_point
427 upper_inds = lower_inds + 1 + 2*fd_extra_point
429 DO i = lower_inds(1), upper_inds(1)
431 i_pbc = pbc_index(i, lbounds(1), ubounds(1))
432 buffer_i = i - lower_inds(1) + 1 - fd_extra_point
433 DO j = lower_inds(2), upper_inds(2)
434 j_pbc = pbc_index(j, lbounds(2), ubounds(2))
435 buffer_j = j - lower_inds(2) + 1 - fd_extra_point
440 DO ip = 0, num_pe - 1
441 IF (grid%pw_grid%para%bo(1, 1, ip, 1) <= i_pbc - lbounds(1) + 1 .AND. &
442 grid%pw_grid%para%bo(2, 1, ip, 1) >= i_pbc - lbounds(1) + 1 .AND. &
443 grid%pw_grid%para%bo(1, 2, ip, 1) <= j_pbc - lbounds(2) + 1 .AND. &
444 grid%pw_grid%para%bo(2, 2, ip, 1) >= j_pbc - lbounds(2) + 1)
THEN
449 IF (my_rank == data_source)
THEN
450 IF (lower_inds(3) >= lbounds(3) .AND. upper_inds(3) <= ubounds(3))
THEN
452 grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
454 DO k = lower_inds(3), upper_inds(3)
455 k_pbc = pbc_index(k, lbounds(3), ubounds(3))
456 buffer_k = k - lower_inds(3) + 1 - fd_extra_point
457 bcast_buffer(buffer_k) = &
458 grid%array(i_pbc, j_pbc, k_pbc)
463 CALL gid%bcast(bcast_buffer, data_source)
464 grid_buffer(buffer_i, buffer_j, :) = bcast_buffer
466 grid_buffer(buffer_i, buffer_j, :) = grid%array(i_pbc, j_pbc, lower_inds(3):upper_inds(3))
473 subgrid_origin = (lower_inds - lbounds + fd_extra_point)*dr
474 func = trilinear_interpolation(r, grid_buffer(1:2, 1:2, 1:2), subgrid_origin, dr)
479 CALL d_finite_difference(grid_buffer, dr, dgrid)
481 dfunc(i) = trilinear_interpolation(r, dgrid(:, :, :, i), subgrid_origin, dr)
486 DEALLOCATE (grid_buffer)
487 CALL timestop(handle)
488 END SUBROUTINE interpolate_external_potential
497 PURE SUBROUTINE d_finite_difference(grid, dr, dgrid)
498 REAL(kind=
dp),
DIMENSION(0:, 0:, 0:),
INTENT(IN) :: grid
499 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: dr
500 REAL(kind=
dp),
DIMENSION(1:, 1:, 1:, :), &
505 DO i = 1,
SIZE(dgrid, 1)
506 DO j = 1,
SIZE(dgrid, 2)
507 DO k = 1,
SIZE(dgrid, 3)
508 dgrid(i, j, k, 1) = 0.5*(grid(i + 1, j, k) - grid(i - 1, j, k))/dr(1)
509 dgrid(i, j, k, 2) = 0.5*(grid(i, j + 1, k) - grid(i, j - 1, k))/dr(2)
510 dgrid(i, j, k, 3) = 0.5*(grid(i, j, k + 1) - grid(i, j, k - 1))/dr(3)
514 END SUBROUTINE d_finite_difference
525 PURE FUNCTION trilinear_interpolation(r, subgrid, origin, dr)
RESULT(value_at_r)
526 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: r
527 REAL(kind=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: subgrid
528 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: origin, dr
529 REAL(kind=
dp) :: value_at_r
531 REAL(kind=
dp),
DIMENSION(3) :: norm_r, norm_r_rev
533 norm_r = (r - origin)/dr
534 norm_r_rev = 1 - norm_r
535 value_at_r = subgrid(1, 1, 1)*product(norm_r_rev) + &
536 subgrid(2, 1, 1)*norm_r(1)*norm_r_rev(2)*norm_r_rev(3) + &
537 subgrid(1, 2, 1)*norm_r_rev(1)*norm_r(2)*norm_r_rev(3) + &
538 subgrid(1, 1, 2)*norm_r_rev(1)*norm_r_rev(2)*norm_r(3) + &
539 subgrid(1, 2, 2)*norm_r_rev(1)*norm_r(2)*norm_r(3) + &
540 subgrid(2, 1, 2)*norm_r(1)*norm_r_rev(2)*norm_r(3) + &
541 subgrid(2, 2, 1)*norm_r(1)*norm_r(2)*norm_r_rev(3) + &
542 subgrid(2, 2, 2)*product(norm_r)
543 END FUNCTION trilinear_interpolation
553 ELEMENTAL FUNCTION pbc_index(i, lowbound, upbound)
554 INTEGER,
INTENT(IN) :: i, lowbound, upbound
557 IF (i < lowbound)
THEN
558 pbc_index = upbound + i - lowbound + 1
559 ELSE IF (i > upbound)
THEN
560 pbc_index = lowbound + i - upbound - 1
564 END FUNCTION pbc_index
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_cube_to_pw(grid, filename, scaling, silent)
Thin wrapper around routine cube_to_pw.
subroutine, public get_generic_info(gen_section, func_name, xfunction, parameters, values, var_values, size_variables, i_rep_sec, input_variables)
Reads from the input structure all information for generic functions.
This public domain function parser module is intended for applications where a set of mathematical ex...
real(rn) function, public evalf(i, Val)
...
real(kind=rn) function, public evalfd(id_fun, ipar, vals, h, err)
Evaluates derivatives.
subroutine, public finalizef()
...
subroutine, public initf(n)
...
subroutine, public parsef(i, FuncStr, Var)
Parse ith function string FuncStr and compile it into bytecode.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Interface to Maxwell equation solver.
subroutine, public maxwell_solver(maxwell_control, v_ee, sim_step, sim_time, scaling_factor)
Computes the external potential on the grid.
Interface to the message passing library MPI.
Define the data structure for the particle information.
integer, parameter, public pw_mode_local
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Routines to handle an external electrostatic field The external field can be generic and is provided ...
subroutine, public external_c_potential(qs_env, calculate_forces)
Computes the force and the energy due to the external potential on the cores.
subroutine, public external_e_potential(qs_env)
Computes the external potential on the grid.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.