76#include "./base/base_uses.f90"
82 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_dispersion_pairpot'
100 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
105 CHARACTER(len=*),
PARAMETER :: routinen =
'qs_dispersion_pairpot_init'
107 CHARACTER(LEN=2) :: symbol
108 CHARACTER(LEN=default_path_length) :: filename
109 CHARACTER(LEN=default_string_length) :: aname
110 CHARACTER(LEN=default_string_length), &
111 DIMENSION(:),
POINTER :: tmpstringlist
112 INTEGER :: elem, handle, i, ikind, j, max_elem, &
113 maxc, n_rep, nkind, nl, vdw_pp_type, &
115 INTEGER,
DIMENSION(:),
POINTER :: exlist
116 LOGICAL :: at_end, explicit, found, is_available
121 CALL timeset(routinen, handle)
123 nkind =
SIZE(atomic_kind_set)
125 vdw_type = dispersion_env%type
126 SELECT CASE (vdw_type)
131 vdw_pp_type = dispersion_env%type
132 SELECT CASE (dispersion_env%pp_type)
138 CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
142 filename = dispersion_env%parameter_file_name
145 IF (
PRESENT(pp_section))
THEN
149 c_vals=tmpstringlist)
150 IF (trim(tmpstringlist(1)) == trim(symbol))
THEN
152 READ (tmpstringlist(2), *) disp%c6
153 READ (tmpstringlist(3), *) disp%vdw_radii
159 IF (.NOT. found)
THEN
161 CALL dftd2_param(elem, disp%c6, disp%vdw_radii, found)
163 IF (.NOT. found)
THEN
165 INQUIRE (file=filename, exist=is_available)
166 IF (is_available)
THEN
175 IF (trim(aname) == trim(symbol))
THEN
180 disp%vdw_radii = disp%vdw_radii*
bohr
190 disp%defined = .true.
192 disp%defined = .false.
195 IF (.NOT. disp%defined) &
196 CALL cp_abort(__location__, &
197 "Dispersion parameters for element ("//trim(symbol)//
") are not defined! "// &
198 "Please provide a valid set of parameters through the input section or "// &
199 "through an external file! ")
200 CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
209 dispersion_env%max_elem = max_elem
210 dispersion_env%maxc = maxc
211 ALLOCATE (dispersion_env%maxci(max_elem))
212 ALLOCATE (dispersion_env%c6ab(max_elem, max_elem, maxc, maxc, 3))
213 ALLOCATE (dispersion_env%r0ab(max_elem, max_elem))
214 ALLOCATE (dispersion_env%rcov(max_elem))
215 ALLOCATE (dispersion_env%eneg(max_elem))
216 ALLOCATE (dispersion_env%r2r4(max_elem))
217 ALLOCATE (dispersion_env%cn(max_elem))
220 filename = dispersion_env%parameter_file_name
221 CALL dftd3_c6_param(dispersion_env%c6ab, dispersion_env%maxci, filename, para_env)
222 CALL setr0ab(dispersion_env%r0ab, dispersion_env%rcov, dispersion_env%r2r4)
224 CALL seten(dispersion_env%eneg)
226 CALL setcn(dispersion_env%cn)
234 dum = 0.5_dp*dispersion_env%r2r4(i)*real(i,
dp)**0.5_dp
236 dispersion_env%r2r4(i) = sqrt(dum)
239 dispersion_env%k1 = 16.0_dp
240 dispersion_env%k2 = 4._dp/3._dp
248 dispersion_env%k3 = -4._dp
249 dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*
bohr
251 dispersion_env%alp = 14._dp
254 CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
258 disp%defined = .true.
260 disp%defined = .false.
262 IF (.NOT. disp%defined) &
263 CALL cp_abort(__location__, &
264 "Dispersion parameters for element ("//trim(symbol)//
") are not defined! "// &
265 "Please provide a valid set of parameters through the input section or "// &
266 "through an external file! ")
267 CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
270 IF (
PRESENT(pp_section))
THEN
274 ALLOCATE (dispersion_env%cnkind(n_rep))
277 c_vals=tmpstringlist)
278 READ (tmpstringlist(1), *) dispersion_env%cnkind(i)%cnum
279 READ (tmpstringlist(2), *) dispersion_env%cnkind(i)%kind
284 ALLOCATE (dispersion_env%cnlist(n_rep))
287 c_vals=tmpstringlist)
288 nl =
SIZE(tmpstringlist)
289 ALLOCATE (dispersion_env%cnlist(i)%atom(nl - 1))
290 dispersion_env%cnlist(i)%natom = nl - 1
291 READ (tmpstringlist(1), *) dispersion_env%cnlist(i)%cnum
293 READ (tmpstringlist(j + 1), *) dispersion_env%cnlist(i)%atom(j)
301 DO j = 1,
SIZE(exlist)
303 CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp)
304 disp%defined = .false.
308 dispersion_env%nd3_exclude_pair = n_rep
310 ALLOCATE (dispersion_env%d3_exclude_pair(n_rep, 2))
314 dispersion_env%d3_exclude_pair(i, :) = exlist
324 CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
327 disp%defined = .true.
328 CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
333 dispersion_env%max_elem = max_elem
334 dispersion_env%maxc = maxc
335 ALLOCATE (dispersion_env%maxci(max_elem))
336 ALLOCATE (dispersion_env%rcov(max_elem))
337 ALLOCATE (dispersion_env%eneg(max_elem))
338 ALLOCATE (dispersion_env%cn(max_elem))
340 CALL setrcov(dispersion_env%rcov)
342 CALL setcn(dispersion_env%cn)
344 CALL seten(dispersion_env%eneg)
346 dispersion_env%k1 = 16.0_dp
347 dispersion_env%k2 = 4._dp/3._dp
348 dispersion_env%k3 = -4._dp
349 dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*
bohr
350 dispersion_env%alp = 14._dp
352 dispersion_env%cnfun = 3
354 IF (
PRESENT(pp_section))
THEN
361 CALL timestop(handle)
375 TYPE(qs_environment_type),
POINTER :: qs_env
376 TYPE(qs_dispersion_type),
POINTER :: dispersion_env
377 REAL(kind=dp),
INTENT(INOUT) :: energy
378 LOGICAL,
INTENT(IN) :: calculate_forces
379 REAL(kind=dp),
DIMENSION(:),
OPTIONAL :: atevdw
381 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_dispersion_pairpot'
383 INTEGER :: atom_a, handle, iatom, ikind, iw, natom, &
385 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
386 LOGICAL :: atenergy, atex, debugall, use_virial
387 REAL(kind=dp) :: evdw, gnorm
388 REAL(kind=dp),
ALLOCATABLE,
DIMENSION(:) :: atomic_energy
389 REAL(kind=dp),
DIMENSION(3) :: fdij
390 REAL(kind=dp),
DIMENSION(3, 3) :: dvirial, pv_loc, pv_virial_thread
391 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
392 TYPE(atprop_type),
POINTER :: atprop
393 TYPE(cell_type),
POINTER :: cell
394 TYPE(cp_logger_type),
POINTER :: logger
395 TYPE(mp_para_env_type),
POINTER :: para_env
396 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
397 TYPE(virial_type),
POINTER :: virial
403 IF (dispersion_env%type .NE. xc_vdw_fun_pairpot)
THEN
407 CALL timeset(routinen, handle)
409 NULLIFY (atomic_kind_set)
411 CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
412 cell=cell, virial=virial, para_env=para_env, atprop=atprop)
414 debugall = dispersion_env%verbose
417 logger => cp_get_default_logger()
418 IF (
ASSOCIATED(dispersion_env%dftd_section))
THEN
419 unit_nr = cp_print_key_unit_nr(logger, dispersion_env%dftd_section,
"PRINT_DFTD", &
426 atenergy = atprop%energy
429 IF (
PRESENT(atevdw))
THEN
433 IF (unit_nr > 0)
THEN
435 WRITE (unit_nr, *)
" Pair potential vdW calculation"
436 IF (dispersion_env%pp_type == vdw_pairpot_dftd2)
THEN
437 WRITE (unit_nr, *)
" Dispersion potential type: DFT-D2"
438 WRITE (unit_nr, *)
" Scaling parameter (s6) ", dispersion_env%scaling
439 WRITE (unit_nr, *)
" Exponential prefactor ", dispersion_env%exp_pre
440 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3)
THEN
441 WRITE (unit_nr, *)
" Dispersion potential type: DFT-D3"
442 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj)
THEN
443 WRITE (unit_nr, *)
" Dispersion potential type: DFT-D3(BJ)"
444 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd4)
THEN
445 WRITE (unit_nr, *)
" Dispersion potential type: DFT-D4"
449 CALL get_qs_env(qs_env=qs_env, force=force)
450 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
451 IF (use_virial .AND. debugall)
THEN
452 dvirial = virial%pv_virial
455 pv_loc = virial%pv_virial
459 pv_virial_thread(:, :) = 0._dp
461 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
463 IF (dispersion_env%pp_type == vdw_pairpot_dftd2)
THEN
464 CALL calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
465 ELSEIF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
466 dispersion_env%pp_type == vdw_pairpot_dftd3bj)
THEN
467 CALL calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
469 ELSEIF (dispersion_env%pp_type == vdw_pairpot_dftd4)
THEN
470 IF (dispersion_env%lrc)
THEN
471 cpabort(
"Long range correction with DFTD4 not implemented")
473 IF (dispersion_env%srb)
THEN
474 cpabort(
"Short range bond correction with DFTD4 not implemented")
476 IF (dispersion_env%domol)
THEN
477 cpabort(
"Molecular approximation with DFTD4 not implemented")
481 IF (dispersion_env%verbose) iw = cp_logger_get_default_io_unit(logger)
483 IF (atenergy .OR. atex)
THEN
484 ALLOCATE (atomic_energy(natom))
485 CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
486 iw, atomic_energy=atomic_energy)
488 CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw)
492 atevdw(1:natom) = atomic_energy(1:natom)
495 CALL atprop_array_init(atprop%atevdw, natom)
496 atprop%atevdw(1:natom) = atomic_energy(1:natom)
498 IF (atenergy .OR. atex)
THEN
499 DEALLOCATE (atomic_energy)
504 CALL para_env%sum(evdw)
506 IF (unit_nr > 0)
THEN
507 WRITE (unit_nr, *)
" Total vdW energy [au] :", evdw
508 WRITE (unit_nr, *)
" Total vdW energy [kcal] :", evdw*kcalmol
511 IF (calculate_forces .AND. debugall)
THEN
512 IF (unit_nr > 0)
THEN
513 WRITE (unit_nr, *)
" Dispersion Forces "
514 WRITE (unit_nr, *)
" Atom Kind Forces "
518 ikind = kind_of(iatom)
519 atom_a = atom_of_kind(iatom)
520 fdij(1:3) = force(ikind)%dispersion(:, atom_a)
521 CALL para_env%sum(fdij)
522 gnorm = gnorm + sum(abs(fdij))
523 IF (unit_nr > 0)
WRITE (unit_nr,
"(i5,i7,3F20.14)") iatom, ikind, fdij
525 IF (unit_nr > 0)
THEN
527 WRITE (unit_nr, *)
"|G| = ", gnorm
531 dvirial = virial%pv_virial - dvirial
532 CALL para_env%sum(dvirial)
533 IF (unit_nr > 0)
THEN
534 WRITE (unit_nr, *)
"Stress Tensor (dispersion)"
535 WRITE (unit_nr,
"(3G20.12)") dvirial
536 WRITE (unit_nr, *)
" Tr(P)/3 : ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
542 IF (calculate_forces .AND. use_virial)
THEN
543 virial%pv_vdw = virial%pv_vdw + (virial%pv_virial - pv_loc)
546 IF (
ASSOCIATED(dispersion_env%dftd_section))
THEN
547 CALL cp_print_key_finished_output(unit_nr, logger, dispersion_env%dftd_section,
"PRINT_DFTD")
550 CALL timestop(handle)
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
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.
Holds information on atomic properties.
subroutine, public atprop_array_init(atarray, natom)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public caldeweyher2020
integer, save, public grimme2006
integer, save, public caldeweyher2019
integer, save, public caldeweyher2017
integer, save, public goerigk2017
integer, save, public grimme2011
integer, save, public grimme2010
Handles all functions related to the CELL.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Interface to the message passing library MPI.
Definition of physical constants:
real(kind=dp), parameter, public kcalmol
real(kind=dp), parameter, public kjmol
real(kind=dp), parameter, public bohr
Coordination number routines for dispersion pairpotentials.
real(kind=dp) function, public get_cn_radius(dispersion_env)
...
subroutine, public setr0ab(rout, rcov, r2r4)
...
subroutine, public setrcov(rcov)
...
subroutine, public seten(enout)
...
subroutine, public setcn(cnout)
...
Calculation of D2 dispersion.
subroutine, public calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
...
subroutine, public dftd2_param(z, c6, r, found)
...
Calculation of D3 dispersion.
subroutine, public dftd3_c6_param(c6ab, maxci, filename, para_env)
...
subroutine, public calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, unit_nr, atevdw)
...
Calculation of dispersion using pair potentials.
subroutine, public calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, atomic_energy)
...
Calculation of dispersion using pair potentials.
subroutine, public qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
...
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition of disperson types for DFT calculations.
integer, parameter, public dftd2_pp
integer, parameter, public dftd4_pp
integer, parameter, public dftd3_pp
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_pp, 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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
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, zatom, 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_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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
Provides all information about an atomic kind.
type for the atomic properties
Type defining parameters related to the simulation cell.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.